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

const long MAX=2000000;
const long NUM=500;

void files(void);
void maxcheck(void);
void finalcheck(void);
void initialize_filter_coefficients(const double &sr,const double &dt);

void initialize_arrays(void);
void enter_parameters(void);
void delta_time(void);
void response(void);



double a1[NUM],a2[NUM];
double b1[NUM],b2[NUM];
double c2[NUM];

double alpha;

double cosd;
double damp;
double domegadt;
double dt;
double fn[NUM];

long idamp;
long last;

double omega, omegad;

double sind;
double sr;
double t, t1, tb;

double xb[NUM], xbb[NUM];
double x[NUM],xmin[NUM],xmax[NUM],xabs[NUM];

double rdmax[NUM],rdmin[NUM],rdabs[NUM];
double rd[NUM],rdb[NUM],rdbb[NUM];

double y;
double yb[NUM];
double yin[MAX];


const double pi=atan2(0.,-1.);
const double tpi=2.*pi;


long i,j,jlast;

long jnum=0;
long ntotal;
long iunit,ix;

int ire;

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


int mflag=0;

int main( int argc, char *argv[] )	
{

	if(argc > 1 )
	{

		printf("\n\n  argc = %d ", argc );

		printf("\n argv[0] %s ", argv[0] );

		printf("\n argv[1] %s ", argv[1] );

		strcpy (filename[0], argv[1] );



		mflag=1;

	}

	printf("\n quakeSRS.cpp  version 2.5,  May 8, 2010 \n");
	printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

	printf("\n This program calculates the shock response spectrum ");
	printf("\n to base excitation. \n");

	initialize_arrays();

	files();

	delta_time();
  
	enter_parameters();
   
	initialize_filter_coefficients(sr,dt);
 
	response();


	fclose(pFile[1]);
	fclose(pFile[2]);
	fclose(pFile[3]);
	fclose(pFile[4]);



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

	exit(1);


}

void delta_time(void)
{

/********************* determine delta t *****************************/

	i=0;

	last=0;

	ntotal=0;

	double dtmin=1.0e+90;

	double dtmax=-dtmin;





   while( fscanf(pFile[0], "%lf %lf", &t, &yin[i]) > 0 )
   {

      if(i==0)t1=t;


/*      if(i>=2000)break;    */

      if(i>0)
      {
		  if(t<tb)
		  {
			 printf("\n Time error.  Time values must ascend.");

			 return;
		  }

          if(dtmin>(t-tb))
		  {
			  dtmin=(t-tb);
		  }

          if(dtmax<(t-tb))
		  {
			  dtmax=(t-tb);

		  }
     }

      tb=t;

      i++;

   }

   fclose(pFile[0]);



   printf("\n dtmax = %8.4g   dtmin = %8.4g \n",dtmax,dtmin);


   if( dtmin < 0.99 *dtmax)
   {
	   printf("\n Time Step Error:   ");
	   printf("\n The time step must be constant. \n");
	   printf("\n Press any key to exit \n");

	   getch();

	   exit(1);

   }

   if(i > 1)
   {
      dt=(t-t1)/(i-1);

      sr=1./dt;

   }	 

  fclose(pFile[0]);

   last=i;

   ntotal=last;

  printf("\n\n %ld samples read \n",ntotal);

   sr=1./dt;


   printf("\n  dt=%g sec   sr=%g samples/sec \n", dt, sr);


}

void response(void)
{

   long iend;

   for(j=0; j<jlast; j++)  //jlast
   {

	   i=0;

       iend = ntotal + long(1./(dt*fn[j]));  // primary+residual
       
        if(ire==1)
        {
              iend = ntotal + long(2./(dt*fn[j]));  // primary+2*residual
        }
        else
        {
              iend = ntotal;  // primary      
        }      


		while( i < iend  ) // iend
		{

			 if( i< last)
			 {
				 y= yin[i];
			 }
			 else
			 {
				 y=0.;
			 }


			rd[j]=   a1[j]* rdb[j]
					+a2[j]*rdbb[j]
					+c2[j]*yb[j];


			x[j]=    a1[j]* xb[j]
					+a2[j]*xbb[j]
					+b1[j]*y
					+b2[j]*yb[j];

 
			maxcheck();


			rdbb[j]=rdb[j];
			 rdb[j]= rd[j];
			 
			 xbb[j]=xb[j];
			  xb[j]= x[j];


			 if( i< last)
			 {
				 yb[j]= yin[i];

			 }
			 else
			 {
				 yb[j]=0.;
			 }

 		     i++;

		}

   }


   finalcheck();


		printf("\n\n Output Files: \n");
		
		printf("\n  %s       -  relative displacement",filename[1]);
		printf("\n  %s   -  pseudo velocity \n",filename[2]);
		printf("\n  %s   -  pseudo acceleration ",filename[3]);
		printf("\n  %s      -  absolute acceleration ",filename[4]);
        printf("\n  %s       -  positive & negative acceleration \n",filename[5]);

	if( ix==1)
	{	

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

			if( iunit==1)
			{
                
				if( j==0){ printf("\n Units: inch, inch/sec, G \n"); }


				fprintf(pFile[1],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha ); 
 				fprintf(pFile[2],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha*tpi*fn[j] );   
  				fprintf(pFile[3],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha*pow( (tpi*fn[j]),2. )/386.); 
  				fprintf(pFile[4],"%12.6e \t %12.6e\n",fn[j],xabs[j]*alpha/386.); 
  				fprintf(pFile[5],"%12.6e \t %12.6e \t %12.6e\n",fn[j],xmax[j]*alpha/386.,fabs(xmin[j])*alpha/386.); 
			}

			else
			{

				if( j==0 ){ printf("\n Units: mm, mm/sec, m/sec^2 \n"); }

				fprintf(pFile[1],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha ); 
 				fprintf(pFile[2],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha*tpi*fn[j] );   
  				fprintf(pFile[3],"%12.6e \t %12.6e\n",fn[j],rdabs[j]*alpha*pow( (tpi*fn[j]),2. )/1000.);
  				fprintf(pFile[4],"%12.6e \t %12.6e\n",fn[j],xabs[j]*alpha/1000.);
			    fprintf(pFile[5],"%12.6e \t %12.6e \t %12.6e\n",fn[j],xmax[j]*alpha/1000.,fabs(xmin[j])*alpha/1000.); 
			}

		}

	}

	if(ix != 1)
	{

		for(j=jlast-1; j>=0; j--)
		{
			if( iunit==1)
			{

				if( j==0){ printf("\n Units: inch, inch/sec, G \n"); }


				fprintf(pFile[1],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha ); 
 				fprintf(pFile[2],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha*tpi*fn[j] );   
  				fprintf(pFile[3],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha*pow( (tpi*fn[j]),2. )/386.); 
  				fprintf(pFile[4],"%12.6e \t %12.6e\n",1./fn[j],xabs[j]*alpha/386.); 		
  				fprintf(pFile[5],"%12.6e \t %12.6e \t %12.6e\n",1./fn[j],xmax[j]*alpha/386.,fabs(xmin[j])*alpha/386.); 
			}
			else
			{

				if( j==0 ){ printf("\n Units: mm, mm/sec, m/sec^2 \n"); }

				fprintf(pFile[1],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha ); 
 				fprintf(pFile[2],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha*tpi*fn[j] );   
  				fprintf(pFile[3],"%12.6e \t %12.6e\n",1./fn[j],rdabs[j]*alpha*pow( (tpi*fn[j]),2. )/1000.);
  				fprintf(pFile[4],"%12.6e \t %12.6e\n",1./fn[j],xabs[j]*alpha/1000.);	
  	            fprintf(pFile[5],"%12.6e \t %12.6e \t %12.6e\n",1./fn[j],xmax[j]*alpha/1000.,fabs(xmin[j])*alpha/1000.);                 		

			}

		}

    }

   printf("\n\n Please call these file into your own plotting program.\n\n");

}

void files(void)
{
	/******************************** open files *************************/

   	if(mflag==0)
	{

		printf( "\n The base input file must contain two columns: \n\n");
		printf( " time(sec)  accel(G)    \n\n");

		printf( " \n An error message will occur if the file does not exist.\n");
		printf( " Input filename \n\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]);


	strcpy(filename[1], "rel_disp.srs");
	pFile[1]=fopen(filename[1], "w");

	strcpy(filename[2], "pseudo_Velox.srs");
	pFile[2]=fopen(filename[2], "w");

	strcpy(filename[3], "pseudo_Accel.srs");
	pFile[3]=fopen(filename[3], "w");

	strcpy(filename[4], "abs_Accel.srs");
	pFile[4]=fopen(filename[4], "w");
	
	strcpy(filename[5], "pn_Accel.srs");
	pFile[5]=fopen(filename[5], "w");

}
void maxcheck(void)
{
			if( rd[j] > rdmax[j] ) 
			{
				rdmax[j]=rd[j];
			}
			if(rd[j] < rdmin[j])
			{
				rdmin[j]=rd[j];
			}

			if( x[j] > xmax[j] ) 
			{
				xmax[j]=x[j];
			}
			if(x[j] < xmin[j])
			{
				xmin[j]=x[j];
			}

}
void finalcheck(void)
{
   for(j=0; j<jlast; j++)
   {
		rdabs[j]=rdmax[j];
		if(rdabs[j]< fabs(rdmin[j]))
		{ 
			rdabs[j]=fabs(rdmin[j]);
		}

		xabs[j]=xmax[j];
		if(xabs[j]< fabs(xmin[j]))
		{ 
			xabs[j]=fabs(xmin[j]);
		}
   }

}

void initialize_filter_coefficients(const double &sr,const double &dt)
{

	 fn[0]=0.01;

   if( ix==1)
   {

		printf("\n\n Enter starting frequency (Hz) \n");

		scanf("%lf",&fn[0]);
		printf("\n");

   }  


   for(j=1; j<NUM; j++)
   {
	   fn[j]=fn[j-1]*pow(2.,(1./12.));

	   if(fn[j] > sr/4. ){break;}
   }

   jlast=j-1;


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

      omega=2.*pi*fn[j];

      omegad=omega*sqrt(1.-pow(damp,2));

      cosd=cos(omegad*dt);

      sind=sin(omegad*dt);

      domegadt=damp*omega*dt;

      a1[j]=2.*exp(-domegadt)*cosd;
      a2[j]=-exp(-2.*domegadt);

      b1[j]=2.*domegadt;
      b2[j]=omega*dt*exp(-domegadt);
      b2[j]*=( (omega/omegad)*(1.-2.*pow(damp,2))*sind -2.*damp*cosd );

	  c2[j]=-(dt/omegad)*exp(-domegadt)*sind;

   }

}

void initialize_arrays(void)
{

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


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

	  rdabs[i]=0.;
      rdmax[i]=0.;
      rdmin[i]=0.;

	  xabs[i]=0.;
      xmax[i]=0.;
      xmin[i]=0.;
   
	    rd[i]=0.;
	   rdb[i]=0.;
      rdbb[i]=0.;

		 x[i]=0.;
		xb[i]=0.;
       xbb[i]=0.;
	    yb[i]=0.;

	}

}

void enter_parameters(void)
{	

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

	}

	printf("\n Enter output relative displacement unit. \n 1= inch   2= mm \n");

	scanf("%ld",&iunit);



	if(iunit==1)
	{
		alpha=386.;
	}
	else
	{
		alpha=9810.;
	}

	printf("\n Select X-axis Type:  1=natural frequency  2=period \n");
	scanf("%ld",&ix);

    printf("\n Include residual?   ");
    printf("\n   1=yes  2=no \n");
	scanf("%d",&ire);

}