/*
     Vibration response spectrum   
*/


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


#define MAX 20000


void fileoper(void);

void calc(double *level);

void vrs(void);

void interp(void);
void header(void);


double duration;

double level[MAX];

double f[MAX],fi[MAX], a[MAX],s[MAX],v[MAX],d[MAX],ai[MAX];

double df, db, pi,tpi,ra,grms;

double oa, ov, od;

double fstart,flast;

double f_peak=0.;
double f_freq=0.;


int in[5];
int iflag;
int i,ii, ic, imethod;
int j;
int last;
int m, n, nn;


char amp_label[5][20];

long numBytes=300;
char string[300];


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


void main()
{

	   cout.precision(3);

	   iflag=0;

	   pi=atan2(0.,-1.);

       tpi=2.*pi;


	   strcpy(amp_label[1],"acceleration");
	   strcpy(amp_label[2],"velocity");
	   strcpy(amp_label[3],"displacement");

	   in[1]=1;
	   in[2]=2;
	   in[3]=3;


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

		   f[i]=0.;
          fi[i]=0.;

           a[i]=0.;
		  ai[i]=0.;

           s[i]=0.;

           v[i]=0.;
           d[i]=0.;

	   }

	   oa=0.;
	   ov=0.;
	   od=0.;


       cout << "                               \n";
       cout << " VRS.cpp  ver 5.0, May 16, 2008    \n";
       cout << "							   \n";
       cout << " by Tom Irvine                 \n";

       cout << "    Email:  tomirvine@aol.com  \n";
	   cout << "							   \n";


       cout << "This program calculates the acceleration vibration response \n"; 
       cout << "spectrum of a single-degree-of-freedom system to a          \n";
       cout << "power spectral density input to the base of the system.   \n\n";



       cout << "The power spectral density may vary arbitrarily with frequency. \n";
       cout << "Natural frequency is an independent variable.  Damping is fixed.\n\n";



	   ic=1;

	   while(ic==1)
	   {

          cout << " Select data input method \n";

          cout << " 1=keyboard  2=file   3=library PSD     \n";

	      cin  >> imethod;

       

		  if(imethod==1)
		  {    

              cout << " How Many Breakpoints (minimum=2)? \n";
          	  cin  >> n;


              if(n < 2)
			  {
				  cout << "Input error \n";
			  }
          
              cout << " Enter Points (free-format) "<< endl;
              cout << " Freq(Hz)  Level(G^2/Hz)    "<< endl;

          
              for( i=0; i<n; i++)
			  {
				  scanf("%lf %lf",&f[i],&a[i]);
			  }

          }

          if(imethod==2)
		  {

			 printf( " The file must contain two columns: \n");
             printf( " Freq(Hz)  Level(G^2/Hz)    \n");

			 printf( " \n An error message will occur if the file does not exist.\n");
             printf( " Input filename \n");
		

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

     		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");

			}

			printf(" File: %s opened. \n", filename[0]);

			header();

			 i=0;


             double fff,aaa;


//			 while( fscanf(pFile[0], "%lf %lf", &fff, &aaa) > 0 && iflag < 999 )
//			 {

			while( fgets(string,numBytes,pFile[0])> 0  && iflag < 999)
            {

				if(strrchr(string,',' ) != NULL )
				{
					sscanf(string,"%lf, %lf", &fff,&aaa );
				}
				else
				{
					sscanf(string,"%lf %lf", &fff,&aaa );
				}

				if( fff > 1.0e-04)
				{
					f[i]=fff;

					a[i]=aaa;

					i++;
				}
				if(i>=MAX)break;  

			 }
			 fclose( pFile[0]);

          }
		  if(imethod==3)
		  {
			  int ilibrary=0;

			  int Lflag=0;

			  while(Lflag==0)
			  {

				printf("\n Select PSD:        \n"); 

				printf("\n 1=MIL-STD-1540B QTP  ");
				printf("\n 2=MIL-STD-1540B ATP  ");
				printf("\n 3=NAVMAT P9495		");
				printf("\n 4=NASA CEQATR AVT  \n");

				scanf("%d", &ilibrary);


				if(ilibrary==1 || ilibrary==2 )
				{

					f[0]=20.;	
					a[0]=0.0053;

					f[1]=150.;	
					a[1]=0.04;

					f[2]=600.;
					a[2]=0.04;

					f[3]=2000.; 
					a[3]=0.0036;

					n=4;

					Lflag=1;

				}

				if(ilibrary==1)
				{
					a[0]*=4.;
					a[1]*=4.;
					a[2]*=4.;
					a[3]*=4.;

				}

				if(ilibrary==3)
				{

					f[0]=20.;	
					a[0]=0.01;

					f[1]=80.;	
					a[1]=0.04;

					f[2]=350.;
					a[2]=0.04;

					f[3]=2000.; 
					a[3]=0.007;

					n=4;

					Lflag=1;

				}

				if(ilibrary==4)
				{

					f[0]=20.;	
					a[0]=0.01;

					f[1]=80.;	
					a[1]=0.04;

					f[2]=800.;
					a[2]=0.04;

					f[3]=2000.; 
					a[3]=0.01;

					n=4;

					Lflag=1;

				}

			  }
			  i=n;

		  }

		  printf("\n");


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

		  scanf("%lf",&duration);		  


		  if(iflag < 999)
		  {

		     fileoper();

              m=i;


             for(i=0; i<m; i++)
			 {
				if((i >0) && (f[i] < f[i-1]))
				{

					printf(" Error:  frequencies must ascend. \n");

					exit(1);

				}
				if(f[i] < 0.)
				{

					printf(" all input numbers must be > 0. \n");

					printf(" %g  %g \n",f[i],a[i]);

					exit(1);

				} 
				if(a[i] < 0.)
				{

					printf(" all input numbers must be > 0. \n");

					printf(" %g  %g \n",f[i],a[i]);

					exit(1);

				}

				fprintf( pFile[1],"%lf \t %lf \n",f[i],a[i]);

				fprintf( pFile[6],"%12.4g \t %12.4g \n",f[i],a[i]);

			 } 

		     fclose(pFile[1]);
			 fclose(pFile[6]);

	      
			 printf("\n reference 1 \n");

			 ii=1;

		     calc(a);
			 

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

			 interp();

		     oa=grms;

	
			 printf("\n*********Input Parameters************************ \n");

             printf("     Freq        PSD        Slope       Slope  \n");
             printf("     (Hz)      (G^2/Hz)     log-log    (dB/octave) \n");
		     printf("   %6.1lf    %lf \n", f[0],a[0]);    



			 fprintf(pFile[5],"\n*********Input Parameters************************ \n");

             fprintf(pFile[5],"     Freq        PSD        Slope       Slope  \n");
             fprintf(pFile[5],"     (Hz)      (G^2/Hz)     log-log    (dB/octave) \n");
		     fprintf(pFile[5],"   %6.1lf    %lf \n", f[0],a[0]);    



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

                db=( 10.*log10(2.) )*s[i];

                printf("   %6.1lf    %lf    %lf     %lf  \n",f[i+1],a[i+1],s[i],db);

			   fprintf(pFile[5],"   %6.1lf    %lf    %lf     %lf  \n",f[i+1],a[i+1],s[i],db);

			   			 

			   flast=f[i+1];

			 }

			 fstart=f[0];


             for(i=0; i<m; i++)
			 {
                 if( f[i] >= 1.)  // highpass filter
				 {
					v[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),2.));

				 }
				 else
				 {
					 v[i]=1.0e-60;

				 }
			 }

	
			 ii=2;

	         calc(v);

             ov=grms;


	         for(i=0; i<m; i++)
			 {
                if( f[i] >= 1.)  // highpass filter
				{
                   d[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.));

				}
				
				if(f[i]<1.)
				{
					d[i]=1.0e-60;

				}


				if(d[i]<1.0e-80)
				{

						d[i]=1.0e-80;

/*
                        printf(" %8.4g \n",(pow( 386.,2.) ));  
                        printf(" %8.4g \n",a[i]/(pow( ( tpi*f[i] ),4.)));  
                        printf(" %8.4g \n",(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.)));  

*/
						printf("\n warning: f[%ld]=%8.4g  a[%ld]=%8.4g  d[%ld]=%8.4g \n",i,f[i],i,a[i],i,d[i]);  

					    printf("\n press any key to continue \n");

						getch();

				}

			 }

	         for(i=0; i<m; i++)
			 {
				 if(d[i]<1.0e-80)
				 {

					 printf("pre-calc  f[%ld]=%8.4g  a[%ld]=%8.4g  d[%ld]=%8.4g \n",i,f[i],i,a[i],i,d[i]);

					 exit(1);

				 }
			 }

		
			 ii=3;

	         calc(d);

             od=grms;


             fprintf(pFile[5],"\n Overall acceleration = %10.3lf GRMS\n", oa);



             cout << "\n Overall Acceleration = " << oa   << " GRMS "              << endl;

             cout << " Overall Velocity     = " << ov   << " (inch RMS)/sec"     << endl;

             cout << " Overall Displacement = " << od   << " (inch RMS)"         << endl;

             cout << " Overall Displacement = " << od*3 << " (inch 3-sigma)" << endl;                     


             cout << "\n  Input points written to file:  vrs.in "<< endl;


	         vrs();



			 cout << "\n\n Acceleration Output files"   << endl;

			 cout << "      vrs.grp - fn(Hz) and Response(GRMS)    " << endl;

			 cout << "    tsvrs.grp - fn(Hz) and Response(3-sigma) " << endl;

			 cout << "   ns_vrs.grp - fn(Hz) and Response(n-sigma) " << endl;

			 cout << "                where n =sqrt(2*ln(fnT))     " << endl;

			 cout << "   " << endl;

			 cout << "      vrs.out - text output"                   << endl;




			 cout << "\n\n Relative Displacement Output files"        << endl;

			 cout << "      rd_vrs.grp - fn(Hz) and Response(inch RMS)  " << endl;

			 cout << "   ts_rd_vrs.grp - fn(Hz) and Response(inch 3-sigma)  " << endl;

			 cout << "   ns_rd_vrs.grp - fn(Hz) and Response(inch n-sigma)  " << endl;



			 cout << "\n\n n-sigma Output file:  ns.grp"        << endl;


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

            

		//	 cout << "  Another calculation?  1=yes 2=no        "<< endl;

        //    cin >> ic;

			 ic=9999;

 	     //    cout << endl;


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

			 getch();


		  }
	  }
}
void interp(void)
{

	df=1.;

	fi[0]=f[0];

	ai[0]=a[0];



    for(i=1; i < MAX; i++) 
	{

		fi[i]=f[0]+df*i;



//		printf(" i=%d  fi[i]=%lf m=%d f[m]=%lf \n",i,fi[i],m,f[m]);


        if( fi[i] > f[m-1] ){break;}


        for(j=0; j < m-1; j++)
		{
			if( ( fi[i] >= f[j] ) &&  ( fi[i] <= f[j+1] )  )
			{		

				ai[i]=a[j]*pow( ( fi[i] / f[j] ), s[j] );

				break;

			}
		}
	}

	last=i-1;



//	printf("\n last= %d \n", last );

	for( i=0; i<last; i++)
	{
		fprintf( pFile[2], "%lf \t %lf\n", fi[i], ai[i] );

	}
    fclose( pFile[2] );


}
void vrs(void)
{
	double oct=1./24.;
	double c1, c2;
    double damp;
    double fn;
	double inch_rms;
	double rho;
	double sum;
	double t, tdr;
	double upper=10000.;

	int k;


	printf("\n*********Output Parameters************************ \n");

	fprintf(pFile[5],"\n*********Output Parameters***************\n");


	long idamp;


	printf("\n Enter damping format. \n 1= damping ratio   2= Q \n");

	scanf("%ld",&idamp);


	if(idamp == 1 )
	{
	  printf("\n Enter damping ratio (typically 0.05) ");

      scanf("%lf",&damp);

	}
	else
    {

	  printf("\n Enter amplification factor Q (typically 10) ");

      scanf("%lf",&damp);

	  damp = 1./(2.*damp);

	}


	fprintf(pFile[5],"\n damping = %5.3lf \n\n",damp);

	fprintf(pFile[5],"\n     fn       Response     Response");
	fprintf(pFile[5],"\n    (Hz)      (GRMS)       (3 sigma) \n");



	fn=5.;


//    printf("\n last= %d \n",last);



    printf("\n calculating..");

    k=0;


	for(i=0; i < 1000; i++)
	{
		k++;

		if(k==15)
		{
		   printf(".");

		   k=0;
        }


//   absolute acceleration

        sum=0.; 

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

		    rho = fi[j]/fn;

			tdr=2.*damp*rho;
            
			c1=0.;
			c2=0.;

            c1= pow( tdr, 2. );
			c2= pow( ( 1.- pow(rho,2.)), 2.);

			t= (1.+ c1 ) / ( c2 + c1 );

            sum+=t*ai[j]*df;

		}
        grms=sqrt(sum);

		if(grms>f_peak)
		{
			f_peak=grms;
			f_freq=fn;
		}

	    fprintf(pFile[3],"%lf \t %lf\n",fn,grms);
	    fprintf(pFile[4],"%lf \t %lf\n",fn,grms*3.);	


		if(fn>=fstart && fn<=flast)
		{
			double nnf=sqrt(2*log(fn*duration));

			fprintf(pFile[8],"%lf \t %lf\n",fn,grms*nnf);

			fprintf(pFile[10],"%lf \t %lf\n",fn,nnf);
		}


		fprintf(pFile[5],"%10.2lf \t %10.3lf \t %10.3lf\n",fn, grms, grms*3.);

		
//   relative displacement

		sum=0.; 

		for(j=0; j <last; j++) 
        {
			c1=0.;
			c2=0.;

            c1= pow( ( pow(fn,2.)-pow(fi[j],2.) ),2. );

			c2= pow( ( 2.*damp*fn*fi[j]),2.);

			t= 1. / ( c2 + c1 );

			if(j==0){t/=2.;}

            sum+=t*ai[j]*df;

		}

        inch_rms=sqrt(sum)*386./pow(tpi,2.);


	    fprintf(pFile[7],"%lf \t %12.3e\n",fn,inch_rms);
	    fprintf(pFile[11],"%lf \t %12.3e\n",fn,3.*inch_rms);


		if(fn>=fstart && fn<=flast)
		{
			fprintf(pFile[9],"%lf \t %12.3e\n",fn,inch_rms*sqrt(2*log(fn*duration)));

        }

//  frequency advance		

		fn*=pow(2.,oct);

		if(fn > upper){break;}
	}

	fclose( pFile[3] );
	fclose( pFile[4] );
	fclose( pFile[5] );
	fclose( pFile[7] );
	fclose( pFile[8] );
	fclose( pFile[9] );
	fclose( pFile[10] );
	fclose( pFile[11] );

	printf("\n\n peak VRS at %8.4g Hz  %8.4g GRMS \n",f_freq,f_peak);

}
void calc(double *level)
{
		  ra=0.;
          grms=0.;

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

		  }

		  for(i=0; i < m; i++)
		  {
			  if(level[i]<1.0e-90)
			  {
				    printf("* f[%ld]=%8.4g  level[%ld]=%10.4e \n",i,f[i],i,level[i]);

					printf("\n error 1:  %s \n",amp_label[ii]);

					exit(1);

			  }
		  }


//  calculate slopes

          nn=m-1;


		  double lower_f=0.0001;

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

			 if(f[i] < lower_f){f[i]=lower_f;}


			  if(f[i]>f[i+1])
			  {
					printf("\n error 2:  %s\n",amp_label[ii]);

					exit(1);

			  }
			  if(f[i]< lower_f)
			  {
					printf("\n error 3:  %s \n",amp_label[ii]);

					exit(1);

			  }
			  if(f[i+1]< lower_f)
			  {
					printf("\n error 4:  %s \n",amp_label[ii]);

					exit(1);

			  }
			  if(level[i]<1.0e-90)
			  {
				    printf("***** f[%ld]=%8.4g  level[%ld]=%10.4e \n",i,f[i],i,level[i]);

					printf("\n error 5:  %s \n",amp_label[ii]);

					exit(1);

			  }


              s[i]=log10( level[i+1]/ level[i] )/log10( f[i+1]/f[i] );

          }


	      for( i=0; i < nn; i++)
		  {
              if( fabs( s[i]+1.)>1.0e-03 )
			  {	   

                 ra+= ( level[i+1] * f[i+1]- level[i]*f[i])/( s[i]+1.);

			  }

		      else
			  {

                 ra+= level[i]*f[i]*log( f[i+1]/f[i]);

			  }

		  }

          grms=sqrt(ra);

}
void fileoper(void)
{
       strcpy(filename[1], "a.psd");
       pFile[1]=fopen(filename[1], "w");

	   strcpy(filename[2], "a.int");
       pFile[2]=fopen(filename[2], "w");

	   strcpy(filename[3], "vrs.grp");
       pFile[3]=fopen(filename[3], "w");

	   strcpy(filename[4], "tsvrs.grp");
       pFile[4]=fopen(filename[4], "w");

	   strcpy(filename[5], "vrs.out");
       pFile[5]=fopen(filename[5], "w");

	   strcpy(filename[6], "vrs.in");
       pFile[6]=fopen(filename[6], "w");

	   strcpy(filename[7], "rd_vrs.grp");
       pFile[7]=fopen(filename[7], "w");

	   strcpy(filename[8], "ns_vrs.grp");
       pFile[8]=fopen(filename[8], "w");

	   strcpy(filename[9], "ns_rd_vrs.grp");
       pFile[9]=fopen(filename[9], "w");


	   strcpy(filename[10], "ns.grp");

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


	   strcpy(filename[11], "ts_rd_vrs.grp");
       pFile[11]=fopen(filename[11], "w");

}
void header(void)
{
	long numBytes=300;

	char string[300];

	int nmax=0;
	long j=1;

	while(fgets(string,numBytes,pFile[0] )>0)
	{
//		printf("%ld.  %s\n",j,string);

		if( strrchr(string,'A' ) != NULL ){nmax=j;}
		if( strrchr(string,'B' ) != NULL ){nmax=j;}
		if( strrchr(string,'C' ) != NULL ){nmax=j;}
		if( strrchr(string,'D' ) != NULL ){nmax=j;}
		if( strrchr(string,'E' ) != NULL && strstr(string,"E+")==NULL && strstr(string,"E-")==NULL ){nmax=j;}
		if( strrchr(string,'F' ) != NULL ){nmax=j;}
		if( strrchr(string,'G' ) != NULL ){nmax=j;}
		if( strrchr(string,'H' ) != NULL ){nmax=j;}
		if( strrchr(string,'I' ) != NULL ){nmax=j;}
		if( strrchr(string,'J' ) != NULL ){nmax=j;}
		if( strrchr(string,'K' ) != NULL ){nmax=j;}
		if( strrchr(string,'L' ) != NULL ){nmax=j;}
		if( strrchr(string,'M' ) != NULL ){nmax=j;}
		if( strrchr(string,'N' ) != NULL ){nmax=j;}
		if( strrchr(string,'O' ) != NULL ){nmax=j;}
		if( strrchr(string,'P' ) != NULL ){nmax=j;}
		if( strrchr(string,'Q' ) != NULL ){nmax=j;}
		if( strrchr(string,'R' ) != NULL ){nmax=j;}
		if( strrchr(string,'S' ) != NULL ){nmax=j;}
		if( strrchr(string,'T' ) != NULL ){nmax=j;}
		if( strrchr(string,'U' ) != NULL ){nmax=j;}
		if( strrchr(string,'V' ) != NULL ){nmax=j;}
		if( strrchr(string,'W' ) != NULL ){nmax=j;}
		if( strrchr(string,'X' ) != NULL ){nmax=j;}
		if( strrchr(string,'Y' ) != NULL ){nmax=j;}
		if( strrchr(string,'Z' ) != NULL ){nmax=j;}
//
		if( strrchr(string,'a' ) != NULL ){nmax=j;}
		if( strrchr(string,'b' ) != NULL ){nmax=j;}
		if( strrchr(string,'c' ) != NULL ){nmax=j;}
		if( strrchr(string,'d' ) != NULL ){nmax=j;}
		if( strrchr(string,'e' ) != NULL  && strstr(string,"e+")==NULL && strstr(string,"e-")==NULL ){nmax=j;}
		if( strrchr(string,'f' ) != NULL ){nmax=j;}
		if( strrchr(string,'g' ) != NULL ){nmax=j;}
		if( strrchr(string,'h' ) != NULL ){nmax=j;}
		if( strrchr(string,'i' ) != NULL ){nmax=j;}
		if( strrchr(string,'j' ) != NULL ){nmax=j;}
		if( strrchr(string,'k' ) != NULL ){nmax=j;}
		if( strrchr(string,'l' ) != NULL ){nmax=j;}
		if( strrchr(string,'m' ) != NULL ){nmax=j;}
		if( strrchr(string,'n' ) != NULL ){nmax=j;}
		if( strrchr(string,'o' ) != NULL ){nmax=j;}
		if( strrchr(string,'p' ) != NULL ){nmax=j;}
		if( strrchr(string,'q' ) != NULL ){nmax=j;}
		if( strrchr(string,'r' ) != NULL ){nmax=j;}
		if( strrchr(string,'s' ) != NULL ){nmax=j;}
		if( strrchr(string,'t' ) != NULL ){nmax=j;}
		if( strrchr(string,'u' ) != NULL ){nmax=j;}
		if( strrchr(string,'v' ) != NULL ){nmax=j;}
		if( strrchr(string,'w' ) != NULL ){nmax=j;}
		if( strrchr(string,'x' ) != NULL ){nmax=j;}
		if( strrchr(string,'y' ) != NULL ){nmax=j;}
		if( strrchr(string,'z' ) != NULL ){nmax=j;}
//
		if( strrchr(string,'!' ) != NULL ){nmax=j;}
		if( strrchr(string,'@' ) != NULL ){nmax=j;}
		if( strrchr(string,'#' ) != NULL ){nmax=j;}
		if( strrchr(string,'$' ) != NULL ){nmax=j;}
		if( strrchr(string,'%' ) != NULL ){nmax=j;}
		if( strrchr(string,'^' ) != NULL ){nmax=j;}
		if( strrchr(string,'&' ) != NULL ){nmax=j;}
		if( strrchr(string,'*' ) != NULL ){nmax=j;}
		if( strrchr(string,'/' ) != NULL ){nmax=j;}
		if( strrchr(string,'>' ) != NULL ){nmax=j;}
		if( strrchr(string,'<' ) != NULL ){nmax=j;}
		if( strrchr(string,'?' ) != NULL ){nmax=j;}
		if( strrchr(string,'=' ) != NULL ){nmax=j;}
		if( strrchr(string,';' ) != NULL ){nmax=j;}

		j++;
	}
	rewind(pFile[0]);

	printf("\n %ld header lines detected \n\n",nmax);

	for(j=0;j<nmax;j++)
	{
		fgets(string,numBytes,pFile[0] );
		printf(" %s ",string);
	}
}