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


void onethird(void);

void integrate_input(void);

void integrate_output(void);

void apply_weights(void);


void dataentry(void);
void read_data(void);


#define MAX 100000

#define LN 400

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

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

long i,iamp,iunits,imax,j,k,num;
long iweight;

long nn;


double fc[50],wk[50],wd[50],wf[50],wc[50],we[50],wj[50],wb[50],ww[50];

double s[MAX];

double f[MAX],amp[MAX];

double first,last;


int numBytes=320;

char string[320];

void main()
{

	dataentry();

	
	printf("\n reading input data\n");

	read_data();



    printf("\n calculate weights \n");   
	
	onethird();



    printf("\n integrate input PSD \n");
 
	integrate_input();



    printf("\n apply weights \n");

	apply_weights();
  


    printf("\n integrate output PSD \n");

	integrate_output();




    fclose(pFile[1]);

	printf("\n\n Output File: %s  ",filename[1]);

	printf("\n\n Format: freq (Hz)   Weighted PSD ((m/sec^2)^2)/Hz)     \n");




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

}
void dataentry(void)
{
   printf("\n  psd_weighted - version 1.0, December 28, 2004 ");
   printf("\n  by Tom Irvine ");
   printf("\n  Email;  tomirvine@aol.com \n");

   printf("\n This program converts a PSD to a weighted format per ISO 2631. \n");
   printf("\n The PSD may have either a constant bandwidth or one-third octave format. \n");

   printf("\n Enter filename for input PSD \n");
   scanf("%s",&filename[0]);
   pFile[0]=fopen(filename[0], "rb");


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


   iunits = 0;

   while( iunits != 1 && iunits !=2 )
   {
		printf("\n\n Enter input units;\n  1=G^2/Hz  2= (m/sec^2)^2/Hz \n");
		scanf("%ld",&iunits);
   }
  
   iweight = 0;

   while( iweight != 1 &&  iweight != 1 &&   iweight != 2 &&   iweight != 3 &&   iweight != 4 &&  iweight != 5 &&  iweight != 6   &&  iweight != 7  )
   {
		printf("\n\n Enter weighting by choice number: \n");
		printf("\n  1= Wk  z-axis,      seat surface"); 
		printf("\n  2= Wd  x or y-axis, seat surface");
		printf("\n  3= Wf"); 				   
		printf("\n  4= Wc");  
		printf("\n  5= We"); 
		printf("\n  6= Wj");
		printf("\n  7= Wb    \n\n");

		scanf("%ld",&iweight);
   }
  


   printf("\n Enter filename for output PSD \n");
   scanf("%s",&filename[1]);
   pFile[1]=fopen(filename[1], "w");


   strcpy(filename[2],"wint.out");
   pFile[2]=fopen(filename[2], "w");


}
void read_data(void)
{
/*
	for(i=0; i< num; i++)
	{
		fgets(string,numBytes,pFile[0] );
	}
*/

	j=0;

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

//		printf("\n  %lf %lf", f[j], amp[j] );

		if(iunits==1)
		{ 
			amp[j]*=pow(9.81,2.); 
		}
        			
		j++;

		if(j==MAX){ break;}
	}
	nn=j;

	fclose(pFile[0]);

}
void onethird(void)
{

	fc[0]=0.02; 
	wk[0]=-90.; 
	wd[0]=-90.; 
	wf[0]=-32.33; 
	wc[0]=-90.; 
	we[0]=-90.; 
	wj[0]=-90.; 
	wb[0]=-90.; 
	
	
	fc[1]=0.025; 
	wk[1]=-90.; 
	wd[1]=-90.; 
	wf[1]=-28.48;
	wc[1]=-90.;  
	we[1]=-90.;  
	wj[1]=-90.;  
	wb[1]=-90.; 
	
	
	fc[2]=0.0315;
	wk[2]=-90.; 
	wd[2]=-90.; 
	wf[2]=-24.47; 
	wc[2]=-90.;  
	we[2]=-90.;  
	wj[2]= -90.;  
	wb[2]=-90.; 
	
	
	fc[3]=0.04;
	wk[3]=-90.; 
	wd[3]=-90.; 
	wf[3]=-20.25; 
	wc[3]=-90.;  
	we[3]=-90.;  
	wj[3]=-90.;  
	wb[3]=-90.; 
	
	
	fc[4]=0.05;
	wk[4]=-90.; 
	wd[4]=-90.; 
	wf[4]=-16.10; 
	wc[4]=-90.;  
	we[4]=-90.;  
	wj[4]=-90.;  
	wb[4]=-90.; 
	
	
	fc[5]=0.063;
	wk[5]=-90.; 
	wd[5]=-90.; 
	wf[5]=-11.49; 
	wc[5]=-90.;  
	we[5]= -90.;  
	wj[5]=-90.;  
	wb[5]=-90.; 
	
	
	fc[6]=0.08;
	wk[6]=-90.; 
	wd[6]=-90.; 
	wf[6]=-6.73; 
	wc[6]= -90.;  
	we[6]=-90.;  
	wj[6]=-90.;  
	wb[6]=-90.; 
	
	
	fc[7]=0.1;
	wk[7]=-30.11; 
	wd[7]=-24.09; 
	wf[7]=-3.16; 
	wc[7]=-24.11;  
	we[7]=-24.08;  
	wj[7]=-30.18;  
	wb[7]=-32.04; 
	
	
	fc[8]=0.125;
	wk[8]=-26.26; 
	wd[8]=-20.24; 
	wf[8]=-0.96; 
	wc[8]=-20.25;  
	we[8]=-20.22;  
	wj[8]=-26.32;  
	wb[8]=-28.20; 
	
	
	fc[9]=0.16;
	wk[9]=-22.05; 
	wd[9]=-16.01; 
	wf[9]=0.05; 
	wc[9]=-16.03; 
	we[9]=-15.98;  
	wj[9]=-22.11;  
	wb[9]=-23.98; 
	
	
	fc[10]=0.2;
	wk[10]=-18.33; 
	wd[10]=-12.28; 
	wf[10]=-0.07; 
	wc[10]=-12.03;  
	we[10]=-12.23;  
	wj[10]=-18.38;  
	wb[10]=-20.23; 
	
	
	fc[11]=0.25;
	wk[11]=-14.81; 
	wd[11]=-8.75; 
	wf[11]=-1.37; 
	wc[11]=-8.78; 
	we[11]=-8.67;  
	wj[11]=-14.86;  
	wb[11]=-16.71; 
	
	
	fc[12]=0.315;
	wk[12]=-11.60; 
	wd[12]=-5.52; 
	wf[12]=-4.17; 
	wc[12]=-5.56;  
	we[12]=-5.41;  
	wj[12]=-11.65;  
	wb[12]=-13.51; 
	
	
	fc[13]=0.4;
	wk[13]=-9.07; 
	wd[13]=-2.94; 
	wf[13]=-8.31; 
	wc[13]=-3.01;  
	we[13]=-3.01;  
	wj[13]=-9.10;  
	wb[13]=-10.98; 
	
	
	fc[14]=0.5;
	wk[14]=-7.57; 
	wd[14]=-1.38; 
	wf[14]=-13.00; 
	wc[14]=-1.48;  
	we[14]=-1.29;  
	wj[14]=-7.60;  
	wb[14]=-9.53; 
	
	
	fc[15]=0.63;
	wk[15]=-6.77; 
	wd[15]=-0.50; 
	wf[15]=-18.69; 
	wc[15]=-0.64; 
	we[15]=-0.55;  
	wj[15]=-6.78;  
	wb[15]=-8.71; 
	
	
	fc[16]=0.8;
	wk[16]=-6.43; 
	wd[16]=-0.07; 
	wf[16]=-25.51; 
	wc[16]=-0.24;  
	we[16]=-0.53;  
	wj[16]=-6.42;  
	wb[16]=-8.38; 
	
	
	fc[17]=1.;
	wk[17]=-6.33; 
	wd[17]=0.10; 
	wf[17]=-32.57; 
	wc[17]=-0.08;  
	we[17]=-1.11;  
	wj[17]=-6.30;  
	wb[17]=-8.29; 
	
	
	fc[18]=1.25;
	wk[18]=-6.29; 
	wd[18]=0.07; 
	wf[18]=-40.02; 
	wc[18]=0.00; 
	we[18]=-2.25;  
	wj[18]=-6.28;  
	wb[18]=-8.27; 
	
	
	fc[19]=1.6;
	wk[19]=-6.12; 
	wd[19]=-0.28; 
	wf[19]=-48.47; 
	wc[19]=0.06;  
	we[19]=-3.99;  
	wj[19]=-6.32;  
	wb[19]=-8.07; 
	
	
	fc[20]=2.;
	wk[20]=-5.49; 
	wd[20]=-1.01; 
	wf[20]=-56.19; 
	wc[20]=0.10;  
	we[20]=-5.82;  
	wj[20]=-6.34;  
	wb[20]=-7.60; 
	
	
	fc[21]=2.5;
	wk[21]=-4.01; 
	wd[21]=-2.20; 
	wf[21]=-63.93; 
	wc[21]=0.15;  
	we[21]=-7.77;  
	wj[21]=-6.22;  
	wb[21]=-6.13; 
	
	
	fc[22]=3.15;
	wk[22]=-1.90; 
	wd[22]=-3.85; 
	wf[22]=-71.96; 
	wc[22]=0.19;  
	we[22]=-9.81;  
	wj[22]=-5.62;  
	wb[22]=-3.58; 
	
	
	fc[23]=4.;
	wk[23]=-0.29; 
	wd[23]=-5.82; 
	wf[23]=-80.26; 
	wc[23]=0.20;  
	we[23]=-11.93;  
	wj[23]=-4.04;  
	wb[23]=-1.02; 
	
	
	fc[24]=5.;
	wk[24]=0.33;  
	wd[24]=-7.76;  
	wf[24]=-90.; 
	wc[24]=0.11;  
	we[24]=-13.91;  
	wj[24]=-2.01;  
	wb[24]=0.21; 
	
	
	fc[25]=6.3;
	wk[25]=0.46;  
	wd[25]=-9.81;  
	wf[25]=-90.;  
	wc[25]=-0.23;  
	we[25]=-15.94;  
	wj[25]=-0.48;  
	wb[25]=0.46; 
	
	
	fc[26]=8.;
	wk[26]=0.31;  
	wd[26]=-11.93;  
	wf[26]=-90.;  
	wc[26]=-1.00;  
	we[26]=-18.03;  
	wj[26]=0.15;  
	wb[26]=0.21; 
	
	
	fc[27]=10.;
	wk[27]=-0.10;  
	wd[27]=-13.91;  
	wf[27]=-90.;  
	wc[27]=-2.20;  
	we[27]=-19.98;  
	wj[27]=0.26;  
	wb[27]=-0.23; 
	
	
	fc[28]=12.5;
	wk[28]=-0.89;  
	wd[28]=-15.87;  
	wf[28]=-90.;  
	wc[28]=-3.79;  
	we[28]=-21.93;  
	wj[28]=0.22;  
	wb[28]=-0.85; 
	
	
	fc[29]=16.;
	wk[29]=-2.28;  
	wd[29]=-18.03;  
	wf[29]=-90.;  
	wc[29]=-5.82;  
	we[29]=-24.08;  
	wj[29]=0.16;  
	wb[29]=-1.83; 
	
	
	fc[30]=20.;
	wk[30]=-3.93;  
	wd[30]=-19.99;  
	wf[30]=-90.;  
	wc[30]=-7.77;  
	we[30]=-26.02;  
	wj[30]=0.10;  
	wb[30]=-3.00; 
	
	
	fc[31]=25.;
	wk[31]=-5.80;  
	wd[31]=-21.94;  
	wf[31]=-90.;  
	wc[31]=-9.76;  
	we[31]=-27.97;  
	wj[31]=0.06;  
	wb[31]=-4.44; 
	
	
	fc[32]=31.5;
	wk[32]=-7.86;  
	wd[32]=-23.98;  
	wf[32]=-90.;  
	wc[32]=-11.84;  
	we[32]=-30.01;  
	wj[32]=0.00;  
	wb[32]=-6.16; 
	
	
	fc[33]=40.;
	wk[33]=-10.05;  
	wd[33]=-26.13;  
	wf[33]=-90.;  
	wc[33]=-14.02;  
	we[33]=-32.15;  
	wj[33]=-0.08;  
	wb[33]=-8.11; 
	
	
	fc[34]=50.;
	wk[34]=-12.19;  
	wd[34]=-28.22;  
	wf[34]=-90.;  
	wc[34]=-16.13;  
	we[34]=-34.24;  
	wj[34]=-0.24;  
	wb[34]=-10.09; 

	
	
	fc[35]=63.;
	wk[35]=-14.61;  
	wd[35]=-30.60;  
	wf[35]=-90.;  
	wc[35]=-18.53;  
	we[35]=-36.62;  
	wj[35]=-0.62;  
	wb[35]=-12.43; 
	
	
	fc[36]=80.;
	wk[36]=-17.56;  
	wd[36]=-33.53;  
	wf[36]=-90.;  
	wc[36]=-21.47;  
	we[36]=-39.55;  
	wj[36]=-1.48;  
	wb[36]=-15.34; 
	
	
	fc[37]=100.;
	wk[37]=-21.04;  
	wd[37]=-36.99;  
	wf[37]=-90.;  
	wc[37]=-24.94;  
	we[37]=-43.01;  
	wj[37]=-3.01;  
	wb[37]=-18.72; 
	
	
	fc[38]=125.;
	wk[38]=-25.35;  
	wd[38]=-41.28;  
	wf[38]=-90.;  
	wc[38]=-29.24;  
	we[38]=-47.31;  
	wj[38]=-5.36;  
	wb[38]=-23.00; 
	
	
	fc[39]=160.;
	wk[39]=-30.91;  
	wd[39]=-46.84;  
	wf[39]=-90.;  
	wc[39]=-34.80;  
	we[39]=-52.86;  
	wj[39]=-8.78;  
	wb[39]=-28.56; 
	
	
	fc[40]=200.;
	wk[40]=-36.38;  
	wd[40]=-52.30;  
	wf[40]=-90.;  
	wc[40]=-40.26;  
	we[40]=-58.33;  
	wj[40]=-12.30;  
	wb[40]=-34.03; 
	
	
	fc[41]=250.;
	wk[41]=-42.0;  
	wd[41]=-57.97;  
	wf[41]=-90.;  
	wc[41]=-45.92;  
	we[41]=-63.99;  
	wj[41]=-16.03;  
	wb[41]=-39.69; 
	
	
	fc[42]=315.;
	wk[42]=-48.00;  
	wd[42]=-63.92;  
	wf[42]=-90.;  
	wc[42]=-51.88;  
	we[42]=-69.94;  
	wj[42]=-19.98;  
	wb[42]=-46.65; 
	
	
	fc[43]=400.;
	wk[43]=-54.20;  
	wd[43]=-70.12;  
	wf[43]=-90.;  
	wc[43]=-58.08;  
	we[43]=-76.14;  
	wj[43]=-24.10;  
	wb[43]=-51.84; 

   
    strcpy(filename[10],"wk.out");
    strcpy(filename[11],"wd.out");
    strcpy(filename[12],"wf.out");
    strcpy(filename[13],"wc.out");
    strcpy(filename[14],"we.out");
    strcpy(filename[15],"wj.out");
    strcpy(filename[16],"wb.out");


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

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

		fprintf(pFile[10]," %12.4lf  %12.4lf \n",fc[i],wk[i]);
		fprintf(pFile[11]," %12.4lf  %12.4lf \n",fc[i],wd[i]);
		fprintf(pFile[12]," %12.4lf  %12.4lf \n",fc[i],wf[i]);
		fprintf(pFile[13]," %12.4lf  %12.4lf \n",fc[i],wc[i]);
		fprintf(pFile[14]," %12.4lf  %12.4lf \n",fc[i],we[i]);
		fprintf(pFile[15]," %12.4lf  %12.4lf \n",fc[i],wj[i]);
		fprintf(pFile[16]," %12.4lf  %12.4lf \n",fc[i],wb[i]);

	}

	for(i=0; i<=43; i++)  // for PSD
	{
		wk[i]=pow(10.,(wk[i]/10.));
		wd[i]=pow(10.,(wd[i]/10.));
		wf[i]=pow(10.,(wf[i]/10.));
		wc[i]=pow(10.,(wc[i]/10.));
		we[i]=pow(10.,(we[i]/10.));
		wj[i]=pow(10.,(wj[i]/10.));
		wb[i]=pow(10.,(wb[i]/10.));
	}

	for(i=0; i<=43; i++)  // for PSD
	{

		if(iweight==1)
		{ 
			ww[i]=wk[i]; 
		}
		if(iweight==2)
		{ 
			ww[i]=wd[i]; 
		}
		if(iweight==3)
		{ 
			ww[i]=wf[i]; 
		}
		if(iweight==4)
		{ 
			ww[i]=wc[i]; 
		}
		if(iweight==5)
		{ 
			ww[i]=we[i]; 
		}
		if(iweight==6)
		{ 
			ww[i]=wj[i]; 
		}
		if(iweight==7)
		{ 
			ww[i]=wb[i]; 
		}

	}
//
	   first = 0.1;

	   for(i=0; i<43; i++)
	   {
			if(ww[i] > 1.0e-8)
			{
				first=fc[i];
				break;
			}
	   }

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

//		    printf(" %12.4g  %12.4g \n",fc[i],ww[i]);

			if(ww[i] > 1.0e-8)
			{
				last=fc[i];
			}
	   }

//	   getch();


//  Calculate slopes between weights

	double a,b;

       for(i=0; i< 43; i++)
	   {
          a=log(ww[i+1])-log(ww[i]);
          b=log(fc[i+1])-log(fc[i]);
          s[i]=a/b;

//		  printf("%lf\n",s[i]);
       }




}
void integrate_input(void)
{
	
	double slope[MAX];

	double ra,grms;

	ra=0.;
	grms=0.;

	for(i=0; i < MAX; i++)
	{
		slope[i]=0.;
	}

	printf("\n number of input points = %ld \n",nn);
		  

//  calculate slopes

    if(f[0] < 1.0e-04){ f[0] = 1.0e-04; }

	for( i=0; i < nn-1; i++)
	{
		if( amp[i] < 1.0e-60 ||  amp[i+1] < 1.0e-60 )
		{
			printf("\n input amplitude error \n");
			exit(1);
		}
		if( f[i] < 1.0e-60 ||  f[i+1] < 1.0e-60 )
		{
			printf("\n input frequency error \n");
			exit(1);
		}



		      s[i]=log( amp[i+1]/amp[i] )/log( f[i+1]/f[i] );

              if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			  {	   
                 ra+= ( amp[i+1] * f[i+1]- amp[i]*f[i])/( s[i]+1.);
			  }
		      else
			  {
                 ra+= amp[i]*f[i]*log( f[i+1]/f[i]);
			  }

	}
    grms=sqrt(ra);

	printf("\n Overall Input Level = %12.4g (m/sec^2)RMS ",grms);
	
	printf("\n                     = %12.4g GRMS\n\n",(grms/9.81));


}
void integrate_output(void)
{
	
	double slope[MAX];

	double ra,grms;

	ra=0.;
	grms=0.;

	for(i=0; i < MAX; i++)
	{
		slope[i]=0.;
	}
		  

//  calculate slopes

	for( i=0; i < nn-1; i++)
	{


	    if( amp[i] < 1.0e-60 ||  amp[i+1] < 1.0e-60 )
		{
			printf("\n input amplitude error \n");
			exit(1);
		}
		if( f[i] < 1.0e-60 ||  f[i+1] < 1.0e-60 )
		{
			printf("\n input frequency error \n");
			exit(1);
		}


		if( f[i] >= first && f[i+1] <= last )
		{

			s[i]=log( amp[i+1]/amp[i] )/log( f[i+1]/f[i] );

			if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			{	   
                 ra+= ( amp[i+1] * f[i+1]- amp[i]*f[i])/( s[i]+1.);
			}
			else
			{
                 ra+= amp[i]*f[i]*log( f[i+1]/f[i]);
			}
		}
	}
    grms=sqrt(ra);

	printf("\n Overall Weighted Output Level = %12.4g (m/sec^2)RMS \n",grms);
	

}
void apply_weights(void)
{
	int iflag=0;

	double ss;

	double ratio, rr, wint;

	for(i=0; i<nn; i++)
	{
        if( f[i] > last ){break;}

		iflag=0;

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

            rr = (f[i] - fc[j])/fc[j];

			if( fabs( rr )  < 0.01)
			{
				amp[i]*=ww[j];

				fprintf(pFile[2]," %12.4g  %12.4g \n",f[i],ww[j]);

				iflag=1;
				break;
			}

			if( j < 43 )
			{
				if(   f[i] > fc[j] && f[i] < fc[j+1] )
				{

					ss = log(ww[j+1]/ww[j])/log(fc[j+1]/fc[j]);  

					ratio = f[i]/fc[j];

					wint= ww[j]*pow(ratio,ss); 

			        fprintf(pFile[2]," %12.4g  %12.4g \n",f[i],wint);

                    if( wint > ww[j] && wint > ww[j+1] )
					{

						printf("\n error:\n f=%8.4g  fc=%8.4g  ss=%12.4g  wint=%12.4g \n",f[i],fc[j],ss,wint);
					
					    exit(1);
					}


					amp[i]*=wint;

					iflag=1;
					break;
				}
            }

		}

		if( iflag==1 && f[i]>=first && f[i]<=last)
		{
			fprintf(pFile[1],"\n %10.2lf \t %12.4e",f[i],amp[i]);
//			 printf("\n %10.2lf \t %12.4e",f[i],amp[i]);
		}
	
	}

//	printf("\n Start Freq =%12.4g Hz  End Freq =%12.4g Hz \n",first,last);

}