real_multiplier[m_i] = gainScalar_re;
                        imag_multiplier[m_i] = gainScalar_im;
                     }
                     /* gains aren't scalar */
                     else {
                        ksnm = ksn*m_i+k*s*n_i+k*s_i+ k_i;
                        kunm = kun*m_i+k*u*n_i+k*u_i+k_i;
                        a = *(re_sq_G_BS + ksnm);
                        c = *(re_sq_G_MS + kunm);
                        if(both_real)
                           real_multiplier[m_i] = a*c;
                        else {
                           if (!bsreal) 
                              b = *(im_sq_G_BS + ksnm);
                           if (!msreal) 
                              d = *(im_sq_G_MS + kunm);
                           complex_multiply(a, b, c, d,
                              &real_multiplier[m_i], &imag_multiplier[m_i]);
                        }
                     }
                     exp_helper[m_i] = kds * cos_table[abs((*(aod+knm)-halfpi)*r2p)&divider]
                        + *(phase+knm) + kdu * cos_table[abs((*(aoa+knm)-halfpi)*r2p)&divider];

                     /* next value should be calculated accurately, because it is the main source of
                     * error when multiplied with large time sample values!
                     */
                     exp_t_coeff[m_i] =  kv * cos(*(aoa+knm)-theta_v[k_i]);	

                     if (tn >1) /* output pphase can be evaluated only if more than one time sample */
                        *(output_SubPathPhase+knm) = exp_t_coeff[m_i]* (*(ts+k*(tn-1)+k_i)+delta_t);
                  }  
                  
                  m_count = 0;
                  for (l_i=0;l_i<l;l_i++ ) {
                     knl = k*nl_i+k_i;
                     usnltk_i = usnl*tn*k_i + u*s*nl_i + u*s_i + u_i;
                     /* cycling the time vector */
                     for (t_i=0; t_i<tn; t_i++) {
                        real_sum = 0;
                        imag_sum = 0;
                        kt = k*t_i + k_i;
                        /* cycling m, now no calculations with m dependable variables is required */
                        m_i = m_count;
                        sumloop_with_table(m_i,m_count+(int)ln[l_i],exp_helper,exp_t_coeff, real_multiplier, imag_multiplier,
                           *(ts+kt), lm, &real_sum, &imag_sum, cos_table, r2p,divider, halfpi);
                        /* setting the output values and multiplying them with proper coefficient*/
                        usnltk = usnltk_i + usnl*t_i;
                        *(re_h+usnltk) = real_sum * *(sq_Pn+knl) * one_over_sq_ln[l_i];
                        *(im_h+usnltk) = imag_sum* *(sq_Pn+knl)* one_over_sq_ln[l_i];
                     }
                     m_count = m_count+(int)ln[l_i];
                     nl_i++;
                  }
               }
            }
         }
      }
      /* freeing memory of temporary tables */
      free(cos_table);
   }

   /* sin_look_up not used */
   else {

      /* cycling links */
      for(k_i=0; k_i<k; k_i++) {
         if(tn>1)
            delta_t = *(ts+k+k_i)-*(ts+k_i);
         kv = k_CONST * v[k_i];
         /* cycling u, u is a given constant*/
         for(u_i=0; u_i < u; u_i++) {
            kdu = k_CONST * *(d_u+u_i);
            /* cyckling s, s is a given constant */
            for(s_i=0;s_i< s; s_i++) {
               kds = k_CONST * *(d_s+s_i);
               /* running through paths, n is a given constant*/
               nl_i = 0;
               for(n_i=0;n_i<n;n_i++) {
                  kn = k*n_i +k_i;
                  /* calculating m-variant t-invariant values */
                  for(m_i=0; m_i<m; m_i++) {                     
                     /* calculation of pointer parameters */
                     knm = k*n*m_i+kn;

                     /* calculation of sqrt(G_BS)*sqrt(G_MS): */
                     if (GainsAreScalar) {
                        real_multiplier[m_i] = gainScalar_re;
                        imag_multiplier[m_i] = gainScalar_im;
                     }
                     /* gains aren't scalar */
                     else {
                        ksnm = ksn*m_i+k*s*n_i+k*s_i+ k_i;
                        kunm = kun*m_i+k*u*n_i+k*u_i+k_i;
                        a = *(re_sq_G_BS + ksnm);
                        c = *(re_sq_G_MS + kunm);
                        if(both_real)
                           real_multiplier[m_i] = a*c;
                        else {
                           if (!bsreal) 
                              b = *(im_sq_G_BS + ksnm);
                           if (!msreal) 
                              d = *(im_sq_G_MS + kunm);
                           complex_multiply(a, b, c, d,
                              &real_multiplier[m_i], &imag_multiplier[m_i]);
                        }
                     }

                     /* calculation of t-invariant part of exp(j...) in [1] 5.4 */
                     exp_helper[m_i] = kds * sin(*(aod+knm)) + *(phase+knm) +
                        kdu * sin(*(aoa+knm));

                     /* the coefficient just in front of t in exp(j...)*/
                     exp_t_coeff[m_i] = kv * cos(*(aoa+knm)-theta_v[k_i]);

                     if (tn >1) /* output pphase can be evaluated only if more than one time sample */
                     *(output_SubPathPhase+knm) = exp_t_coeff[m_i]* (*(ts+k*(tn-1)+k_i)+delta_t);
                  }

                  m_count = 0;
                  for (l_i=0;l_i<l;l_i++ ) {
                     knl = k*nl_i+k_i;
                     usnltk_i = usnl*tn*k_i + u*s*nl_i + u*s_i + u_i;
                     /* cycling the time vector */
                     for (t_i=0; t_i<tn; t_i++) {
                        real_sum = 0;
                        imag_sum = 0;
                        kt = k*t_i + k_i;

                        /* cycling m, now no calculations with m dependable variables is required */
                        m_i = m_count;
                        sumloop(m_i,m_count+(int)ln[l_i],exp_helper,exp_t_coeff, real_multiplier, imag_multiplier,
                           *(ts+kt), lm, &real_sum, &imag_sum);

                        /* setting the output values and multiplying them with proper coefficient*/
                        usnltk = usnltk_i + usnl*t_i;
                        *(re_h+usnltk) = real_sum * *(sq_Pn+knl) * one_over_sq_ln[l_i];
                        *(im_h+usnltk) = imag_sum* *(sq_Pn+knl)* one_over_sq_ln[l_i];
                     }
                     m_count = m_count+(int)ln[l_i];
                     nl_i++;
                  }
               }
            }
         }
      }
   }
   /* freeing memory of temporary tables */
   free(real_multiplier);
   free(imag_multiplier);
   free(exp_t_coeff);
   free(exp_helper);
   free(one_over_sq_ln);
   
   return 0;
}

/**
 * Function scm_pol_sum Creates array for Spatial channel model sums with cross-polarization.
 * The parameter angles must be in radians!
 * @param sin_look_up_points if nonzero, then look-up table for sin/cos is used. Use -1 for default number of points.
 * @param u the number of MS elements
 * @param s the number of BS elements
 * @param n number of multipaths
 * @param l number of midpaths
 * @param m number of subpaths for each n multipath
 * @param k the number of individual links
 * @param *re_X_BS_v real part of BS antenna V-pol component response (size [k][s][n][m])
 * @param *im_X_BS_v imaginary part of BS antenna V-pol component response (size [k][s][n][m])
 * @param *re_X_BS_h real part of BS antenna H-pol component response (size [k][s][n][m])
 * @param *im_X_BS_h imaginary part of BS antenna H-pol component response (size [k][s][n][m])
 * @param *re_X_MS_v real part of MS antenna V-pol component response (size [k][u][n][m])
 * @param *im_X_MS_v imaginary part of MS antenna V-pol component response (size [k][u][n][m])
 * @param *re_X_MS_h real part of MS antenna H-pol component response (size [k][u][n][m])
 * @param *im_X_MS_h imaginary part of MS antenna H-pol component response (size [k][u][n][m])
 * @param k wavenumber (2*pi/lambda)
 * @param *d_s distance of BS antenna s from ref. antenna (s=1), (size [s])
 * @param *d_u distance of MS antenna u from ref. antenna (u=1), (size [u])
 * @param *aod Angel of Departure (size [k][n][m])
 * @param *aoa Angel of Arrival (size [k][n][m])
 * @param *phase_v_v Phase offset of the mth subpath of the nth path between vertical components of BS and MS (size = [k][n][m])
 * @param *phase_v_h Phase offset of the mth subpath of the nth path between vertical components of BS and horizontal components of MS (size = [k][n][m])
 * @param *phase_h_v Phase offset of the mth subpath of the nth path between horizontal components of BS and vertical components of MS (size = [k][n][m])
 * @param *phase_h_h Phase offset of the mth subpath of the nth path between horizontal components of BS and MS (size = [k][n][m])
 * @param *sq_r_n1 square root of power ratio of (v-h)/(v-v) (size [k][n])
 * @param *sq_r_n2 square root of power ratio of (h-v)/(v-v) (size [k][n])
 * @param *v magnitude of the MS velocity vector (size [k])
 * @param *theta_v angle of the MS velocity vector (size [k])
 * @param *ts vectors of time samples (size [k][tn])
 * @param tn number of time samples
 * @param *sq_Pn square root of Pn (size [k][n*l])
 * @param *ln number of subpaths per midpath (size [l])
 * @param *lm subpath indexing order for midpaths, Matlab indexing: 1-m (size [m])
 * @param GainsAreScalar has value 1 if gain is scalar, zero if not
 * @param *re_h pointer to real part output matrice h, must be initialized outside scm function (size [u][s][n][t][k])
 * @param *im_h pointer to imag part output matrice h, must be initialized outside scm function (size [u][s][n][t][k])
 * @param *output_SubPathPhase (size [k][n][m])
 * @return 0 if success, 1 if bad arguments
 */
int scm_pol_sum(const long int sin_look_up_points, const int u, const int s, const int n, const int l, const int m, const int k,
   const double *re_X_BS_v, const double *im_X_BS_v, const double *re_X_BS_h, const double *im_X_BS_h, 
   const double *re_X_MS_v, const double *im_X_MS_v, const double *re_X_MS_h, const double *im_X_MS_h, 
   const double k_CONST, const double *d_s, const double *d_u, const double *aod, const double *aoa,
   const double *phase_v_v, const double *phase_v_h, const double *phase_h_v, const double *phase_h_h,
   const double *r_n1, const double *r_n2,	const double *v, const double *theta_v, const double *ts,
   const int tn, const double *sq_Pn, const double *ln, const double *lm, int GainsAreScalar,
   double *re_h, double *im_h, double *output_SubPathPhase) {

   int i, n_i, nl_i, l_i, m_i, u_i, s_i, t_i, k_i, m_count;	/* running index variables */
   /* notation conversion variables from [][] to *(p[0]+var) */
   long int ksnm, kunm, knm, usnltk, kn, knl, nl, kt, usnl, usnltk_i, ksn, kun;

   int bsreal, msreal, both_real;
   double a1_re, a1_im, a2_re, a2_im, c1_re, c1_im, c2_re, c2_im; /* variables for multiplications*/
   double b11_re, b11_im, b12_re, b12_im, b21_re, b21_im, b22_re, b22_im;
   double a1c1_re, a1c2_re, a2c1_re, a2c2_re;
   double tmp; /* temporary variable to lighten calculations */
   double kds, kdu, kv, real_sum, imag_sum;
   double *real_multiplier;
   double *imag_multiplier;
   double *exp_t_coeff, *exp_helper;
   double delta_t;
   double *cos_table;
   double r2p;
   double halfpi;
   long int num_of_points, divider, pows_of_two;
   double *one_over_sq_ln;
   
   double sq_r_n1, sq_r_n2;     /* (Jari, Apr 17, 2005) */
   
   real_multiplier = (double*)malloc(m*sizeof(double));
   imag_multiplier = (double*)calloc(m,sizeof(double));
   exp_t_coeff = (double*)malloc(m*sizeof(double));
   exp_helper = (double*)malloc(m*sizeof(double));

   /* check if input gains are not complex */
   bsreal = 0;
   msreal = 0;
   both_real = 0;
   if (im_X_BS_v == NULL)
      bsreal = 1;
   if (im_X_MS_v == NULL)
      msreal = 1;
   if(bsreal && msreal)
      both_real = 1;

   if (GainsAreScalar) {
      if (both_real) {
         a1c1_re = *re_X_BS_v * *re_X_MS_v;
         a2c1_re = *re_X_BS_h * *re_X_MS_v;
         a1c2_re = *re_X_BS_v * *re_X_MS_h;
         a2c2_re = *re_X_BS_h * *re_X_MS_h;
      }
   }

   nl = n*l;
   usnl = u*s*nl;
   ksn = k*s*n;
   kun = k*u*n;

   /* calculate coefficient terms */
   one_over_sq_ln = (double*)malloc(l*sizeof(double));
   for(i=0;i<l;i++) {
      one_over_sq_ln[i] = 1/sqrt(ln[i]); 
   }

   /* checks if look-up table is used for sin/cos */
   if (sin_look_up_points) {
      if (sin_look_up_points == -1)
         num_of_points = DEFAULT_LUT_POINTS;
      else {
         pows_of_two = 2;
         while (pows_of_two < sin_look_up_points)
            pows_of_two = pows_of_two*2;
         num_of_points = pows_of_two;
         if (pows_of_two != sin_look_up_points)
            printf("Warning: Number of LU points is not a power-of-2: size changed to %li.\n", num_of_points);
      }

      divider =  num_of_points -1;
      r2p = num_of_points/(2*PI);
      halfpi = PI*0.5;

      cos_table = (double*)malloc(num_of_points*sizeof(double));
      /* calculate the table */
      for (i=0;i<num_of_points;i++) 
         cos_table[i] = cos((i+0.5)/r2p);

      /* cycling links */
      for(k_i=0; k_i<k; k_i++) {
         if(tn>1)
            delta_t = *(ts+k+k_i)-*(ts+k_i);
         kv = k_CONST * v[k_i]; /* value depends only on k */
         /* u is a given constant */
         for (u_i=0; u_i<u; u_i++) {
            kdu = k_CONST * *(d_u+u_i); /* value depends only on u */
            /* s is a given constant */
            for (s_i=0; s_i<s; s_i++) {
               kds = k_CONST * *(d_s+s_i); /* value depends only on s */
               /* running through paths, n is a given constant */
               nl_i=0;
               for (n_i=0; n_i<n; n_i++) {
                  kn = k*n_i +k_i;
                  /* m is a given constant */
                  for (m_i=0; m_i<m; m_i++) {
                     /* calculation of pointer parameters */
                     knm = k*n*m_i+kn;
                     
                     /* Jari (Apr 17, 2005): changed power normalization */
                     sq_r_n1 = sqrt(*(r_n1+kn));
                     sq_r_n2 = sqrt(*(r_n2+kn));
                     
                     b12_re = sq_r_n1 * cos_table[abs((*(phase_v_h+knm))*r2p)&divider];
                     b12_im = sq_r_n1 * cos_table[abs((*(phase_v_h+knm)-halfpi)*r2p)&divider];
                     b11_re = cos_table[abs((*(phase_v_v+knm))*r2p)&divider];
                     b11_im = cos_table[abs((*(phase_v_v+knm)-halfpi)*r2p)&divider];                        
                     b21_re = sq_r_n2 * cos_table[abs((*(phase_h_v+knm))*r2p)&divider];
                     b21_im = sq_r_n2 * cos_table[abs((*(phase_h_v+knm)-halfpi)*r2p)&divider];
                     b22_re = cos_table[abs((*(phase_h_h+knm))*r2p)&divider];
                     b22_im = cos_table[abs((*(phase_h_h+knm)-halfpi)*r2p)&divider];
                     


                     /* Jari (Apr 17, 2005): Xpd independent power, not used in this version */
                     /*
                     tmp = 1/sqrt(1+*(r_n1+kn));
                     b12_re = tmp * cos_table[abs((*(phase_v_h+knm))*r2p)&divider];
                     b12_im = tmp * cos_table[abs((*(phase_v_h+knm)-halfpi)*r2p)&divider];
                     tmp = sqrt(*(r_n1+kn))*tmp;
                     b11_re =  tmp * cos_table[abs((*(phase_v_v+knm))*r2p)&divider];
                     b11_im = tmp * cos_table[abs((*(phase_v_v+knm)-halfpi)*r2p)&divider];
                     tmp = 1/sqrt(1+*(r_n2+kn));
                     b21_re = tmp * cos_table[abs((*(phase_h_v+knm))*r2p)&divider];
                     b21_im = tmp * cos_table[abs((*(phase_h_v+knm)-halfpi)*r2p)&divider];
                     tmp = sqrt(*(r_n2+kn))*tmp;
                     b22_re = tmp * cos_table[abs((*(phase_h_h+knm))*r2p)&divider];
                     b22_im = tmp * cos_table[abs((*(phase_h_h+knm)-halfpi)*r2p)&divider];
                     */

                     if(both_real) { 
                        if (GainsAreScalar) {
                           real_multiplier[m_i] = a1c1_re *b11_re + a2c1_re * b21_re +