/*

     Vibration response spectrum   

*/



#include <math.h>

#include <iostream.h>

#include <string.h>

#include <stdlib.h>

#include <stdio.h>

#include <conio.h>



#define MAX 5000



void fileoper(void);

void calc(double *level);

void vrs(void);

void interp(void);

void fatigue_file(void);





double level[MAX];

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

double df, db,ra,grms;

double oa, ov, od;



double fe,dfact;



double Q;



long fe_index;



int iflag;

int i, ic, imethod;

int j;

int last;

int m, n, nn;



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

double tpi=2.*pi;



double duration;



double damp;





FILE *pFile[10];

char filename[10][FILENAME_MAX];



void main()

{

	   cout.precision(3);



	   iflag=0;





       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.;



       printf("											\n");

       printf(" vrs_fatigue.cpp  ver 1.1, June 20, 2008	\n");

       printf("											\n");

       printf(" by Tom Irvine							\n");

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

	   printf("											\n");

       



       printf(" This program calculates the acceleration vibration response \n"); 

       printf(" spectrum of a single-degree-of-freedom system to a          \n");

       printf(" power spectral density input to the base of the system.   \n\n");



       printf(" The power spectral density may vary arbitrarily with frequency. \n");



       printf(" Natural frequency is an independent variable.  Damping is fixed.\n\n");



       

	   printf(" This program also implements the fatigue calculation     \n");

	   printf(" per Kent Hardy & Mike Marchak, ME File 030-267,         \n");             

	   printf(" Use of Vibration Response Spectrum (VRS) and Fatigue     \n");

	   printf(" Equivalence Methods for Qualification Status Comparisons, 2004 \n\n");	   

	   





	   ic=1;

	   while(ic==1)

	   {

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

          cout << " 1=keyboard  2=file       \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]);

			  }

          }

          else

		  {



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



			 if(pFile[0] == NULL)

			 {

				 printf("Input file error \n");

				 iflag=999;

			 }



			 i=0;



             double fff,aaa;



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

			 {

				   if( fff > 1.0e-04)

				   {

						f[i]=fff;

						a[i]=aaa;



						i++;

                   }



		

                  if(i>=MAX)break;  

			 }

			 fclose( pFile[0]);



          }

		  cout << " " << endl;

		  

		  

		  if(iflag < 999)

		  {

		     fileoper();



              m=i;



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

			 {

             if((i >0) && (f[i] < f[i-1]))

			 {

                cout << " Error:  frequencies must ascend. \n";

             }

             if(f[i] < 0.)

			 {

                cout << " all input numbers must be > 0.\n";

                cout << f[i] << a[i];

			 } 

             if(a[i] < 0.)

			 {

                cout << " all input numbers must be > 0.\n";

                cout << f[i] << a[i];

			 }

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

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



			 } 

		     fclose(pFile[1]);

			 fclose(pFile[6]);





			cout << " Enter duration (sec) \n";

			cin  >> duration;





			cout << " Select fatigue exponent \n";

			cout << " 1=4.0    2=6.4       \n";

			cin  >> fe_index;



			if(fe_index==1)

			{

				dfact = 8.;

				fe=4.0;

			}

			else

			{

				dfact = 65.93;

				fe=6.4;

			}





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

	}



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





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



    

	fatigue_file();

	



		     calc(a);

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

			 }

          

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

			 {

                 if( f[0] >= 1.)

				 {

                 v[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),2.));

				 }       

			 }

	         calc(v);

             ov=grms;



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

			 {

                if( f[0] >= 1.)

				{

                   d[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.));

				}       

			 }

	         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 three sigma)" << endl;                     



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



	         vrs();



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

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

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

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





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

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





			 cout << " \n Please call the *.txt files into your own plotting program.\n";

            

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

        //    cin >> ic;

			 ic=9999;

 	     //    cout << endl;







			 printf("\n Fatigue response spectrum output file:\n");

			 printf("\n        %s - fn(Hz) and Fatigue (R) \n\n",filename[8]);







			 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 %10.5g\n", fi[i], ai[i] );

	}

    fclose( pFile[2] );



}



void vrs(void)

{



	double oct=1./24.;

	double c1, c2;

    

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

	



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

        

	

	    fprintf(pFile[3],"%lf %10.5g\n",fn,grms);

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

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



		double sigma2=pow(grms,2.);



		double two_sigma2=2.*sigma2;





		double rfd = (fn*duration/sigma2);

		       rfd*=dfact;



		double exponent=(fe+1.)/2.;



        double term=pow((2./two_sigma2),exponent);


		double denom=2.*term;

			   rfd/=denom;



		       rfd*=sqrt(pi*two_sigma2);



	    fprintf(pFile[8],"%lf \t %10.5g\n",fn,rfd);   // fatigue file



		



//   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 %12.3e\n",fn,inch_rms);





//  frequency advance		



		fn*=pow(2.,oct);

		if(fn > upper){break;}

	}

	fclose( pFile[3] );

	fclose( pFile[4] );

	fclose( pFile[5] );

	fclose( pFile[7] );

}

	   

void calc(double *level)

{

		  ra=0.;

          grms=0.;



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

		  {

			  s[i]=0.;

		  }

		  



//  calculate slopes



          nn=m-1;



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

		  {

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

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

          }



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

		  {

              if(s[i] < -1.0001 ||  s[i] > -0.9999 )

			  {	   

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

			  }

		      else

			  {

                 ra+= a[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.txt");

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



	   strcpy(filename[4], "tsvrs.txt");

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

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



}

void fatigue_file(void)

{

	//



  	   char temp1[10],temp2[10];





	   printf("\n  Q=%lf \n",Q);

	   printf("\n fe=%lf \n",fe);



	   sprintf (temp1,"%ld", long(Q) );

	   

	   sprintf (temp2,"%ld", long(fe) );



	   if(fe_index==1)

	   {

	   }

	   else

	   {

			strcat(temp2,"p4" );

	   }





	   strcpy(filename[8],"random_Q");

	   strcat(filename[8],temp1);

	   strcat(filename[8],"_b");

	   strcat(filename[8],temp2);



	   strcat(filename[8],".txt");



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



//

}