#include <stdio.h>
#include <math.h>
#include <string.h>
#include <complex.h>
using namespace std;

const int Nf = 128*2*2;   //Number of points in the frequency spectrum
const int Mf = Nf;
const int Nt = 1024;
const int NN = 10*Nt;
const int N1 = 128;
const int Nx = 2*N1;
const int maxSize = 15;
const int nColsOut = 15;
const double pi = 3.14159265358979323846;

struct constants
{
   int N;
   int N1;
   int Nm;
   int Nt;
   double oml;
   double omh;
   double dpract;
   double beta1;
   double Hs;
   double Te;
   double ti;
   double tf1;
   double tcut;
   double tend;
};

struct output
{
   double    x[Nx];
   double   he[Nx];
   double   ho[Nx];
   double  etap[NN];
   double etaRc[NN];
   double  fclr[NN];
   double  fcla[NN];
   double  pwra[NN];
   double  pwrc[NN];
   double velnc[NN];
   double    tp[NN];
};



int crest_wave_profile_mbari_multi1_pared( int rowsAlpha, double alpha[], int rows, int cols,
                                           int processedScan[maxSize][maxSize] ); 

int prop_impresp_eta_30s_pred_act_fixedtime_mbari_relt_a_wdispV2( struct output *, double alpha[],
                                                                  struct constants *);
#ifdef __cplusplus
extern "C" {
#endif
double      Fresnel_Sine_Integral( double x );
double      Fresnel_Cosine_Integral( double x );
#ifdef __cplusplus
}
#endif

//float ran0(long *idum)


int main()
{
//
   int error = 0;
   double alpha[Mf];
   
//   
// Read in input array:
   const char *CrestsCsv = "CrestDist222253.csv";
   const char *ImpulseRespCsv = "impresps_mbari.csv";
   const char *AlphaCsvOut = "Alphas.csv";
   const char *BuoyCsvOut = "BuoyOut.csv";
   const char *SeaSurfaceCsvOut = "SeaSurfaceOut.csv";
   FILE *FidCrests   = fopen(CrestsCsv, "r");
   if( NULL == FidCrests )
   {
      printf("predict: Error %d opening file %s.\n", errno, CrestsCsv);
      error= -1;
   }

   FILE *FidImpulse  = fopen(ImpulseRespCsv, "r");
   if( NULL == FidImpulse )
   {
      printf("predict: Error %d opening file %s.\n", errno, ImpulseRespCsv);
      error= -1;
   }
   FILE *FidAlphaOut = fopen(AlphaCsvOut,"w");
   if( NULL == FidAlphaOut )
   {
      printf("predict: Error %d opening file %s.\n", errno, AlphaCsvOut);
      error= -1;
   }

   FILE *FidBuoyOut = fopen(BuoyCsvOut,"w");
   if( NULL == FidBuoyOut )
   {
      printf("predict: Error %d opening file %s.\n", errno, BuoyCsvOut);
      error= -1;
   }

   FILE *FidSeaSurfaceOut = fopen(SeaSurfaceCsvOut,"w");
   if( NULL == FidSeaSurfaceOut )
   {
      printf("predict: Error %d opening file %s.\n", errno, SeaSurfaceCsvOut);
      error= -1;
   }

   if( error ) return error;
   
   int buflen = 1024;
   char buf[buflen];
   char *bufp = buf;
   int Mc=0, Mt=0;
   int num;
   int processedScan[maxSize][maxSize];
   int sret;
   char *token;
   
   while( fgets( buf,  buflen, FidCrests) != NULL && Mc < maxSize )
   {
      Mt = 0;
      bufp = buf;
      printf("%s\n", buf);
      //while( (sret = sscanf(buf, "%d,", &processedScan[Mc][Mt] )) && Mt < maxSize )
      while( (token = strtok( bufp, ",")) != 0 )
      {
         processedScan[Mc][Mt] = atoi(token);
         printf("%d = ps[%d][%d], returned %d\n",processedScan[Mc][Mt], Mc, Mt, sret );
         Mt++;
         bufp = NULL;
      }
      Mc++;
   }
   printf(" Mc = %d, Mt = %d\n", Mc,Mt);

   if( Mc > maxSize )
   {
      printf("Mc = %d exceeds maximum array limit.\n",Mc);
      return -1;
   }
   if( Mt > maxSize )
   {
      printf("Mc = %d exceeds maximum array limit.\n",Mt);
      return -1;
   }

   error = crest_wave_profile_mbari_multi1_pared( Mf, alpha, Mc, Mt, processedScan );
   
   for( int i = 0; i < Mf; i++ )
   {
      fprintf(FidAlphaOut, "%g\n", alpha[i]);
   }

//
// Setup constants for the function below:

   constants params, *pptr;
   pptr = &params;

   params.N1     = N1;
   params.Nm     = 8;
   params.Nt     = Nt;
   params.oml    = 0.18;
   params.omh    = 1.57;
   params.dpract = 475.0;
   params.beta1  = .74;
   params.Hs     = 1.83;
   params.Te     = 10.0;
   params.ti     = 0.0;
   params.tf1    = 5859.7376;
   params.N      = 4096;
   params.tcut   = 500.0;
   params.tend   = 2400.0;

   struct output out;


   //for( int i=0; i<Mf; i++ ) printf("Alpha[%d] = %g\n",i, alpha[i]);


   error = prop_impresp_eta_30s_pred_act_fixedtime_mbari_relt_a_wdispV2( &out, alpha, pptr );
   
//
// Save it to a csv file:

   int i;
   for(i=0; i<Nx; i++)
   {
      fprintf(FidBuoyOut,"%g, %g, %g,\n", out.x[i], out.he[i], out.ho[i]);
   }
   for(i=0; i<NN; i++)
   {
      fprintf(FidSeaSurfaceOut,"%g, %g, %g, %g, %g, %g, %g, %g, \n", out.etap[i], out.etaRc[i],
              out.fclr[i], out.fcla[i], out.pwra[i], out.pwrc[i], out.velnc[i], out.tp[i]);
   }

   
   fclose(FidCrests);
   fclose(FidImpulse);
   fclose(FidAlphaOut);
   fclose(FidBuoyOut);
   fclose(FidSeaSurfaceOut);

   
   return error;
}

int prop_impresp_eta_30s_pred_act_fixedtime_mbari_relt_a_wdispV2( struct output *out, double alpha[],
                                                                  struct constants *pp)
{
//function [etap, velnc, tp, tpred, x] = prop_impresp_eta_30s_pred_act_fixedtime_mbari_relt_a_wdispV2( alpha, params )
// This is the correct version. Both phases and magnitudes agree, as shown
// by predicted and computed wave records,and significant wave height
// comparisons. May 29, 2016.
// Impulse repsponse function for linear wave propagation, up-wave-down-wave
// causalization; fixed-time, snapshot measurement
//
//warning('off','all')
  double g = 9.806; 
  double rho = 1.025e3;
//
  int Nf = pp->N1/2;
//int Nx = pp->N1*2;
  double omeg[Nf], k[Nf], vg[Nf];
  double x[Nx], xx[2*Nx], hlx[2*Nx];
  int i;
  int error = 0;



  for( i=0; i<Nf; i++ )
  {
     omeg[i]=0.;
     k[i]=0.;
     vg[i]=0.;
  }

  for( i=0; i<Nx; i++ )
  {
     x[i]=0.;
  }

  for( i=0; i<2*Nx; i++ )
  {
     xx[i]=0.;
     hlx[i]=0.;
  }

//
  double delom = (pp->omh - pp->oml)/Nf;


//
// Deep water approximation
//
  for(i=0; i<Nf; i++)
  {
    omeg[i] = pp->oml + i*delom;
    k[i]    = omeg[i]*omeg[i]/g;
    vg[i]   = 0.5*omeg[i]/k[i];
  }
//
  double vgmin, vgmax;
  vgmin = vg[0];
  vgmax = vg[0];
  for(i=0; i<Nf; i++ )
  {
     if(vg[i] < vgmin) vgmin = vg[i];
     if(vg[i] > vgmax) vgmax = vg[i];
  }
  double df = 0.0;
//
//
//To avoid dealing with small x values in denominators, set xl = 50.0 (can
//be arbitrary)
//

  double tpred = pp->dpract/vgmax;
//
  double p = tpred;
//
  double dR = vgmin*tpred;
  double xl = dR*1;
  double xzero = xl;
  double xR = vgmax*tpred + xl;
  double xh = xR - dR;
  double ell = xh - xl;
  double xpos = ell+xl;
  double delx = (xpos - xzero)/Nx;
//
  double hl[Nx], he[Nx], ho[Nx];

  for(i=0; i<Nx; i++)
  {
     hl[i] = 0.;
     he[i] = 0.;
     ho[i] = 0.;
  }

  for(i=0; i<Nx; i++)  x[i] = xl + i*delx;

  for(i=0; i<Nx; i++)
  {
     xx[i] = -xl - i*delx;
     xx[i+Nx] = x[i];
  }

  double arg, cofx1, cofx2, cofx, csx;
  
  for(i=0; i<Nx; i++)
  {
     arg = g*p*p/(4*(xR-x[i]));
     cofx1 = sqrt(pi*g/(2*(xR-x[i])));
     cofx2 = p/(xR-x[i]);
     cofx = cofx1*cofx2;
     csx = cos(arg) + sin(arg);
     he[i] = cofx*csx/(2*pi); // correct as of May 28
  }
  
//

  double arg1, arg2, a1, a2, b1, b2, ho1;
  double rootPi2 = sqrt(pi/2.);
  for(i=0; i<Nx; i++)
  {
    arg1 = g*p*p/(4*(xR-x[i]));
    arg2 = sqrt(g/(xR-x[i]))*p/2;
    a1 = cos(arg1);
    b1 = sin(arg1);
//    [a2,b2] = fresnelCS(arg2);
    a2 = Fresnel_Cosine_Integral( rootPi2*arg2 );
    b2 = Fresnel_Sine_Integral( rootPi2*arg2 );
    cofx1 = sqrt(pi*g/(2*(xR-x[i])));
    cofx2 = p/(xR-x[i]);
    cofx = cofx1*cofx2;
    ho1 = cofx*(a1*a2 + b1*b2);
    ho[i] = ho1/pi;
  }

  for(i=0; i<Nx; i++)
  {
     hlx[i] = 0.0;
  }

//
  for(i=0; i<Nx; i++)
  {
     hlx[i+Nx] = he[i] + ho[i];
  }
// %
// use up-wave measurements over a time period to evaluate eta at xB at t;
// for a range of t values
// Generate eta at x = xR
//
// wave amplitude in spectrum
  const int Mf = Nf;

// 2-parameter spectrum based on Newman 1977; Modified Pierson Moskowitz spectrum
  double enom = 2*pi/pp->Te;
  double Ue = pow(  (4/5)*pp->beta1, 0.25  ) * g / enom;
  double multcof = 2*131.5*pow(pp->Hs,2)/pow(pp->Te,4);  
  double multarg = 1054/pow(pp->Te,4);  

  double ommin, ommax;
  ommin = omeg[0];
  ommax = omeg[0];
  for(i=0; i<Nf; i++ )
  {
     if(omeg[i] < ommin) ommin = omeg[i];
     if(omeg[i] > ommax) ommax = omeg[i];
  }

  double spec[Mf], amp[Mf], om[Mf], phase[Mf], etaf[Mf];

  for(i=0; i<Mf; i++)
  {
     spec[i] = 0.;
     amp[i] = 0.;
     om[i] = 0.;
     phase[i] = 0.;
     etaf[i] = 0.;
  }

//
  double dom = omeg[2]-omeg[1];       

// for(i=0; i<Mf; i++)
// {
//    phase[i] = 2.*pi*ran0(1);
// }
//
  for(i=0; i<Mf; i++)
  {
     phase[i] = alpha[i];     
  }

//
//compi = complex(0,1);            
//
  double scof, st;
  for(i=0; i<Mf; i++)
  {
     om[i] = omeg[i];           
     k[i] = pow(omeg[i],2)/g;        
     scof = multcof/pow(om[i],5);    
     st = exp(-multarg/pow(om[i],4));
     spec[i] = scof*st;
  }
//
// End P-M calculation (?).   
//
// spectral moment m_0
  double mom0 = 0.0;
  double sum = 0.0;
  for(i=0; i<Mf; i++)
  {
    sum = sum + spec[i];
  }
  double amom0 = sum*dom/pi;
//

//needed for Pierson Moskowitz spectral density, amplitude calculation 
  for(i=0; i<Mf; i++)
  {
    amp[i] = sqrt(2*spec[i]*dom); // use if simply summing varous frequency components to form
                                  //irregular signal, not inv FT as here.
//   amp[i] = sqrt(spec[i]/2);    // Inverse Fourier transformation with division by delt accounts for dom here.
  }
//
//needed for Pierson Moskowitz spectral density calculation 

//
  std::complex<double> etaf1[Mf];
  std::complex<double> compi(0., 1.);
  std::complex<double> czero(0.,0.);

  for(i=0; i<Mf; i++)  etaf1[i] = czero;
//
//
  for(i=0; i<Mf; i++)
  {
    om[i] = omeg[i];
    k[i]  = pow(omeg[i],2)/g;
//   amp[i] = sqrt(spec[i]*dom/2); // use if simply summing varous frequency components to form
                                   //irregular signal, not inv FT as here.
//   amp[i] = 0.5*sqrt(spec[i]); % Inverse Fourier transformation with division by delt accounts for dom here.
    etaf1[i] = amp[i]*exp(compi*phase[i]);    
  }

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  double etam[Nx][10*Nt], etap[10*Nt] ;
  std::complex<double> etaps[10*Nt];
  std::complex<double> etamxt;
  double t[10*Nt], tp[10*Nt];
  int j, jf;
  
//needed for wave-elevation map based on phases from crestline locations and along predominant direction 
  for(i=0; i<Nx; i++)
  {
     for(j=0; j<10*Nt; j++)
     {
        //etam[i][j] = czero;
        etam[i][j] = 0.;
     }
  }
  for(i=0; i<10*Nt; i++)
  {
     etap[i]  = 0.;
     etaps[i] = czero;
     t[i]     = 0.;
     tp[i]    = 0.;
  }

  double delt = (pp->tf1  -  pp->ti)/(10*Nt);
//
  double etam0, etam01;
  for(i=0; i<10*Nt; i++)
  {
     t[i] = pp->ti + (i+1)*delt;
     //
     etam0 = 0.0;
     for(j=0; j<Nx; j++)
     {
        etam0 = 0.0;
        for(jf=0; jf<Mf; jf++)
        {
           etamxt = exp(-compi*(k[jf]*x[j]-om[jf]*t[i]));
           etam01 = real(etaf1[jf]*etamxt);
           etam0 = etam0 + etam01;
        } // jf
        etam[j][i] = etam0;
     } //j
  } //i
//

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// wave elevation prediction based on etam and hlx
//
//etai = zeros(Nx,10*Nt,1);
  double etap0, etap0j;
  int jj;
  for(i=0; i<10*Nt; i++)
  {
    tp[i] = t[i] + p;
    etap0 = 0.0;
     for(j=0; j<Nx; j++)
     {
        jj = j + Nx;
        etap0j = hlx[jj]*etam[j][i];
        etap0 = etap0+etap0j;
     } //j
    etap[i] = etap0*delx;
  } // i
//
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
  double etaRm0, etaRm01, etaRc[10*Nt], etaRcs[10*Nt];
  for(i=0; i<10*Nt; i++)
  {
    etaRm0 = 0.0;
    for(jf=0; jf<Mf; jf++)
    {
        etamxt = exp(-compi*(k[jf]*xR - om[jf]*tp[i]));
        etaRm01 = real(etamxt*etaf1[jf]);
        etaRm0 = etaRm0 + etaRm01;
    } //jf
    etaRc[i] = etaRm0;
  } //i
//
//
//  int NN = 10*Nt;
//


  int n;
  for( n=0; n<NN; n++ )
  {
     tp[n] = (double) n*delt + p;
     tp[n] = t[n] + p;
  } //n
//
  double tcf = 30.0;
//

  int Ncf = tcf/delt;
  int N2 = pp->N * 2;
//
//
  double hrct[N2], hactt[N2], hactrt[N2], hldat[N2], hldrt[N2], hdct[N2];
  for( i=0; i<N2; i++ )
  {
     hrct[i]   = 0.;
     hactt[i]  = 0.;
     hactrt[i] = 0.;
     hldat[i]  = 0.;
     hldrt[i]  = 0.;
     hdct[i]   = 0.;
  }
//
  double vro[NN], vract[NN], vractr[NN], vracta[NN],velnc[NN];
  double fcla[NN], fclr[NN], pwrc[NN],pwra[NN], etap1[NN];
     
  for( i=0; i<NN; i++ )
  {
     vro[i] = 0.;
     vract[i] = 0.;
     vractr[i] = 0.;
     vracta[i] = 0.;
     velnc[i] = 0.;
//
     fcla[i] = 0.;
     fclr[i] = 0.;
//
     pwrc[i] = 0.;
     pwra[i] = 0.;

     etap1[i] = 0.;
  }
//
//   
//Read in impulse response array:
  const char *ImpulseRespCsv = "impresps_mbari.csv";
  FILE *FidImpulse  = fopen(ImpulseRespCsv, "r");
  int buflen = 1024;
  char buf[buflen];
  char *bufp = buf;
  char *token;
  int nRows = 0;
  int maxRows = N2;

  double temp1, temp2, temp3, temp4, temp5, temp6;
  int scannum;

  while( fgets( buf,  sizeof(buf), FidImpulse) != NULL && nRows < maxRows)
  {
//   printf("%s\n", buf);
     scannum = sscanf(buf, "%le,%le,%le,%le,%le,%le",
            &temp1,
            &temp2,
            &temp3,
            &temp4,
            &temp5,
            &temp6);
//   printf("sscanf read %d numbers.\n",scannum);
     hrct[nRows] = temp1;
     hactt[nRows] = temp2;
     hactrt[nRows] = temp3;
     hldat[nRows] = temp4;
     hldrt[nRows] = temp5;
     hdct[nRows] = temp6;
//   printf("impresps = [%g %g %g %g %g %g]\n", hrct[nRows], hactt[nRows], hactrt[nRows],
//          hldat[nRows], hldrt[nRows], hdct[nRows]);
     nRows++;
  }
  printf("Read %d rows from %s.\n", nRows, ImpulseRespCsv);


//for( i=0; i<nRows; i++ )
//   printf("impresps = [%g %g %g %g %g %g]\n", hrct[i], hactt[i], hactrt[i], hldat[i], hldrt[i], hdct[i]);


/*
load impresps_mbari.csv;
fname1 = impresps_mbari;
hrct = fname1(:,1);
hactt = fname1(:,2);
hactrt = fname1(:,3);
hldat = fname1(:,4);
hldrt = fname1(:,5);
hdct = fname1(:,6);
*/


  
//
//
  int m;
  double sumvo, sumvrt, sumvat, sumfla, sumflr, sumdo;
  int N = pp->N;
  for( n=0; n<NN; n++ )
  {
    etap0 = 0.0;
    for( j=0; j<Nx; j++ )
    {
       jj = j + Nx;
       etap0j = hlx[jj]*etam[j][n];
       etap0 = etap0+etap0j;
    } //j
    etap1[n] = etap0*delx;
    //
    if( n >= 3*Ncf )
    {
        sumvo = 0.0;
        sumvrt = 0.0;
        sumvat = 0.0;
        sumfla = 0.0;
        sumflr = 0.0;
        sumdo = 0.0;
        //
        for( m=0; m<3*Ncf; m++)
        {
           /*
            sumvo  = sumvo  + hrct[N-Ncf+m]  *etap1[n-m+1];
            sumvrt = sumvrt + hactrt[N-Ncf+m]*etap1[n-m+1];
            sumvat = sumvat + hactt[N-Ncf+m] *etap1[n-m+1];
            sumfla = sumfla - hldat[N-Ncf+m] *etap1[n-m+1];
            sumflr = sumflr + hldrt[N-Ncf+m] *etap1[n-m+1];
            sumdo  = sumdo  + hdct[N-Ncf+m]  *etap1[n-m+1];
           */

            sumvo  = sumvo  + hrct[N-Ncf+m]  *etap1[n-m];
            sumvrt = sumvrt + hactrt[N-Ncf+m]*etap1[n-m];
            sumvat = sumvat + hactt[N-Ncf+m] *etap1[n-m];
            sumfla = sumfla - hldat[N-Ncf+m] *etap1[n-m];
            sumflr = sumflr + hldrt[N-Ncf+m] *etap1[n-m];
            sumdo  = sumdo  + hdct[N-Ncf+m]  *etap1[n-m];


        }
        //
        vro[n] = sumvo*delt;
        vractr[n] = sumvrt*delt;
        vracta[n] = sumvat*delt;
        fcla[n] = sumfla*delt;
        fclr[n] = sumflr*delt;
        velnc[n] = sumdo*delt;
    }
    vract[n] = vractr[n] + vracta[n];
  }//n
//

//
//
  for( n=0; n<NN; n++ )
  {
    pwrc[n] = fclr[n]*vract[n];
    pwra[n] = fcla[n]*vract[n];
  }
//
  double sumeta = 0.0;
  int NN1 = 0;
  for( n=0; n<NN; n++ )
  {
     if( fabs(etap[n]) > 0)
     {
        sumeta = sumeta + pow(etap[n],2);
        NN1 = NN1 + 1;
     }
  }
  double fNN = NN1;
  double Hs1 = 4.0*sqrt(sumeta/fNN);
//
  double sumeta1 = 0.0;
  NN1 = 0;
  for( n=0; n<NN; n++ )
  {
     if( abs(etam[100][n]) > 0 )
     {
        sumeta1 = sumeta1 + pow(etam[100][n],2);
        NN1 = NN1 + 1;
     }
  }
  fNN = NN1;
  double Hs3 = 4.0*sqrt(sumeta1/fNN);
//
  for( n=0; n<NN; n++ )
  {
     etaRc[n] = 0.;
     etaRcs[n] = 0.;
  }
//
  for( i=0; i<NN; i++ )
  {
    etaRm0 = 0.0;
    for( jf=0; jf<Mf; jf++ )
    {
       etamxt = exp(-compi*(k[jf]*xR - om[jf]*tp[i]));
       etaRm01 = real(etamxt*etaf1[jf]);
       etaRm0 = etaRm0 + etaRm01;
    } //jf
    etaRc[i] = etaRm0;
  } //i

//
  sumeta = 0.0;
  NN1 = 0;
  for( n=0; n<NN; n++ )
  {
     if( fabs(etaRc[n]) > 0 )
     {
        sumeta = sumeta + pow(etaRc[n],2);
        NN1 = NN1 + 1;
     }
  }
  fNN = NN1;
  double Hs2 = 4.0*sqrt(sumeta/fNN);
//
//  etaps = etap1;
//  etaRcs = etaRc;
//
//
//

  int Ncut = pp->tcut/delt;
  int Nend = pp->tend/delt;
  double sumpc = 0.0;
  double sumpam = 0.0;
  for( n=0; n<NN; n++ )
  {
     if( (n > Ncut) && (n < Nend) )
     {
        sumpc = sumpc + pwrc[n];
        sumpam = sumpam + pwra[n];
     }
  }
//
  double  Nsum = Nend-Ncut;
  double  pwrcav = sumpc/Nsum;
//
  for(i=0; i<Nx; i++ )
  {
     out->x[i]     = x[i];
     out->he[i]    = he[i];
     out->ho[i]    = ho[i];
  }
  for(i=0; i<NN; i++ )
  {
     out->etap[i]  = etap[i];
     out->etaRc[i] = etaRc[i];
     out->fclr[i]  = fclr[i];
     out->fcla[i]  = fcla[i];
     out->pwra[i]  = pwra[i];
     out->pwrc[i]  = pwrc[i];
     out->velnc[i] = velnc[i];
     out->tp[i]    = tp[i];
  }
return error;
}







int crest_wave_profile_mbari_multi1_pared( int Mf, double alpha[], int Mc, int Mt, int processedScan[maxSize][maxSize] ) 
{
   int error = 0;
   double g = 9.806;
   double rho = 1.025e3;
//
// alpha has dimension Mf.
//   
// Mc is the number of crests considered in each scan.
// Mt is the number of scans.
// processedScan is Mc x Mt
//
   double oml = 0.18;     //Low and high cut-off frequencies, radian/sec. 
   double omh = 1.57;
   double delom;
   //int Mc = 3;         //Number of crests
   //int Mt = 8;         //Number of scans  
   double delt = 1.26;    //time between scans + 1/2 the time needed for 1 scan.
   int mf, mc, mt;

/*

   size_t Mc = mxGetM(prhs[0]);
   size_t Mt = mxGetN(prhs[0]);
   mexPrintf("Mc = %d, Mt = %d\n", Mc, Mt);
   mxDouble *processedScan = mxGetDoubles(prhs[0]);
//
// Create intermediate variables:   
   mxDouble *omeg = mxGetDoubles(mxCreateDoubleMatrix(Mf,1,mxREAL));
   mxDouble *k    = mxGetDoubles(mxCreateDoubleMatrix(Mf,1,mxREAL));
   mxDouble *tm   = mxGetDoubles(mxCreateDoubleMatrix(Mt,1,mxREAL));
   mxDouble *xc   = mxGetDoubles(mxCreateDoubleMatrix(Mc,Mt,mxREAL));
   mxDouble *yc   = mxGetDoubles(mxCreateDoubleMatrix(Mc,Mt,mxREAL));
   mxDouble *Ac   = mxGetDoubles(mxCreateDoubleMatrix(Mf,1,mxREAL));
   mxDouble *Bc   = mxGetDoubles(mxCreateDoubleMatrix(Mf,1,mxREAL));
   mxDouble *C    = mxGetDoubles(mxCreateDoubleMatrix(Mf,1,mxREAL));
//   
// Create output array:
   plhs[0] = mxCreateDoubleMatrix(Nf, 1,  mxREAL);
   mxDouble *alpha = mxGetPr(plhs[0]);
   int i;
*/

   double omeg[Mf];
   double k[Mf];
   double tm[Mt];
   double xc[Mc][Mt];
   double yc[Mc][Mt];
   double Ac[Mf];
   double Bc[Mf];
   double C[Mf];
// double alpha[Nf];


   for(int i=0; i<Nf; i++) omeg[i] = 0.;

   delom = (omh-oml)/( (double) Nf );

   for( mf = 0; mf < Mf; mf++ )
   {
      omeg[mf] = oml + mf*delom;
      printf("omeg[%d] = %.4f\n",mf,omeg[mf]);
   }

   for( mf = 0; mf < Mf; mf++ )
   {
      k[mf] = omeg[mf]*omeg[mf]/g;
      printf("k[%d] = %.4f\n",mf,k[mf]);
   }

   for( mt = 0; mt < Mt; mt++ ) tm[mt] = 0.;

   for( mt = 0; mt < Mt; mt++ )  tm[mt] = mt*delt;

//
// fname1 = processedScan;

//xc = zeros(Mc,Mt,1);
//yc = zeros(Mc,Mt,1);

   int i;

   for( mc = 0; mc < Mc; mc++ )
      for( mt = 0; mt < Mt; mt++ )
         xc[mc][mt] = 0.;

   //for( i = 0; i<Mc*Mt; i++ )  xc[i] = 0.;

   for( mc = 0; mc < Mc; mc++ )
      for( mt = 0; mt < Mt; mt++ )
         yc[mc][mt] = 0.;

   //for( i = 0; i<Mc*Mt; i++ )  yc[i] = 0.;

   
 for( mc = 0; mc < Mc; mc++ )
 {
    for( mt = 0; mt < Mt; mt++ )
    {
       xc[mc][mt] = -processedScan[mc][mt];
       yc[mc][mt] = .5;
    }
 }

   // for( i = 0; i < Mc*Mt; i++ )
   // {
   //    xc[i] = -processedScan[i];
   //    yc[i] = .5;
   // }

   
   for( mf = 0; mf < Mf; mf++ )
   {
      Ac[mf] = 0.;
      Bc[mf] = 0.;      // set to unity
      alpha[mf] = 0.;   // phases
      C[mf] = 0.;       // amplitudes
   }

//

   size_t index;
   
   for( int n = 0; n < Mf; n++ )
   {
      double sumnum = 0.0;
      double sumden = 0.0;
      double sumnt, sumdt, cterm, sterm;
      int m;
      for( m = 0; m<Mc; m++ )
      {
         sumnt = 0.0;
         for( mt = 0; mt<Mt;  mt++ )
         {
            //index = SUB2IND_2D(Mc,Mt,m,mt);
            cterm = cos(k[n]*xc[m][mt] - omeg[n]*tm[mt]);
            //cterm = cos(k[n]*xc[index] - omeg[n]*tm[mt]);
            sumnt = sumnt + cterm;
         } //mt
         sumnum = sumnum + sumnt/Mt;
      } //m
      //
      for( m = 0; m<Mc; m++ )
      {
         sumdt = 0.0;
         for( mt = 0; mt<Mt;  mt++ )
         {
            //index = SUB2IND_2D(Mc,Mt,m,mt);
            sterm = sin(k[n]*xc[m][mt] - omeg[n]*tm[mt]);
            //sterm = sin(k[n]*xc[index] - omeg[n]*tm[mt]);
            sumdt = sumdt + sterm;
         } //mt
         sumden = sumden + sumdt/Mt;
      }
      //
      Ac[n] = sumnum/sumden;
      alpha[n] = atan2(1.0,Ac[n]);

   }

   return 0;
}