}
                     }
                     /* all is real */
                     else {
                        a = *(re_G_BS+ks) * *(re_G_MS+ku);
                        b = 0.0;
                     }

                     /* calculation of t-invariant part of exp(j...) in [1] 5.4 */
                     exp_helper = k_CONST * (*(d_s+s_i) * sin(theta_BS[k_i]) + 
                        *(d_u+u_i) * sin(theta_MS[k_i])) + phase_los[k_i];

                     /* product of previous two */
                     complex_multiply(a, b, cos(exp_helper), sin(exp_helper),
                        &real_multiplier, &imag_multiplier);

                     usntk_i = usn*tn*k_i + u*s*n_i + u*s_i + u_i;
                     for(t_i=0; t_i<tn; t_i++) {
                        usntk = usntk_i + usn*t_i;
                        kt = k*t_i + k_i;
                        angle = exp_t_coeff * *(ts+kt);
                        complex_multiply(real_multiplier, imag_multiplier,
                           cos(angle), sin(angle),
                           &re_LOS, &im_LOS);
                        *(re_h + usntk) += sq_K * sq_1_over_K_plus_1 * re_LOS;
                        *(im_h + usntk) += sq_K * sq_1_over_K_plus_1 * im_LOS;
                     }
                  }
               }
            }
         }
      }
   }
   return 0;
}



/**
 * This is the mex-function.
 *
 * Input arguments are:
 *
 * FOR GENERAL COEFFICIENTS MODE:
 *
 * prhs[0] = mode, either 1 (GENERAL), 2 (POLARIZED) or 3 (LOS) 
 * prhs[1] = [G_BS], complex BS_antenna gains (size [k][s][n][m])
 * prhs[2] = [G_MS], complex MS_antenna gains (size [k][u][n][m])
 * prhs[3] = [aod], angles of departure (size [k][n][m])
 * prhs[4] = [aoa], angles of arrival (size [k][n][m])
 * prhs[5] = [d_s], distances of BS antenna s from ref. antenna (s=1), (size [s])
 * prhs[6] = [d_u], distance of MS antenna u from ref. antenna (u=1), (size [u])
 * prhs[7] = [phase], phase of the mth subpath of the nth path (size [k][n][m])
 * prhs[8] = [ts], time samples (size [k][tn])
 * prhs[9] = k_CONST, wave number
 * prhs[10] = [v], magnitude of the MS velocity vector (size [k])
 * prhs[11] = [theta_v], angle of the MS velocity vector (size [k])
 * prhs[12] = [sq_Pn], square root of Pn (size [k][n*l])
 * prhs[13] = look_up_points, 0 if look-up table not used, -1 if default, otherwise the number of points
 * prhs[14] = u, number of MS antennas
 * prhs[15] = s, number of BS antennas
 * prhs[16] = n, number of multipaths
 * prhs[17] = l, number of midpaths
 * prhs[18] = m, number of subpaths
 * prhs[19] = k, number of links
 * prhs[20] = tn, number of time samples
 * prhs[21] = GainsAreScalar, nonzero if gains are scalar
 * prhs[22] = [lm], indexing vector for subpaths (size [m])
 * prhs[23] = [ln], number of subpaths per midpath (size [l])
 * Output arguments are 
 * plhs[0] = [h], output array of the channel coefficients (size [u][s][n][tn][k])
 * plhs[1] = [output_SubPathPhase] pointer to output phases (size [k][n][m])
 *
 * FOR POLARIZED MODE:
 *
 * prhs[0] = mode, either 1 (GENERAL), 2 (POLARIZED) or 3 (LOS) 
 * prhs[1] = [X_BS_v] BS antenna V-pol component response (size[k][s][n][m])
 * prhs[2] = [X_BS_h] BS antenna H-pol component response (size[k][s][n][m])
 * prhs[3] = [X_MS_v] MS antenna V-pol component response (size[k][u][n][m])
 * prhs[4] = [X_MS_h] MS antenna H-pol component response (size[k][u][n][m])
 * prhs[5] = [aod], angles of departure (size[k][n][m])
 * prhs[6] = [aoa], angles of arrival (size[k][n][m])
 * prhs[7] = [d_s], distances of BS antenna s from ref. antenna (s=1), (size [s])
 * prhs[8] = [d_u], distance of MS antenna u from ref. antenna (u=1), (size [u])
 * prhs[9] = [phase_v_v], Phase offset of the mth subpath of the nth path between vertical components of BS and MS (size = [k][n][m])
 * prhs[10] = [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])
 * prhs[11] = [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])
 * prhs[12] = [phase_h_h], Phase offset of the mth subpath of the nth path between horizontal components of BS and MS (size = [k][n][m])
 * prhs[13] = [r_n1] power ratio of (v-h)/(v-v) (size [k][n])
 * prhs[14] = [r_n2] power ratio of (h-v)/(v-v) (size [k][n])
 * prhs[15] = [ts], time sample vector (size [k][tn])
 * prhs[16] = k_CONST, wave number
 * prhs[17] = [v], magnitude of the MS velocity vector (size [k])
 * prhs[18] = [theta_v], angle of the MS velocity vector (size [k])
 * prhs[19] = [sq_Pn], square root of Pn (size [k][n*l])
 * prhs[20] = look_up_points, 0 if look-up table not used, -1 if default, otherwise the number of points
 * prhs[21] = u, number of MS antennas
 * prhs[22] = s, number of BS antennas
 * prhs[23] = n, number of multipaths
 * prhs[24] = l, number of midpaths
 * prhs[25] = m, number of subpaths
 * prhs[26] = k, number of links
 * prhs[27] = tn, number of time samples
 * prhs[28] = GainsAreScalar, nonzero if gains are scalar
 * prhs[29] = [lm], indexing vector for subpaths (size [m])
 * prhs[30] = [ln], number of subpaths per midpath (size [l])
 * Output arguments are 
 * plhs[0] = [h], output array of the channel coefficients (size [u][s][n][tn])
 * plhs[1] = [output_SubPathPhase] pointer to output phases (size [n][m])
 *
 * FOR LOS MODE:
 *
 * prhs[0] = mode, either 1 (GENERAL), 2 (POLARIZED) or 3 (LOS) 
 * prhs[1] = [G_BS], complex BS_antenna gains (size [k][s])
 * prhs[2] = [G_MS], complex MS_antenna gains (size [k][u])
 * prhs[3] = theta_BS, angles of departure (size [k])
 * prhs[4] = theta_MS, angles of arrival (size [k])
 * prhs[5] = [d_s], distances of BS antenna s from ref. antenna (s=1), (size [s])
 * prhs[6] = [d_u], distance of MS antenna u from ref. antenna (u=1), (size [u])
 * prhs[7] = phase_los, phase of LOS component (size [k])
 * prhs[8] = [ts], time sample vector (size [k][tn])
 * prhs[9] = k_CONST, wave number
 * prhs[10] = [v], magnitude of the MS velocity vector (size [k])
 * prhs[11] = [theta_v], angle of the MS velocity vector (size[k])
 * prhs[12] = [h_in], input matrices (size [u][s][n][tn][k])
 * prhs[13] = [input_phase_los], (size [k])
 * prhs[14] = [k_FACTOR], (size [k])
 * prhs[15] = u, number of MS antennas
 * prhs[16] = s, number of BS antennas
 * prhs[17] = n, number of multipaths
 * prhs[18] = k, number of links
 * prhs[19] = tn, number of time samples
 * Output arguments are 
 * plhs[0] = [h], output array of the coefficients with LOS terms (size [u][s][n][tn][k])
 * plhs[1] = [output_phase_los] pointer to output phases (size [k])

 *
 */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

   /*pointers for input parameters*/

   /* variables for general mode coefficients */
   double *re_G_BS, *im_G_BS, *re_G_MS, *im_G_MS, *phase;

   /* variables for polarized mode coefficients */
   double *re_X_BS_v, *im_X_BS_v, *re_X_BS_h, *im_X_BS_h, *re_X_MS_v, *im_X_MS_v, *re_X_MS_h, *im_X_MS_h,
      *phase_v_v, *phase_v_h, *phase_h_v, *phase_h_h, *r_n1, *r_n2;

   /* variables for los mode */
   double *k_FACTOR, *theta_BS, *theta_MS, *phase_los;

   /* variables for more or less common parameters */
   double k_CONST;
   double *d_s, *d_u, *aod, *aoa, *sq_Pn, *v, *theta_v;
   double *ts;
   double *lm, *ln;  /* subpath ordering */
   int u, s, n, l, m, tn, k; /* parameters to be extracted from inputs' dimensions */
   int h_dims[5]; /* table for output array dimensions */
   int out_Phase_dims[3];
   long int look_up_points;
   int error, mode, GainsAreScalar;

   /*pointers to output*/
   double *re_h, *im_h, *output_SubPathPhase, *output_los_phase;

   /* check which mode is used */
   mode = (int)mxGetScalar(prhs[0]);


   /*Check if general coefficients are wanted */
   if (mode == GENERAL) {

      if(nrhs != 24)
         mexErrMsgTxt("scm_mex_core error: exactly 24 input arguments required for GENERAL option");

      /* ripping values from input */
      re_G_BS = (double*)mxGetPr(prhs[1]);
      im_G_BS = (double*)mxGetPi(prhs[1]);
      re_G_MS = (double*)mxGetPr(prhs[2]);
      im_G_MS = (double*)mxGetPi(prhs[2]);	
      aod = (double*)mxGetPr(prhs[3]);
      aoa = (double*)mxGetPr(prhs[4]);
      d_s = (double*)mxGetPr(prhs[5]);
      d_u = (double*)mxGetPr(prhs[6]);
      phase = (double*)mxGetPr(prhs[7]);
      ts = (double*)mxGetPr(prhs[8]);
      k_CONST = (double)mxGetScalar(prhs[9]);
      v = (double*)mxGetPr(prhs[10]);
      theta_v = (double*)mxGetPr(prhs[11]);
      sq_Pn = (double*)mxGetPr(prhs[12]);
      look_up_points = (long int)mxGetScalar(prhs[13]);
      u = (int)mxGetScalar(prhs[14]);
      s = (int)mxGetScalar(prhs[15]);
      n = (int)mxGetScalar(prhs[16]);
      l = (int)mxGetScalar(prhs[17]);
      m = (int)mxGetScalar(prhs[18]);
      k = (int)mxGetScalar(prhs[19]);
      tn = (int)mxGetScalar(prhs[20]);
      GainsAreScalar = (int)mxGetScalar(prhs[21]);
      lm = (double*)mxGetPr(prhs[22]);
      ln = (double*)mxGetPr(prhs[23]);

      /* now we know the output dimensions */
      h_dims[0] = u;
      h_dims[1] = s;
      h_dims[2] = n*l;
      h_dims[3] = tn;
      h_dims[4] = k;
      out_Phase_dims[0] = k;
      out_Phase_dims[1] = n;
      out_Phase_dims[2] = m;

      /* setting up output */
      plhs[0] = mxCreateNumericArray(5, h_dims, mxDOUBLE_CLASS, mxCOMPLEX);
      plhs[1] = mxCreateNumericArray(3, out_Phase_dims, mxDOUBLE_CLASS, mxREAL);
      re_h = (double*)mxGetPr(plhs[0]); /* Create a C pointer to real values of the output array */
      im_h = (double*)mxGetPi(plhs[0]); /* Create a C pointer to imaginary values of output */
      output_SubPathPhase = (double*)mxGetPr(plhs[1]);

      /* calculating the values */
      error = scm_sum(look_up_points, u, s, n, l, m, k, re_G_BS, im_G_BS, re_G_MS, im_G_MS,
         k_CONST, d_s, d_u, aod, aoa, phase, v, theta_v, ts, tn, sq_Pn, ln, lm, GainsAreScalar,
         re_h, im_h, output_SubPathPhase);

   } /* end general coefficients mode */

   /* check for polarized mode */
   else if (mode == POLARIZED)	{

      if(nrhs != 31)
         mexErrMsgTxt("scmmex error: exactly 31 input arguments required");

      /* ripping values from input */
      re_X_BS_v = (double*)mxGetPr(prhs[1]);
      im_X_BS_v = (double*)mxGetPi(prhs[1]);
      re_X_BS_h = (double*)mxGetPr(prhs[2]);
      im_X_BS_h = (double*)mxGetPi(prhs[2]);	
      re_X_MS_v = (double*)mxGetPr(prhs[3]);
      im_X_MS_v = (double*)mxGetPi(prhs[3]);
      re_X_MS_h = (double*)mxGetPr(prhs[4]);
      im_X_MS_h = (double*)mxGetPi(prhs[4]);
      aod = (double*)mxGetPr(prhs[5]);
      aoa = (double*)mxGetPr(prhs[6]);
      d_s = (double*)mxGetPr(prhs[7]);
      d_u = (double*)mxGetPr(prhs[8]);
      phase_v_v = (double*)mxGetPr(prhs[9]);
      phase_v_h = (double*)mxGetPr(prhs[10]);
      phase_h_v = (double*)mxGetPr(prhs[11]);
      phase_h_h = (double*)mxGetPr(prhs[12]);
      r_n1 = (double*)mxGetPr(prhs[13]);
      r_n2 = (double*)mxGetPr(prhs[14]);
      ts = (double*)mxGetPr(prhs[15]);
      k_CONST = (double)mxGetScalar(prhs[16]);
      v = (double*)mxGetPr(prhs[17]);
      theta_v = (double*)mxGetPr(prhs[18]);
      sq_Pn = (double*)mxGetPr(prhs[19]);
      look_up_points = (long int)mxGetScalar(prhs[20]);
      u = (int)mxGetScalar(prhs[21]);
      s = (int)mxGetScalar(prhs[22]);
      n = (int)mxGetScalar(prhs[23]);
      l = (int)mxGetScalar(prhs[24]);
      m = (int)mxGetScalar(prhs[25]);
      k = (int)mxGetScalar(prhs[26]);
      tn = (int)mxGetScalar(prhs[27]);
      GainsAreScalar = (int)mxGetScalar(prhs[28]);
      lm = (double*)mxGetPr(prhs[29]);
      ln = (double*)mxGetPr(prhs[30]);

      /* now we know the output dimensions */
      h_dims[0] = u;
      h_dims[1] = s;
      h_dims[2] = n*l;
      h_dims[3] = tn;
      h_dims[4] = k;
      out_Phase_dims[0] = k;
      out_Phase_dims[1] = n;
      out_Phase_dims[2] = m;

      /* setting up output */
      plhs[0] = mxCreateNumericArray(5, h_dims, mxDOUBLE_CLASS, mxCOMPLEX);
      plhs[1] = mxCreateNumericArray(3, out_Phase_dims, mxDOUBLE_CLASS, mxREAL);
      re_h = (double*)mxGetPr(plhs[0]); /* Create a C pointer to real values of the output array */
      im_h = (double*)mxGetPi(plhs[0]); /* Create a C pointer to imaginary values of output */
      output_SubPathPhase = (double*)mxGetPr(plhs[1]);

      /* calculating the values */
      error = scm_pol_sum(look_up_points, u, s, n, l, m, k, re_X_BS_v, im_X_BS_v, re_X_BS_h, im_X_BS_h, 
         re_X_MS_v, im_X_MS_v, re_X_MS_h, im_X_MS_h, k_CONST, d_s, d_u, aod, aoa,
         phase_v_v, phase_v_h, phase_h_v, phase_h_h,	r_n1, r_n2,
         v, theta_v, ts, tn, sq_Pn, ln, lm, GainsAreScalar, re_h, im_h, output_SubPathPhase);
   } /* end polarized mode */


   /* check los mode */
   else if (mode == LOS) {

      if(nrhs != 20)
         mexErrMsgTxt("scmlosmex error: exactly 20 input arguments required");

      /* ripping values from input */
      re_G_BS = (double*)mxGetPr(prhs[1]);
      im_G_BS = (double*)mxGetPi(prhs[1]);
      re_G_MS = (double*)mxGetPr(prhs[2]);
      im_G_MS = (double*)mxGetPi(prhs[2]);
      theta_BS = (double*)mxGetPr(prhs[3]);
      theta_MS = (double*)mxGetPr(prhs[4]);
      d_s = (double*)mxGetPr(prhs[5]);
      d_u = (double*)mxGetPr(prhs[6]);
      phase_los = (double*)mxGetPr(prhs[7]);
      ts = (double*)mxGetPr(prhs[8]);
      k_CONST = (double)mxGetScalar(prhs[9]);
      v = (double*)mxGetPr(prhs[10]);
      theta_v = (double*)mxGetPr(prhs[11]);
      re_h = (double*)mxGetPr(prhs[12]);
      im_h = (double*)mxGetPi(prhs[12]);
      output_los_phase = (double*)mxGetPr(prhs[13]);
      k_FACTOR = (double*)mxGetPr(prhs[14]);
      u = (int)mxGetScalar(prhs[15]);
      s = (int)mxGetScalar(prhs[16]);
      n = (int)mxGetScalar(prhs[17]);
      k = (int)mxGetScalar(prhs[18]);
      tn = (int)mxGetScalar(prhs[19]);

      /* setting up output */
      plhs[0] = (mxArray*)prhs[12];
      plhs[1] = (mxArray*)(prhs[13]);

      /* calculating the values */
      error = scm_los(u, s, n, k, re_G_BS, im_G_BS,
         re_G_MS, im_G_MS, k_CONST, d_s, d_u, theta_BS, theta_MS,
         phase_los, v, theta_v, ts, tn, k_FACTOR, re_h, im_h, output_los_phase);

      } /* end los mode */

}