#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>


#define STAT(y,sum,sum2)   \
	sum+=y;                \
	sum2+=pow(y,2);        \
    if(fabs(maxabs)<fabs(y)){maxabs=fabs(y);}\

#define MAX 1000000

void printfile(void);
void freq_amp(void);

unsigned long i;
unsigned long limit;


double dt;

double x, x1;
double sr;
double tlimit;

double freq;

long nv=0;

long nzc;

double y;

double sum=0.;
double sum2=0.;

double ave, sd, grms;

double maxabs=0.;

double peak=0.;

double ybefore;

int j;

int icheck, numfiles;

double ff[MAX];
double pp[MAX];


char filename[10][FILENAME_MAX];
FILE *pFile[10];


int main( int argc, char *argv[] )
{

	printf("\n sine_sweep_spectral.cpp   ver 1.0  March 1, 2010 \n");
	printf("\n by Tom Irvine   Email: tomirvine@aol.com  \n");


   numfiles =4;

   printf("\n  argc= %d \n",argc);


   if(argc > 1 )
   {
		strcpy (filename[0], argv[1] );
   }

   if(argc == 1 )
   {
/*
      printf("\n  Command-line mode."
	  " User failed to enter input file and limit.\n");
      exit(1);

*/

		printf("\n The input file must be: time(sec) and amplitude(units) ");
		printf("\n Format is free, but no header lines allowed.\n");
		printf("\n Please enter the input filename. \n");


		scanf("%s",filename[0]);

//	   printf("\n Enter segment duration (sec) \n");

//       scanf ("%lf", &tlimit);

   }

   if(argc == 2 )
   {

//      tlimit = .2;

	   printf("\n Enter segment duration (sec) \n");

       scanf ("%lf", &tlimit);

   }


   if(argc == 3 )
   {
      sscanf (argv[2], "%lf", &tlimit);
   }

/*

	printf("\n\n  argc = %d ", argc );
	printf("\n argv[0] %s ", argv[0] );
	printf("\n argv[1] %s ", argv[1] );

*/

	printf("\n\n");

	pFile[0] = fopen(filename[0], "rb");

	while(pFile[0] == NULL )
	{

		printf("\n Failed to open file: %s \n", filename[0]);
		printf("\n Please enter the input filename: \n");

		scanf("%s",filename[0]);

		pFile[0] = fopen(filename[0], "rb");

	}

   strcpy (filename[1], "s.dat");
   strcpy (filename[2], "ss.dat");
   strcpy (filename[3], "sa.dat");
   strcpy (filename[4], "sp.dat");
   strcpy (filename[5], "time_freq.dat");
   strcpy (filename[6], "time_amp.dat");
   strcpy (filename[7], "freq_amp.dat");
   strcpy (filename[8], "time_freq_amp.dat");

   for(i=1;i<=8;i++)
   {
      pFile[i] = fopen(filename[i], "w");
   }


   for(j=0; j < 2000; ++j)
   {

       icheck = fscanf(pFile[0], "%lf %lf", &x, &y);

       if(  icheck <= 0 )
       {

			printf("\n %d  %lf %lf  %d", j, x, y, icheck);

			printf("\n  data error ");

			break;

       }               
       if(j==0)
       {
			x1=x;
       }
   }


   if(j<=1)
   {
		printf("\n data error j=%d", j);
		exit(1);
   }
   else
   {
		printf(" \n j = %ld \n",j);
   }

   dt=(x-x1)/(j-1);

   sr=(1./dt);

   limit=long(sr*tlimit);

   printf("\n dt= %f  sr= %f ", dt, sr);

//   printf("\n time segment = %f  samples per segment = %ld\n", tlimit, limit);


//   printf("\n\n begin rewind \n");

   fseek(pFile[0], 0L, SEEK_SET);     /*  rewind */

//   printf("\n\n end rewind \n");

   i=0L;
   nzc=0L;

//   printf("\n\n begin STAT scan \n");

    ybefore=0.;

	while( fscanf(pFile[0], "%lf %lf", &x, &y) > 0 )
	{
           
           
        if(i>1)
        {
             if((y*ybefore)<0.)
             { 
                 nzc++;
             }  
        }  
    
           
		i+=1;

		STAT(y, sum, sum2);

		if(fabs(y)>peak){peak=fabs(y);}

		ybefore=y;
		
		if(nzc==20)
		{
			printfile();
		}
   }

//   printf("\n i=%ld limit=%ld ",i,limit);  
   
   printf("\n\n Output files: \n");


   for(i=5;i<=8;i++)
   {
      printf("  %s \n",filename[i]);
   }

   freq_amp();


   for(i=0;i<=8;i++)
   {
      fclose(pFile[i]);
   }

	printf("\n\n Please call the output files into your own plotting program.");

	printf("\n Press any key to exit.");
	getch();

}
void printfile(void)
{
		double xp;

		if( i >= 1)
		{
			ave=sum/i;

			grms=sqrt(sum2/i);


			if( (pow(grms,2) - pow(ave,2) ) >=0.)
			{		  
				sd = sqrt( pow(grms,2) - pow(ave,2) );
			} 
			else
			{
				sd=0.;
			}
		}
		else
		{
			printf("\n i error \n");
		}

		xp=(x+x1)/2.;

//		printf("\n samples= %ld ", i);  

        freq=(nzc/2.)/(fabs(x-x1));





// printf("\n %12.6e \t %12.6e \t %12.6e \t %12.6e \t %12.6e ", xp, ave, sd, grms, maxabs);

fprintf(pFile[1], "%12.6e \t  %12.6e \t %12.6e \t %12.6e \t %12.6e \n",xp, ave, sd, grms, maxabs );



		fprintf(pFile[2], "%12.6e \t %12.6e \n", xp, sd );
		fprintf(pFile[3], "%12.6e \t %12.6e \n", xp, ave );
		fprintf(pFile[4], "%12.6e \t %12.6e \n", xp, peak );
		
		double peak2=sd*sqrt(2);
		
		if(peak>peak2)
		{peak2=peak;}
		
		fprintf(pFile[5], "%12.6e \t %12.6e \n", xp, freq );
		fprintf(pFile[6], "%12.6e \t %12.6e \n", xp, peak2 );
//		fprintf(pFile[7], "%12.6e \t %12.6e \n", freq, peak2 );		
		fprintf(pFile[8], "%12.6e \t %12.6e \t %12.6e\n", xp, freq,peak2 );	
		
        ff[nv]=freq;
        pp[nv]=peak2;
        nv++;


		peak=0.;
		maxabs=0.;
		sum=0.;
		sum2=0.;

		i=0;
		
		nzc=0;

		x1=x;
}
void freq_amp()
{
        long last=nv-1;
   
   double temp;
   
   printf("\n last=%ld \n",last);
   
   for(i=0;i<=(last-1);i++)
   {
      for(j=i+1;j<=last;j++)
      {
           if(ff[i]>ff[j])
           {
               temp=ff[i];
               ff[i]=ff[j];
               ff[j]=temp;
               
               temp=pp[i];
               pp[i]=pp[j];
               pp[j]=temp;                                       
          }
      }                    
   }
   
   i=1;
   fprintf(pFile[7], "%12.6e \t %12.6e \n", ff[i],pp[i] );	
         
   for(i=0;i<=last;i++)
   {
      if(ff[i]>ff[i-1])
      {                 
      fprintf(pFile[7], "%12.6e \t %12.6e \n", ff[i],pp[i] );	                       
      }
   }

}