a1c2_re*b12_re + a2c2_re*b22_re;
                           imag_multiplier[m_i] = a1c1_re*b11_im + a2c1_re*b21_im +
                              a1c2_re*b12_im + a2c2_re*b22_im;
                        }
                        else { /* aren't scalar, but real*/
                           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;
                           a1_re = *(re_X_BS_v+ksnm);
                           a2_re = *(re_X_BS_h+ksnm);
                           c1_re = *(re_X_MS_v+kunm);
                           c2_re = *(re_X_MS_h+kunm);
                           real_multiplier[m_i] = a1_re*c1_re*b11_re + a2_re*c1_re*b21_re + 
                              a1_re*c2_re*b12_re + a2_re*c2_re*b22_re;
                           imag_multiplier[m_i] = a1_re*c1_re*b11_im + a2_re*c1_re*b21_im +
                              a1_re*c2_re*b12_im + a2_re*c2_re*b22_im;
                        }
                     }
                     else { /* aren't real*/
                        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;
                        a1_re = *(re_X_BS_v+ksnm);
                        a2_re = *(re_X_BS_h+ksnm);
                        c1_re = *(re_X_MS_v+kunm);
                        c2_re = *(re_X_MS_h+kunm);
                        if (bsreal) {
                           a1_im= 0.0;
                           a2_im = 0.0;
                        }
                        else {
                           a1_im = *(im_X_BS_v+ksnm);
                           a2_im = *(im_X_BS_h+ksnm);
                        }
                        if (msreal) {
                           c1_im = 0.0;
                           c2_im = 0.0;
                        }
                        else {
                           c1_im = *(im_X_MS_v+kunm);
                           c2_im =  *(im_X_MS_h+kunm);
                        }
                        matrix_multiply(a1_re,a1_im,a2_re,a2_im,c1_re,c1_im,c2_re,c2_im,
                           b11_re,b11_im,b12_re,b12_im,b21_re,b21_im,b22_re,b22_im,
                           &real_multiplier[m_i], &imag_multiplier[m_i]);
                     }

                     /* time invariant exponent term */
                     exp_helper[m_i] = kds * cos_table[abs((*(aod+knm)-halfpi)*r2p)&divider]
                        + kdu * cos_table[abs((*(aoa+knm)-halfpi)*r2p)&divider];

                     /* the coefficient just in front of t in exp(j...)
                     * 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++;
                  }
               }
            }
         }
      }
      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]; /* 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(*(phase_v_h+knm));
                     b12_im = sq_r_n1 * sin(*(phase_v_h+knm));
                     b11_re = cos(*(phase_v_v+knm));
                     b11_im = sin(*(phase_v_v+knm));
                     b21_re = sq_r_n2 * cos(*(phase_h_v+knm));
                     b21_im = sq_r_n2 * sin(*(phase_h_v+knm));
                     b22_re = cos(*(phase_h_h+knm));
                     b22_im = sin(*(phase_h_h+knm));

                     
                     
                     /* Jari (Apr 17, 2005): Xpd independent power normalization, not used in this version */
                     /*
                     tmp = 1/sqrt(1+*(r_n1+kn));
                     b12_re = tmp * cos(*(phase_v_h+knm));
                     b12_im = tmp * sin(*(phase_v_h+knm));
                     tmp = sqrt(*(r_n1+kn))*tmp;
                     b11_re =  tmp * cos(*(phase_v_v+knm));
                     b11_im = tmp * sin(*(phase_v_v+knm));
                     tmp = 1/sqrt(1+*(r_n2+kn));
                     b21_re = tmp * cos(*(phase_h_v+knm));
                     b21_im = tmp * sin(*(phase_h_v+knm));
                     tmp = sqrt(*(r_n2+kn))*tmp;
                     b22_re = tmp * cos(*(phase_h_h+knm));
                     b22_im = tmp * sin(*(phase_h_h+knm));
                     */

                     if(both_real) { 
                        if (GainsAreScalar) {
                           real_multiplier[m_i] = a1c1_re *b11_re + a2c1_re * b21_re +
                              a1c2_re*b12_re + a2c2_re*b22_re;
                           imag_multiplier[m_i] = a1c1_re*b11_im + a2c1_re*b21_im +
                              a1c2_re*b12_im + a2c2_re*b22_im;
                        }
                        else { /* aren't scalar*/
                           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;
                           a1_re = *(re_X_BS_v+ksnm);
                           a2_re = *(re_X_BS_h+ksnm);
                           c1_re = *(re_X_MS_v+kunm);
                           c2_re = *(re_X_MS_h+kunm);
                           real_multiplier[m_i] = a1_re*c1_re*b11_re + a2_re*c1_re*b21_re + 
                              a1_re*c2_re*b12_re + a2_re*c2_re*b22_re;
                           imag_multiplier[m_i] = a1_re*c1_re*b11_im + a2_re*c1_re*b21_im +
                              a1_re*c2_re*b12_im + a2_re*c2_re*b22_im;
                        }
                     }
                     else { /* aren't real*/
                        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;
                        a1_re = *(re_X_BS_v+ksnm);
                        a2_re = *(re_X_BS_h+ksnm);
                        c1_re = *(re_X_MS_v+kunm);
                        c2_re = *(re_X_MS_h+kunm);
                        if (bsreal) {
                           a1_im= 0.0;
                           a2_im = 0.0;
                        }
                        else {
                           a1_im = *(im_X_BS_v+ksnm);
                           a2_im = *(im_X_BS_h+ksnm);
                        }
                        if (msreal) {
                           c1_im = 0.0;
                           c2_im = 0.0;
                        }
                        else {
                           c1_im = *(im_X_MS_v+kunm);
                           c2_im =  *(im_X_MS_h+kunm);
                        }
                        matrix_multiply(a1_re,a1_im,a2_re,a2_im,c1_re,c1_im,c2_re,c2_im,
                           b11_re,b11_im,b12_re,b12_im,b21_re,b21_im,b22_re,b22_im,
                           &real_multiplier[m_i], &imag_multiplier[m_i]);
                     }

                     /* time invariant exponent term */
                     exp_helper[m_i] = kds * sin(*(aod+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_los creates a line of sight component.
 * The parameter angles must be in radians!
 * @param u the number of MS elements
 * @param s the number of BS elements
 * @param n number of paths
 * @param k the number of individual links
 * @param *re_G_BS real part of square root of BS Gain coefficients (size [k][s])
 * @param *im_G_BS imaginary part of BS Gain coefficients (size [k][s])
 * @param *re_G_MS real part of MS Gain coefficients (size [k][u])
 * @param *im_G_MS imaginary part of MS Gain coefficients (size [k][u])
 * @param k_CONST 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 *theta_BS Angel of Departure for LOS (size [k])
 * @param *theta_MS Angel of Arrival for LOS (size [k])
 * @param phase_los phase of LOS component (size [k])
 * @param *v magnitude of the MS velocity vector (size [k])
 * @param *theta_v angle of the MS velocity vector (size [k])
 * @param *ts vector of time samples (size [k][tn])
 * @param tn number of time samples
 * @param *k_FACTOR the Ricean K factors (size [k])
 * @return 0 if success, 1 if bad arguments
 * @param *re_h pointer to real part of output los (size[u][s][t])
 * @param *im_h pointer to imag part output matrice h, re_h has to be initialized outside scm function
 * @param *output_los_phase output LOS phase (size [k])

 * @return 0 if success, 1 if bad arguments
 */
int scm_los(const int u, const int s, const int n, const int k, const double *re_G_BS, const double *im_G_BS,
   const double *re_G_MS, const double *im_G_MS,	const double k_CONST,
   const double *d_s, const double *d_u, const double *theta_BS, const double *theta_MS,
   const double *phase_los, const double *v, const double *theta_v, const double *ts, const int tn,
   double *k_FACTOR, double *re_h, double *im_h, double *output_los_phase) {

   int t_i, u_i, s_i, k_i, n_i;	/* running index variables */
   long int usntk_i, usn, usntk, ks, ku, kt; 	/* notation conversion variables from [][] to *(p[0]+var) */
   double delta_t;
   double a, b, angle;
   double real_multiplier, imag_multiplier, re_LOS, im_LOS;
   double exp_helper; /* part of coefficient for exp(j...) in the sum */
   double exp_t_coeff;
   double sq_1_over_K_plus_1, sq_K;


   for(k_i=0;k_i<k;k_i++) {
      /* check for LOS component */
      if (k_FACTOR[k_i] != 0) {
         
         /* the coefficients for terms */
         sq_1_over_K_plus_1 = sqrt(1/(k_FACTOR[k_i]+1));
         sq_K = sqrt(k_FACTOR[k_i]);

         /* the coefficient just in front of t in exp(j...)
         doesn't depend on u or s */
         exp_t_coeff = k_CONST * v[k_i] * cos(theta_MS[k_i]-theta_v[k_i]);

         /* output phase can be evaluated only if more than one time sample */
         if (tn >1) {
            delta_t = *(ts+k*1+k_i)-*(ts+k_i);
            *(output_los_phase+k_i) = exp_t_coeff * (*(ts+k*(tn-1)+k_i)+delta_t);
         }
         usn = u*s*n;

         /* u is a given constant*/
         for(u_i=0; u_i < u; u_i++) {
            /* s is a given constant */
            for(s_i=0;s_i< s; s_i++) {

               for(n_i=0; n_i<n; n_i++) {
                  
                  usntk_i = usn*tn*k_i + u*s*n_i + u*s_i + u_i;
                  /* for all components */
                  for(t_i=0; t_i<tn; t_i++) {
                     usntk = usntk_i + usn*t_i;
                     *(re_h + usntk) = sq_1_over_K_plus_1 * *(re_h + usntk);
                     *(im_h + usntk) = sq_1_over_K_plus_1 * *(im_h + usntk);
                  }

                  /* for LOS components */
                  if (n_i == 0) {

                     ks = k*s_i + k_i;
                     ku = k*u_i + k_i;

                     /* Multiplying gain components (t-invariant) */
                     /* check G_BS and G_MS for being real valued */
                     if (im_G_BS != NULL || im_G_MS != NULL) {
                        /* Both are imaginary */
                        if (im_G_BS != NULL && im_G_MS != NULL) {
                           complex_multiply(*(re_G_BS+ks), *(im_G_BS+ks),
                              *(re_G_MS+ku), *(im_G_MS+ku), &a, &b);
                        }
                        /* only G_BS is complex */
                        else if (im_G_BS != NULL && im_G_MS == NULL) {
                           complex_multiply(*(re_G_BS+ks), *(im_G_BS+ks),
                              *(re_G_MS+ku), 0.0, &a, &b);
                        }
                        /* only G_MS is complex */
                        else {
                           complex_multiply(*(re_G_BS+ks), 0.0,
                              *(re_G_MS+ku), *(im_G_MS+ku), &a, &b);