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

void displace(void);
void velox(void);
void accel(void);

void initial(void);
void files(void);

void coefficients(void);
void enter_data(void);

void header(void);


#define CORE \
x[j]=0.;\
if(i<last)\
{\
yy=yin[i];\
}\
else\
{\
yy=0.;\
}\
x[j]=  a1[j]* xb[j]\
+a2[j]*xbb[j]\
+b1[j]*yy\
+b2[j]*yb[j];\
if(x[j] > xmax[j])\
{\
xmax[j]=x[j];\
tmax[j]=t;\
}\
if(x[j] < xmin[j])\
{\
xmin[j]=x[j];\
tmin[j]=t;\
}\
xbb[j]=xb[j];\
xb[j]= x[j];\
yb[j]= yy;\



#define MAX 1000000

double a1[1],a2[1];
double b1[1],b2[1];
double cosd;
double damp;
double domegandt,domegan;
double dt;
double fn[1];

long idamp;
long last;

double mass;
double omegan, omegad;

double sind;
double sr;
double t, t1, tb;
double tmax[1],tmin[1];
double xb[1], xbb[1];
double x[1],xmax[1],xmin[1];
double y,yy;
double yb[1];


double yin[MAX];

double velocity[MAX];
double displacement[MAX];

double stiffness;

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

double ccc;

double max_tf,min_tf;


long i,iend,j,jnum;

int iunits;


int ire;


FILE *pFile[6];
char filename[6][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 arbit_f.cpp  version 2.0, October 6, 2011 \n");
	printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

	printf("\n This program calculates the acceleration time history");
	printf("\n response of a single-degree-of-freedom system subjected");
	printf("\n to an input force.");

	enter_data();

    jnum=0;

    files();

    coefficients();

	initial();
	displace();

	initial();
	velox();

	initial();
	accel();


   	printf("\n\n\n Output Files: \n");



	if(iunits==2)
	{
		printf("\n       disp.out  - displacement response       (mm) vs time(sec)");
		printf("\n       velox.out - velocity     response   (mm/sec) vs time(sec)");
		printf("\n       accel.out - acceleration response  (m/sec^2) vs time(sec) \n");
		printf("\n trans_force.out - transmitted force  (N) vs time(sec) \n");

	}
	else
	{

		printf("\n       disp.out  - displacement response     (inch) vs time(sec)");

		printf("\n       velox.out - velocity     response (inch/sec) vs time(sec)");

		printf("\n       accel.out - acceleration response        (G) vs time(sec) \n");

		printf("\n trans_force.out - transmitted force          (lbf) vs time(sec) \n");

	}

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

	exit(1);
 }
 void displace(void)
 {

/*********************** initialize filter coefficients ***************/
	i=0;
	j=0;

	b1[j]=0.;
	b2[j]=(dt/(mass*omegad))*exp(-domegandt)*sind;

/********************** response calculation *******************/

	x[0]=0.;
	fprintf(pFile[1],"%12.6e %12.6e\n",t1,x[0] );

	displacement[0]=x[0];


   i=1;


   while( i < iend  )
   {

	    CORE


		if(iunits==2)
		{
			fprintf(pFile[1],"%12.6e %12.6e\n",dt*i+t1,x[j]*1000. );
		}
		else
		{
			fprintf(pFile[1],"%12.6e %12.6e\n",dt*i+t1,x[j] );
		}


	    displacement[i+1]=x[j];

		i++;
   }



   if(iunits==2)
   {
		printf("\n\n maximum disp = %12.3g mm",xmax[0]*1000.);
		printf("\n minimum disp = %12.3g mm",xmin[0]*1000.);

   }
   else
   {

		printf("\n\n maximum disp = %12.3g inch",xmax[0]);

		printf("\n minimum disp = %12.3g inch",xmin[0]);

   }


 }
 void velox(void)
 {

	double tforce;

/*********************** initialize filter coefficients ***************/
	i=0;
	j=0;

	b1[j]=dt/mass;
	b2[j]=dt*exp(-domegandt)/mass;
	b2[j]*=( (-domegan/omegad)*sind - cosd );

	max_tf=0.;
	min_tf=0.;

/********************** response calculation *******************/

	x[0]=0.;
	fprintf(pFile[2],"%12.6e %12.6e\n",t1,x[0] );   // ref 1

	velocity[0]=x[0];


	tforce=stiffness*displacement[0]+ccc*velocity[0];

	fprintf(pFile[4],"%12.6e %12.6e\n",t1,tforce);


	if(tforce>max_tf)
	{
		max_tf=tforce;
	}
	if(tforce<min_tf)
	{
		min_tf=tforce;
	}




	i=1;



	double ttt=t1;

	double tb=t1;

   while( i < iend  )
   {

	    CORE



        ttt=dt*i+t1;


		if(ttt<=tb)
		{

			printf("\n\n error  i=%ld  ttt=%8.4g   tb=%8.4g\n",i,ttt,tb);

			exit(1);

		}



		tb=ttt;


		if(iunits==2)
		{
			fprintf(pFile[2],"%12.6e %12.6e\n",ttt,x[j]*1000. );   // ref 2
		}

		else
		{

			fprintf(pFile[2],"%12.6e %12.6e\n",ttt,x[j] );   // ref 3

		}




	    velocity[i+1]=x[j];


		tforce=stiffness*displacement[i+1]+ccc*velocity[i+1];

		fprintf(pFile[4],"%12.6e %12.6e\n",ttt,tforce);

		if(tforce>max_tf)
		{
			max_tf=tforce;
		}
		if(tforce<min_tf)
		{
			min_tf=tforce;
		}

		i++;

   }



   if(iunits==2)
   {
		printf("\n\n maximum velocity = %12.3g mm/sec",xmax[0]*1000.);
		printf("\n minimum velocity = %12.3g mm/sec",xmin[0]*1000.);

   }
   else

   {

		printf("\n\n maximum velocity = %12.3g inch/sec",xmax[0]);

		printf("\n minimum velocity = %12.3g inch/sec",xmin[0]);

   }


 }
 void accel(void)
 {
/*********************** initialize filter coefficients ***************/
	i=0;
	j=0;

	double actual,amax,amin;

	amax=0.;
	amin=0.;

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


	b1[j]/=mass;
	b2[j]/=mass;

	b1[j]*=-1.;
	b2[j]*=-1.;


/********************** response calculation *******************/

	x[0]=0.;
//	fprintf(pFile[3],"%12.6e %12.6e\n",t1,x[0] );


	double ms=0.;

   while( i < iend  )
   {

	    CORE

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

		actual = x[j]+yy/mass;

		if(actual > amax){amax=actual;}
		if(actual < amin){amin=actual;}


		if(iunits==1)
		{

			actual/=386.;

		}
		fprintf(pFile[3],"%12.6e %12.6e\n",dt*i+t1,actual );



		ms+=pow(actual,2.);

		i++;
   }



   double rms = sqrt(ms/double(i));





   if(iunits==2)
   {
		printf("\n\n maximum accel = %12.3g m/sec^2",amax);
		printf("\n minimum accel = %12.3g m/sec^2",amin);

		printf("\n overall accel = %12.3g m/sec^2",rms);

		printf("\n\n maximum transmitted force = %12.3g N",max_tf);
		printf("\n minimum transmitted force = %12.3g N",min_tf);
   }

   else
   {

	    amax/=386.;
		amin/=386.;


		printf("\n\n maximum accel = %12.3g G",amax);
		printf("\n minimum accel = %12.3g G",amin);

		printf("\n overall accel = %12.3g G",rms);

		printf("\n\n maximum transmitted force = %12.3g lbf",max_tf);

		printf("\n minimum transmitted force = %12.3g lbf",min_tf);

   }

 }
 void initial(void)
 {
/******************************** initialize arrays ******************/
   for(i=0; i< 1; ++i)
   {
      xmax[i]=0.;
      xmin[i]=0.;
      tmax[i]=0.;
      tmin[i]=0.;
		 x[i]=0.;
		xb[i]=0.;
       xbb[i]=0.;
	    yb[i]=0.;
   }
 }
void files(void)
{
	/******************************** open files *************************/

	if(mflag==0)
	{



		if(iunits==2)
		{

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

			printf( " time(sec)  force(N)    \n");

		}
		else
		{

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

			printf( " time(sec)  force(lbf)    \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(" Failed to open file: %s \n", filename[0]);
      exit(1);
	}
	else
	{
      printf(" File: %s opened. \n", filename[0]);
	}

	strcpy(filename[1], "disp.out");
	strcpy(filename[2], "velox.out");
	strcpy(filename[3], "accel.out");
	strcpy(filename[4], "trans_force.out");


	for(i=1; i<=4; i++)
	{
		pFile[i]=fopen(filename[i], "w");
	}
}
void coefficients(void)
{
	/********************* determine delta t *****************************/


   header();


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

      if(i>0)
      {
	    if(t<tb)
	    {
	       printf("\n Time error.  Time values must ascend.");
	       exit(1);
	    }
      }
      tb=t;
      i++;
	  if(i==MAX){break;}
   }
   fclose(pFile[0]);

   last=i-1;
   printf("\n\n last=%ld",last);

   dt=(t-t1)/double(last);
   sr=1./dt;
   printf("\n\n  damp=%g ",damp);
   printf("\n  dt=%g sec   sr=%g samples/sec \n", dt, sr);


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

	j=0;

	omegan=2.*pi*fn[j];
	omegad=omegan*sqrt(1.-pow(damp,2));
	cosd=cos(omegad*dt);
	sind=sin(omegad*dt);
	domegandt=damp*omegan*dt;
	domegan=damp*omegan;

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

}
void enter_data(void)
{

	printf("\n\n Select units:  1=English  2=metric \n");
	scanf("%d",&iunits);


	/********************* intialize frequency array *********************/

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



	if(iunits==2)
	{
		printf("\n Enter mass (kg) \n");
		scanf("%lf",&mass);
	}
	else
	{

		printf("\n Enter mass (lbm) \n");

		scanf("%lf",&mass);

		mass/=386.;

	}

// calculate stiffness...

	double omegan=tpi*fn[0];

	stiffness=mass*pow(omegan,2);


	if(iunits==2)
	{
		printf("\n stiffness = %10.4g N/m \n",stiffness);
	}
	else
	{
		printf("\n stiffness = %10.4g lbf/in \n",stiffness);
	}

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

	ccc=2.*damp*omegan*mass;

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


}

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

}