#include <stdio.h>
#include <math.h>
#include <float.h>
#include <io.h>
#include <dos.h>
#include <stdlib.h>
#include <process.h>
#include <string.h>

#define slag 0.23               /* lag between the acquisition of the wind and scalar data, in seconds*/
#define samplerate 10.0         /* samples per second */
#define boom 0                  /* boom angle for the sonic anemometer */
#define nch 6                   /* number of channels */

/* ************************************************************

                raw_flux_grass.C 

                The program processes RAW data files and despikes the data. It produces *.flx, *.nrt and *.err files
                that can be processed in the normal manner with programs that merge met, soil and flux data.

                The input raw data are not adjusted for lags between velocity and scalar variables. So we call the subroutine
                lag_adjust.
                
                
                2/15/2001 

                Dennis Baldocchi
                Ecosystem Science Division
                Department of Environmental Science, Policy and Management
                151 Hilgard Hall
                University of California, Berkeley
                Berkeley, CA 94720
                baldocchi@nature.berkeley.edu

                Vaira Grassland


                debug notes

                1/6/2001
                new file names

                12/4/2000

                too much v data was despiked at night when var v was small and
                atmospheric stratification was stable. Need to up the threshold

                Also debug running variances. There were errors in the calcs
                

************************************************************ */


/* declare subroutine names */

        void filename();
        void indata();         /* reads the 2 byte binary integer files from *.raw */
        void lag_adjust();     /* adjust the raw data for the lag in the data */                                                                                                         
        void despike();        /* removes spikes from the data stream */
        void reynolds_ave();   /* computes flux statistics using Reynolds averaging */
        void output();         /* outputs the data */                                                                                                      
        void zero();           /* zero the arrays */                                                                                                         


     /* declare structures */


struct output
{
/* output variables */

        long int endtime, sample;
        int end_day, end_hour, end_hourmin, end_minute;

        float ww,uu,vv,uw,uv,wv;
        float vari[nch],x_bar[nch],wx[nch],ux[nch],vx[nch],wxr[nch];
        float theta,eta,u_bar,w_bar,v_bar;
        float uur,vvr,wwr,uwr,wvr,u_star;
        float skew[nch],kurt[nch];
        float wdir;

        long int spkcnt[nch];
        float spkave[nch];

} out;


struct input
{
/* input variables */

        short int intdata[nch];

        float fltdata[19000][nch], newdata[19000];
} in;
                                                                                                                                    

struct filename
{
/* structure for filenames */
    
        char flxbuff[39], nrtbuff[39], errbuff[39], spkbuff[39];
}name;


/* declare global variables */

        int handle, handle1, day, fcheck;                                                                                                                         
        long int daytime;
        float timeconst;
        char fname[39], tname[8];
       

/*                                                                                                                                  
File pointers                                                                                                                       
*/                                                                                                                                  
                                                                                                                                    
        FILE *fptr1, *fptr3, *fptr4, *fptr5, *fptr6, *fptr7, *fptr9;

        
        main ()
{
        int hourmin, hour, minute, oldday;
        int j,i;

      
        fcheck=0;
        oldday=190;

     
/*

open list of raw filenames

This file is created from the most recent raw files. It is in ascii form as:

d:\vaira\raw\g3101100.raw
d:\vaira\raw\g3101130.raw

etc......

file296.raw  *
file310.raw  *
file324.raw  *
file338.raw  *
file361.raw  *

*/

        fptr1=fopen("c:\\vaira_2000\\raw\\file361.raw","rt");

/*
    INPUT DATA FROM DISK, Read consequtive files to construct hour average
*/

        handle1=fileno(fptr1);         
                                    
/* while file.raw is open */

        while(! feof(fptr1))
{

/* read RAW file */

        fscanf(fptr1,"%s",fname);

/* strip the string to find the daytime */

        j=0;
        
        for(i=19;i<=25; i++)
        {
        tname[j]=fname[i];
        j++;
        }

        daytime=atoi(tname);

        day=daytime/10000;  
        hour=daytime-(daytime/10000)*10000;  
        hour=(hour/100);  

        printf("%2s\n",tname);
        daytime=atoi(tname);
        hourmin=daytime-(daytime/10000)*10000;  
        hour=(hourmin/100);
        minute=hourmin-hour*100;

/* compute end time */

        out.end_day=day;
        out.end_hour=hour;

        if (minute < 30)
        out.end_minute=30;
        else
        out.end_minute=0;

        if(out.end_minute==0)
        out.end_hour+=1;

        if(out.end_hour >=24)
        {out.end_day += 1;
        out.end_hour=0;
        out.end_minute=0;
        }

        out.end_hourmin=out.end_hour*100+out.end_minute;
        
        out.endtime=out.end_day*10000+out.end_hourmin;

/* create new file if hourmin is 0030 */

        if(out.end_hour==0 && out.end_minute == 30) fcheck = 0;  /* new day and new file */

/* create new file if newday but some other hour, as if there was an interruption between days */

        if (out.end_day > oldday && out.end_hourmin > 30)
        fcheck = 0; 
        
        oldday=out.end_day;
        
      

        filename();
   
/*
    read data from raw files

    Remember element 0 is channel 1, etc
*/

        indata();

/* adjust the time series for the lag between scalars and velocity measurements */

        lag_adjust();

/* examine data files and remove spikes */

        timeconst=600.;
        
        despike();

/* compute turbulent statistics and flux covariances with Reynolds averaging */
     
        reynolds_ave();

/* write to output files */

         output();
       
         zero();
       
/* end of file.raw loop */
}

           
 /* END OF MAIN */
 
}

void filename()
{

        int i;
        
        char fileflx[35], filenrt[35], fileerr[35], filespk[35];
        char daybuff[4];

		itoa(day,daybuff,10);
	
	 		
		
		/* define flux file name as vryr_day.flx, x is system or expt number  */
 
     	
		strcpy(fileflx,"c:\\vaira_2000\\flux\\vr00_");
        strcat(fileflx,daybuff);
        memmove(name.flxbuff,fileflx,29);
        strcat(name.flxbuff, ".flx_dspk");
        printf("%s %s \n",fileflx, name.flxbuff); 


 /* define flux file name as vryr_day.nrt  */
 
       strcpy(filenrt,"c:\\vaira_2000\\flux\\vr00_");
       strcat(filenrt,daybuff);
       memmove(name.nrtbuff,filenrt,29);
       strcat(name.nrtbuff,".nrt_dspk");
       printf("%s  %s\n",filenrt,name.nrtbuff);

 /* define flux file name as  vryr_day.err  */
 
        strcpy(fileerr,"c:\\vaira_2000\\flux\\vr00_");
        strcat(fileerr,daybuff);
        memmove(name.errbuff,fileerr,29);
        strcat(name.errbuff,".err_dspk");
        printf("%s  %s\n",fileerr,name.errbuff);

/* define file name for quantifying spikes as vryr_day.spk */

        strcpy(filespk,"c:\\vaira_2000\\flux\\vr00_");
        strcat(filespk,daybuff);
        memmove(name.spkbuff,filespk,29);
        strcat(name.spkbuff,".spk");
        printf("%s  %s\n",filespk,name.spkbuff);
       



/* output data */


        i = access(name.flxbuff,0);    /* check if file exists         */
        if(i < 0) fcheck = 0;   /* set set flag for newfile     */
        
        if(fcheck==0)
        {
        if((fptr4=fopen(name.flxbuff,"w"))==NULL)
        printf("can't open flxbuff\n");              

        if((fptr5=fopen(name.nrtbuff,"w"))==NULL)
        printf("can't open nrtbuff\n");
        
        if((fptr6=fopen(name.errbuff,"w"))==NULL)
        printf("can't open errbuff\n"); 

        if((fptr7=fopen(name.spkbuff,"w"))==NULL)
        printf("can't open spkbuff\n");
        
        fclose(fptr4);
        fclose(fptr5);
        fclose(fptr6);
        fclose(fptr7);

        fcheck=1;
        }

return;
}

void indata()
{
        int j;
        long int freebyte, ii;

        ii=0;

        if((fptr9=fopen(fname,"rb")) != NULL)
    
        handle=fileno(fptr9);

        freebyte=filelength(handle);


        while (! feof(fptr9)) 
        {
        fread(in.intdata,2,nch,fptr9);

/*
     convert wind velocity from cm/s to m/s 
*/
        for(j=0; j< 3;j++)
        in.fltdata[ii][j]=(float)in.intdata[j]/100.;

/*
      convert integer data from millivolts to floating point volts 
*/
       in.fltdata[ii][3]=(float)in.intdata[3]/100.;
       in.fltdata[ii][4]=(float)in.intdata[4]/1000.;
       in.fltdata[ii][5]=(float)in.intdata[5]/1000.;
        
/*
    w,u and v are channels 0, 1 and 2, respectively. 
    T, Q, CO2 are channels 3, 4, 5
*/
       
        ii++;
 }

        fclose(fptr9);

        out.sample=ii;
      
return;
/* end of indata */
}



void lag_adjust()
{

/* adjust the raw data set for the lag between the sonic anemometer and LICOR LI-7500 data, which is 0.23 seconds */

/* it may be suitable to re-write the code and use pointers to shift the arrays */
    
        int i,j,lag;

        lag=slag*samplerate;  

        for(j=4;j< nch;j++)
        {
                for(i=1;i<=out.sample;i++)
                {
                if(i+lag <out.sample)
                in.newdata[i]=in.fltdata[i+lag][j];
                }

                for(i=1;i<=out.sample;i++)
                in.fltdata[i][j]=in.newdata[i];
        }
            

         
        
return;
}


void despike()
{
        
	
	// This program despikes the raw turbulence data. We had spikes due to lost bits
	// when the operating system performs some background operations and due to the 
	// new LI 7500


		int j, i, iflag[nch];

        int cnt[nch];  /* counts for the running mean variance */
                
        float alpha, onealpha, alpha_fill;
        float runmean[nch], runvar[nch], runstd[nch], runfill[nch], sumvar[nch];
        float delx, fl_delx;

        // compute values for digital filter coefficients. Using a long, 600s, running mean to detect spikes 
        
        alpha=exp(-1./(samplerate*timeconst));
        onealpha=1-alpha;

        // use a short, 10 s, running mean for filling the data; this gives the replacement a more contemporary value 
        
        alpha_fill=exp(-1./(samplerate*10.));
        
        
		// initialize means and counters

		for(j=0; j<nch; j++)
        {
        runmean[j]=0;
        runvar[j]=0;
        sumvar[j]=0;
        runstd[j]=0;
        cnt[j]=100;
        runfill[j]=0;
        out.spkcnt[j]=0;
        out.spkave[j]=0;
        }
        
        // initialize running mean with first 100 values 

        for (j=0; j<nch; j++)
        {
        for (i=1;i<=100;i++)
        runmean[j] += in.fltdata[i][j];
        }       

        for (j=0; j<nch; j++)
        {
         runmean[j] /= 100.;
         runfill[j]=runmean[j];
        }

  // initialize running variance and standard deviation with first 100 values 

        for (j=0; j<nch; j++)
        {
                for (i=1;i<=100;i++)
                {
                delx=in.fltdata[i][j]-runmean[j];
                sumvar[j] += delx*delx;
                }
        }     

        for (j=0; j<nch; j++)
        {
        runvar[j] =sumvar[j]/100.;
        runstd[j] = pow(runvar[j], 0.5);
        }

/*

look for spikes and replace with running mean define spikes as turbulence values exceeding 4 std dev and
scalar values equal or near zero, as when bits are dropped (this criterion does not work for w)

*/   

   for (j=0; j<nch; j++)
   {       
        for(i=1;i<out.sample;i++)
        {

        iflag[j]=0;                
        cnt[j] +=1;
/*
   if sonic values peg above 9.8 m s-1 or a threshold std dev replace them. We have to be careful
   not to throw out small values when w or v are near zero  Hence the threshold of 2 m/s is set.
   Otherwise I was throwing out several thousand v values that were not spiking

   I am also counting spikes and quantifying the mean value of the spikes.
*/

                if(j<3)
                {
                        if(in.fltdata[i][j] > 9.8)
                        {
                        out.spkcnt[j]++;
                        out.spkave[j] +=in.fltdata[i][j];
                        in.fltdata[i][j]=runfill[j];
                        iflag[j]=1;
                        }

				        if (iflag[j] != 1)
						{
						delx=in.fltdata[i][j]-runmean[j];
						fl_delx=fabs(delx);

						// is the fluctuation > six times the std dev?

						if (fl_delx/runstd[j] > 6)
						{out.spkcnt[j]++;
                        out.spkave[j] +=in.fltdata[i][j];
                        in.fltdata[i][j]=runfill[j];
                        iflag[j]=1;}
						}
				   
                } /* end of loop j = 0 to 2 */      

/*
        looking for sonic temperature spikes, with wet sensors getting bad data
        or when u,v,w are bad too. Originally I was not rejecting enough T data and
        computed a lot of bad H's. Changed threshold from 6 to 4.
        
*/
                if (j==3)
                {
                        if(fabs(in.fltdata[i][j]) > fabs(runmean[j]+4.*runstd[j]))
                        {
                        out.spkcnt[j]++;
                        out.spkave[j] +=in.fltdata[i][j];
                        in.fltdata[i][j]=runfill[j];
                        iflag[j]=1;
                        }
                }
                
/* check if q and c are dropping bits (as the LI7500 was observed to have been doing) or spiking out of the 5 volt range */

                if(j>3)
                {
                        if(in.fltdata[i][j] < 0.1 || in.fltdata[i][j] > 4.98)
                        {
                        out.spkcnt[j]++;
                        out.spkave[j] +=in.fltdata[i][j];
                        in.fltdata[i][j]=runfill[j];
                        iflag[j]=1;
                        }

                        if(fabs(in.fltdata[i][j]) > fabs(runmean[j]+6.*runstd[j]))
                        {
                        out.spkcnt[j]++;
                        out.spkave[j] +=in.fltdata[i][j];
                        in.fltdata[i][j]=runfill[j];
                        iflag[j]=1;
                        }
                }

/* update running means, variances and standard deviations for data without spike */
        
   
        if (iflag[j]==1)
        cnt[j]--;
        
        if (iflag[j]==0)
                {
        runmean[j] = runmean[j]*alpha + in.fltdata[i][j]*onealpha;
        delx=in.fltdata[i][j] - runmean[j];

        sumvar[j] += delx*delx;

        runvar[j] = sumvar[j]/cnt[j];
        
        runstd[j]=pow(runvar[j],0.5);
        runfill[j]=runfill[j]*alpha_fill+in.fltdata[i][j]*(1.-alpha_fill);
                }
     
        } /* next i */

} /* next j */


/* evaluate mean values of the spikes */

        for (j=0; j<nch; j++)
        {
        out.spkave[j] /=out.spkcnt[j];
        if(out.spkcnt[j]==0)
        out.spkave[j]=runmean[j];
        }
        
return;
/* end of despike */
}


void reynolds_ave()
{

/*
        this subroutine computes turbulence statistics, variances, covariances, etc with 2 dimensional coordinate
        rotation. Standard Reynolds averaging is used rather than digital filters since we are computing stats from
        scratch.
*/

        float  pi,ce2,se2,ct2,st2,ws,tw,ce,se,ct,st;
        float xp[nch],xn;
        float xp2, xp3, xp4;
        int j,i;

        pi=3.14159;

/* initialize variables */
        
        for (j=0;j<nch;j++)
        {
        out.wx[j]=0;
        out.ux[j]=0;
        out.vx[j]=0;
        out.wxr[j]=0;
        out.x_bar[j]=0;
        out.vari[j]=0;
        out.skew[j]=0;
        out.kurt[j]=0;
        }

        out.ww=0;
        out.uu=0;
        out.vv=0;
        out.uw=0;
        out.uv=0;
        out.wv=0;

/* compute despiked means       */

        for (i=0; i < out.sample; i++)
        {
                for (j = 0; j < nch; j++)
                out.x_bar[j] += in.fltdata[i][j];         
        }


                xn=(float)out.sample;
                
/* compute turbulent statistics, variance, skewness, kurtosis and covariances using arithmetic average */

                for(j=0;j<nch;j++)
                out.x_bar[j] /= xn;
                
                for(i=1;i<out.sample;i++)
                {
                        for(j=0;j<nch;j++)
                        {
                        xp[j]=in.fltdata[i][j]-out.x_bar[j];
                        xp2=xp[j]*xp[j];
                        xp3=xp2*xp[j];
                        xp4=xp3*xp[j];
                        
                        out.vari[j] += xp2;   /* compute variance   */
                        out.skew[j] +=xp3;    /* skewness           */
                        out.kurt[j] += xp4;   /* kurtosis           */


                        out.wx[j] += (in.fltdata[i][0]-out.x_bar[0])*xp[j];  /* wx covariance */
                        out.ux[j] += (in.fltdata[i][1]-out.x_bar[1])*xp[j];  /* ux covariance */
                        out.vx[j] += (in.fltdata[i][2]-out.x_bar[2])*xp[j];  /* vx covariance */
                        }
                

                out.ww += xp[0]*xp[0];        /* velocity variances and covariances */
                out.uu += xp[1]*xp[1];
                out.vv += xp[2]*xp[2];
                out.uw += xp[0]*xp[1];
                out.uv += xp[1]*xp[2];
                out.wv += xp[0]*xp[2];
                }

 /*   compute means of statistics */
                
                for (j=0;j<nch;j++)
                {
                out.vari[j] /= xn;
                out.skew[j] /= (xn*out.vari[j]*sqrt(out.vari[j]));
                out.kurt[j] /= (xn*out.vari[j]*out.vari[j]);    
                out.wx[j] /= xn;
                out.ux[j] /= xn;
                out.vx[j] /= xn;
                }


                out.ww /= xn;
                out.uu /= xn;
                out.vv /= xn;
                out.uw /= xn;
                out.uv /= xn;
                out.wv /= xn;

/*  compute statistics and angles for 2 dimensional coordinate rotation */

                ws = sqrt(out.x_bar[1]*out.x_bar[1] + out.x_bar[2]*out.x_bar[2]);
                tw = sqrt(out.x_bar[0]*out.x_bar[0] + ws*ws);
                ce = out.x_bar[1]/ws;               /* cos eta      */
                se = out.x_bar[2]/ws;               /* sin eta      */
                ct = ws/tw;                     /* cos theta    */
                st = out.x_bar[0]/tw;               /* sin theta    */
                out.theta = 180.*acos(ct)/pi;                /* vertical       */
                out.eta = 180.*acos(ce)/pi;                  /* horizontal     */
                ce2 = ce*ce;
                se2 = se*se;
                ct2 = ct*ct;
                st2 = st*st;

                out.wdir =boom - (atan2(out.x_bar[2],out.x_bar[1]))*180.0/pi;
                if(out.wdir > 360.) out.wdir -= 360.0;
                if(out.wdir < 0.0) out.wdir += 360.0;

                out.u_bar = out.x_bar[1]*ct*ce + out.x_bar[2]*ct*se + out.x_bar[0]*st;
                out.v_bar = out.x_bar[2]*ce - out.x_bar[1]*se;
                out.w_bar = out.x_bar[0]*ct - out.x_bar[1]*st*ce - out.x_bar[2]*st*se;

                out.uur = out.uu*ct2*ce2 + 2*out.uv*ct2*ce*se + 2*out.uw*ct*ce*st +
                      out.vv*ct2*se2 + 2*out.wv*ct*se*st + out.ww*st2;

                out.vvr = out.vv*ce2 + out.uu*se2 - 2*out.uv*ce*se;

                out.wwr = out.ww*ct2 + out.uu*st2*ce2 + out.vv*st2*se2 -
                      2*out.uw*ct*st*ce - 2*out.wv*ct*st*se + 2*out.uv*st2*ce*se;

                out.uwr = out.uw*ce*(ct2-st2) + out.wv*se*(ct2-st2) + out.ww*st*ct -
                      out.uu*st*ct*ce2 - out.vv*st*ct*se2 - 2*out.uv*st*ct*se*ce;
                out.wvr = 0;

                if (out.uwr < 0)
                   out.u_star = sqrt(-out.uwr);
                else
                   out.u_star = 9999;

                for (i = 3; i < nch; i++)
                out.wxr[i] = out.wx[i]*ct - out.ux[i]*st*ce - out.vx[i]*st*se;
return;
/* end of Reynolds_ave */
}

                                                                                                                                                      

void output()
{
        int i;
                        
        /*
        write 'flx' file
        */
                fptr4 = fopen(name.flxbuff,"a");
                fprintf(fptr4,"%3d %02d%02d %4.0f ", out.end_day,out.end_hour,
                                              out.end_minute,out.wdir);
                fprintf(fptr4,"%4.2f %4.3e ",out.u_bar,out.uwr);
                for(i = 3; i < nch; i++)
                   fprintf(fptr4,"%5.3f ", out.x_bar[i]);
                fprintf(fptr4,"%4.3e %4.3e %4.3e",out.uur,out.vvr,out.wwr);
                for(i = 3; i < nch; i++)
                   fprintf(fptr4," %4.3e",out.vari[i]);
                for (i = 3; i < nch; i++)
                    fprintf(fptr4," %4.3e",out.wxr[i]);
                fprintf(fptr4," %u\n",out.sample);
                fclose(fptr4);
        /*
        write unrotated file
        */
                fptr5 = fopen(name.nrtbuff,"a");
                fprintf(fptr5,"%3d %02d%02d %4.0f ",out.end_day,out.end_hour,
                                              out.end_minute,out.wdir);
                fprintf(fptr5,"%5.3f %5.3f %5.3f ",out.x_bar[0],
                out.x_bar[1],out.x_bar[2]);
                fprintf(fptr5," %4.3e %4.3e %4.3e %4.3e %4.3e %4.3e ",
                  out.ww,out.uu,out.vv,out.uw,out.uv,out.wv);
                for (i = 3; i < nch; i++)
                    fprintf(fptr5," %4.3e %4.3e %4.3e %4.3e ",out.wx[i],out.ux[i],
                    out.vx[i],out.vari[i]);
                fprintf(fptr5," %d\n",out.sample);
                fclose(fptr5);
        /*
        write 'err' file
        */
                fptr6 = fopen(name.errbuff,"a");
                fprintf(fptr6,"%3d %02d%02d %4.0f ",out.end_day,out.end_hour,
                                              out.end_minute,out.wdir);
                for (i = 0; i < nch; i++)
                    fprintf(fptr6," %4.3e ",out.skew[i]);
                for (i = 0; i < nch; i++)
                    fprintf(fptr6," %4.3e ",out.kurt[i]);
                fprintf(fptr6," %4.0f %4.0f ", out.theta,out.eta);
                fprintf(fptr6," %d\n",out.sample);
                fclose(fptr6);
               
      /*
      write 'spk file
      */

                fptr7=fopen(name.spkbuff,"a");
                fprintf(fptr6,"%3d %02d%02d ",out.end_day,out.end_hour,
                                              out.end_minute);
                for (i = 0; i < nch; i++)
                    fprintf(fptr7," %d ",out.spkcnt[i]);
                for (i = 0; i < nch; i++)
                    fprintf(fptr7," %f ",out.spkave[i]);

                fprintf(fptr7," \n");
                fclose(fptr7);
return;
/* end of output */
}


void zero()
{
 
      
return;
}