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


void read_data(void);
void psd_sum(void);
void psd_ave(void);
void freq_check(void);


FILE *pFile[5];
char filename[5][25];

double freq[400][2000];
double  amp[400][2000];

long num[400];

double df[400];

double psd[2000];


long i,j,jmax,k;

long imax;


void main()
{

	for(i=0; i<2000; i++)
	{
         psd[i]=0.;
	}

	printf("\n\n   psd_ave.cpp,  ver 1.0,  by Tom Irvine \n");


	printf("\n\n   This program calculates the average 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();


	strcpy(filename[0],"files.in");
    pFile[0]=fopen(filename[0],"rb");


   if(pFile[0] == NULL )
   {
      printf(" Failed to open file: %s \n", filename[0]);
      fclose(pFile[0]);
      exit(1);
   }


	strcpy(filename[2],"file_check.txt");
    pFile[2]=fopen(filename[2],"w");

	
	strcpy(filename[3],"average_psd.txt");
    pFile[3]=fopen(filename[3],"w");



	read_data();

	rewind(pFile[0]);


    freq_check();

	psd_sum();

	psd_ave();



	printf("\n\n Calculation complete. \n\n");


	printf("\n\n   Please open %s to check file status. \n",filename[2]);

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

	
	printf("\n   The output file is: average_psd.txt \n");
	
	printf("\n\n Press any key to exit.");
	

	getch();

}
void read_data(void)
{
	i=0;

    while(1)
	{
   
		if( fscanf(pFile[0],"%s",filename[1]) <= 0 ){break;}
        
        pFile[1]=fopen(filename[1],"rb");




		if(pFile[1] == NULL )
		{
			printf(" Failed to open file: %s \n", filename[1]);
			fclose(pFile[1]);
			exit(1);
		}





		j=0;

        while(1)
		{
            if( fscanf(pFile[1],"%lf %lf",&freq[i][j],&amp[i][j] ) <= 0 ){break;}

            df[i]=( freq[i][j]-freq[i][0] );

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

            j++;
        }

		num[i]=j;
         df[i]/=double(j);

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


		if( j > jmax){jmax=j;}

		fclose(pFile[1]);

		i++;
	}
//	getch();

	imax=i;
}

void freq_check(void)
{



	fprintf(pFile[2],"\n  Filename         df(Hz)   number of points \n\n");	


	for(i=0; i<imax; i++)
	{
		
        fscanf(pFile[0],"%s",filename[1]);

		fprintf(pFile[2],"   %s      %8.3lf       %ld  \n",filename[1],df[i],num[i]);
	
	}


    fclose(pFile[2]);


	for(k=0; k<jmax; k++)
	{
	
		for(i=0; i<imax-1; i++)
		{

			for(j=i+1; j<imax; j++)
			{

				if(  fabs( freq[i][k]-freq[j][k] ) > 0.01 )
				{
					printf(" Frequency synchronization error \n");

					exit(1);
				}

            }
		}
	}

}

void psd_sum(void)
{

	for(i=0; i<imax; i++)
	{	
		for(j=0; j< jmax; j++)
		{

			psd[j]+=amp[i][j];
 
		}
	}
}
void psd_ave(void)
{

	for(j=0; j< jmax; j++)
	{
		psd[j]/=double(imax);

       fprintf(pFile[3]," %12.4g %12.4g \n",freq[0][j],psd[j] );

	}


}