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



#define MAX 2500000



#define NUM 200



#define LIMIT 18000



void read_filenames(void);

void find_limits(void);

void read_data(void);

void print_file(void);



double freq[NUM][LIMIT],psd[NUM][LIMIT];

double max[NUM];



double f,x;



double ffmin,ffmax;



double ff,aa;



long i,j,k,klast[NUM],kmin,last;



long ke;


char string[50];


char filename[NUM][FILENAME_MAX];
FILE *pFile[NUM];    /*  array of pointers  */


int main(  )
{

	printf("\n\n   psd_max.cpp,  ver 1.8,  November 9, 2010\n ");

	   printf("\n   by Tom Irvine  ");

	   printf("\n   Email: tomirvine@aol.com \n\n");





	printf("\n\n   This program calculates the maximum envelope of a number of PSD functions. \n");

	  printf("\n   The names of the individual files must be in a file called files.in. \n");



	  printf("\n   Each PSD file must have the same number of lines and the same \n   frequency step. \n");	





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

	getch();



     printf("\n ref 1 ");

 

	 for(j=0; j< 500; j++)

	 {

		 

			

		 for(i=0; i< NUM; i++)
		 {

			 freq[i][j]=0.;

			  psd[i][j]=0.;

		 }

	 }
	 
	 for(i=0; i< NUM; i++)
     {	 
	      max[i] = 0.;
     }

     strcpy(filename[0],"files.in");

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


	 j=0;



	 if(pFile[j] == NULL )
	 {

			printf(" Failed to open file: \n %s \n", filename[j]);

			exit(1);

	 }

	 else

	 {

			printf(" Opened File: \n %s \n", filename[j]);

	 }



     printf("\n ref 2 \n");





	 read_filenames();





     printf("\n ref 3 ");



	 printf("\n\n last = %ld  \n",last);





     find_limits();



	 read_data();







	 printf("\n ref 4 ");



     kmin=klast[1];



     for(j=2; j<=last; j++)
	 {

		 if(kmin > klast[j] )
		 { 

			 kmin=klast[j];



			 printf(" j=%ld  kmin=%ld  \n",j,kmin);

		 }

	 }

	 printf("\n 1  kmin = %ld \n",kmin);





     printf("\n ref 5 \n\n");





     ke = kmin;





	 for(j=1; j<=last; j++)
	 {



		 for(i=0; i<=ke; i++)
		 {



//            printf(" psd[%ld][%ld]=%g \n",j,i,psd[j][i]);



			if(psd[j][i] > max[i] )
			{

				max[i] = psd[j][i];

			}

		 

		 }



	 }





     printf("\n ref 6 \n");



//     printf("\n 2  ke = %ld \n",ke);

	 

	 for(i = 0; i <= ke; i++ )
	 {

		 for(j=2; j<=last; j++)
		 {



			if(  fabs( freq[j][i] - freq[j-1][i] )  > 0.01*freq[j][i] )
			{

              

				 printf(" i=%ld \t j=%ld \t %12.4g \t %12.4g \n",i,j,freq[j][i],freq[j-1][i]);



                 printf("\n Synch Error \n");

				 exit(1);

			}

		 

		 }

     }

//     printf("\n 3  ke = %ld \n",ke);





    print_file();



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

	



	getch();



 
}


void find_limits(void)
{

	ffmin=0.;

	ffmax=1.0e+99;



	i=1;



    long ijk;



	printf("\n find_limits \n");



	for(ijk=1; ijk<=last; ijk++)
	{

		double local_ffmax=0.;

   

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



		if(pFile[1] == NULL )
		{

			printf("\n  Failed to open file: %s \n", filename[i]);

	

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

			getch();



			exit(1);

		}

		else
		{

//		     printf("\n open file: %s \n",filename[i]);

		}



		j=0;



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

	

        while(1)

		{

            if( fscanf(pFile[1],"%lf %lf",&ff,&aa ) <= 0 ){break;}



//           printf(" ff=%12.4g  aa=%12.4g \n",ff,aa);



			local_ffmax=ff;



			if( j == 0 )

			{ 

				if(ff > ffmin ){ ffmin = ff;}

			}

            j++;



//			if(j==NUM){break;}



        }

		fclose(pFile[1]);

	

		printf("\n %s \tffmin=%8.4g Hz\tffmax=%8.4g Hz\tj=%ld",filename[i],ffmin,local_ffmax,j);



        if(ff<ffmax){ffmax=ff;}



		i++;

	}

	printf("\n\n composite ffmax=%8.4g \n",ffmax);

}

void read_data(void)

{



	 printf("\n read_data \n");



  

	 for(j=1; j<= last; j++)

	 {

		 k=0;



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

		 



		if(pFile[j] == NULL )

		{

			printf(" Failed to open file: \n %s \n", filename[j]);

			exit(1);

		}

		else

		{

			printf(" Opened File:  %s \n", filename[j]);

		}







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





         while(  fscanf(pFile[j],"%lf %lf",&f,&x) > 0  )  

		 {

		 

//            printf(" f=%g  x=%g \n",f,x);


			 if( f >= ffmin && f <=ffmax)

			 {



				freq[j][k]=f;

				 psd[j][k]=x;





				klast[j]=k;





				k++;

			 }

			 if(k==LIMIT){break;}

		 }

		 fclose(pFile[j] );



//		 printf(" j=%ld  klast=%ld  \n",j,klast[j]);

	 	 

	 }

}

void read_filenames(void)
{



     i=1;

	 

     while( fscanf(pFile[0],"%s",string) > 0  )

	 {



		 strcpy(filename[i],string);



		 printf(" %s \n",string);

			 

         

		 last = i;

	 

		 if(i==NUM){break;}



	     i++;

	 }

	 fclose(pFile[0]);



}

void print_file(void)

{

	

     printf("\n ref 7 ");



	 double rms=0.;

     double ss;



     strcpy(filename[0],"psdmax.out");

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





     strcpy(filename[1],"psdmax.rms");

     pFile[1]=fopen(filename[1],"a");





//     printf("\n 4   ke = %ld \n",ke);



	 for(i=0; i<=ke; i++)

	 {

			fprintf(pFile[0]," %12.4g \t %12.4g \n",freq[1][i],max[i] );

	 }



	 for(i=0; i<ke; i++)

	 {

			if(freq[1][i]<=0.)

			{

				printf("warning: freq[1][%ld]=%8.4g  \n",i,freq[1][i]);

			}

			if(max[i]<=0.)

			{

				printf("warning: max[%ld]=%8.4g      \n",i,max[i]);

			}





			if(i<ke)

			{

				ss=log( max[i+1]/max[i] )/log( freq[1][i+1]/freq[1][i] );

			}



			if(ss < -1.0001 ||  ss > -0.9999 )

			{	   

                 rms+= ( max[i+1]*freq[1][i+1]- max[i]*freq[1][i])/( ss+1.);

			}

			else

			{

                 rms+= max[i]*freq[1][i]*log( freq[1][i+1]/freq[1][i]);

			}



			if(rms>=0 && rms <=1.0e+90)

			{

			}

			else

			{

				printf("\n warning  freq[1][%ld]=%8.4g   ss=%8.4g  rms= %8.4g \n",i,freq[1][i],ss,rms);

			}



 //           printf(" %12.4g \t %12.4g \n",freq[1][i],rms );



	 }



	 rms=sqrt(rms);







	 fprintf(pFile[1],"\n\n overall level =  %8.2lf \n",rms );

	 

	 fclose(pFile[1]);





//	 printf("\n ref 8 ");



	 printf("\n\n overall level = %8.2lf \n",rms );

	 printf("\n   The output file is: psdmax.out \n");

}