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


#define MAXCHECK \
if(y>dmax){dmax=y; tmax=t;}\
if(y<dmin){dmin=y; tmin=t;}\
if(v>vmax){vmax=v; tvmax=t;}\
if(v<vmin){vmin=v; tvmin=t;}\
if(acc>amax){amax=acc; tamax=t;}\
if(acc<amin){amin=acc; tamin=t;}\

void ending(void);
void engine(void);
void printout(void);
void timecalc(void);

double U1,V1;
double U2,V2;

double m,k,damp;
double y,yinitial,ypinitial;

double dmax,dmin;
double vmax,vmin;
double amax,amin;

double dt,dt2,t;
double fn;

double ee, cwdt, swdt;

double a1,a2,a3,acc;
double dur;
double f0;

double yone,y2,y3,y4,y5;
double v,f;
double ypi,yi;

double dy1,dy2,dy3,dy4,dy5;

double domegan;

double omegat,omegant,omegadt;
double omega2,omegan2,omegad2;


double omega, omegan, omegad;

double tmin,tmax,tvmin,tvmax,tamin,tamax;

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


long i,ijk,j,jnum;
long idamp,last,limit,nt;
long iopt,ioct;
long iu;


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

int main()
{

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

	printf("\n fhsine.cpp  version 1.4  January 16, 2011 \n");
	printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

	printf("\n This program calculates the displacement time history");
	printf("\n response of a single-degree-of-freedom system subjected");
	printf("\n to a half-sine force shock.");

	printf("\n\n Assume zero initial displacement and zero initial velocity.");

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

	strcpy(filename[1],"disp.out");

	strcpy(filename[2],"velox.out");

	strcpy(filename[3],"accel.out");

	strcpy(filename[4],"force.out");

	strcpy(filename[5],"dsrs.out");

	strcpy(filename[6],"vsrs.out");

	strcpy(filename[7],"asrs.out");

	printf("\n\n");
	for(i=1; i<=7; i++)
	{
		pFile[i]=fopen(filename[i],"w");

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

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

	printf("\n\n Enter analysis option \n");
	printf("\n 1=time history response  2=SRS \n");
	scanf("%ld",&iopt);

	if(iopt==2)
	{
		printf("\n\n Enter frequency spacing option \n");
		printf("\n 1= 1/24 oct  2= constant 10 Hz \n");
		scanf("%ld",&ioct);
	}

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

	if(iu==1)
	{
        printf("\n\n Enter mass (lbm) \n");
        scanf("%lf",&m);
        m/=386.;
	}
    else
    {
        printf("\n\n Enter mass (kg) \n");
        scanf("%lf",&m);
    }



	if(iopt==1)
	{
		printf("\n\n Enter natural frequency (Hz) \n");
		scanf("%lf",&fn);
	}
	else
	{
		printf("\n\n Enter starting frequency (Hz) \n");
		scanf("%lf",&fn);
	}

   	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 amplfication factor Q (typically 10) ");
      scanf("%lf",&damp);
	  damp = 1./(2.*damp);
	}

    if(iu==1)
    {
        printf("\n\n Enter amplitude of half-sine force shock (lbf) \n");
    }
    else
    {
        printf("\n\n Enter amplitude of half-sine force shock (N) \n");
    }
	scanf("%lf",&f0);


	printf("\n\n Enter duration (sec) \n");
	scanf("%lf",&dur);

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

    if(iopt==1)
	{limit=1;}
	else
	{limit=200;}


	dt2 = dur/32.;

    engine();

	ending();

}

void ending(void)
{
	if( iopt==1)
	{
        if(dmax< fabs(dmin)){dmax=fabs(dmin); tmax=tmin;}

        if(iu==1)
        {
            printf("\n\n stiffness k = %12.4g lbf/in",k);

            printf("\n\n maximum displacement = %12.4g in at %8.4f sec",dmax, tmax);
            printf("\n minimum displacement = %12.4g in at %8.4f sec",dmin, tmin);

            printf("\n\n maximum absolute nondimensional displacement kx/F = %12.4f ",k*dmax/f0);

            printf("\n\n\n maximum velocity = %12.4g in/sec at %8.4f sec",vmax, tvmax);
            printf("\n minimum velocity = %12.4g in/sec at %8.4f sec",vmin, tvmin);

            printf("\n\n\n maximum acceleration = %12.4g G at %8.4f sec",amax/386., tamax);
            printf("\n minimum acceleration = %12.4g G at %8.4f sec",amin/386., tamin);

            printf("\n\n\n Output Files: \n");
            printf("\n force.out - force input                 (lbf)  vs time(sec)");
            printf("\n disp.out  - displacement response       (in) vs time(sec)");
            printf("\n velox.out - velocity     response   (in/sec) vs time(sec)");
            printf("\n accel.out - acceleration response  (G) vs time(sec)");
        }
	    else
	    {
            printf("\n\n stiffness k = %12.4g N/m",k);

            printf("\n\n maximum displacement = %12.4g mm at %8.4f sec",1000.*dmax, tmax);
            printf("\n minimum displacement = %12.4g mm at %8.4f sec",1000.*dmin, tmin);

            printf("\n\n maximum absolute nondimensional displacement kx/F = %12.4f ",k*dmax/f0);

            printf("\n\n\n maximum velocity = %12.4g mm/sec at %8.4f sec",1000.*vmax, tvmax);
            printf("\n minimum velocity = %12.4g mm/sec at %8.4f sec",1000.*vmin, tvmin);

            printf("\n\n\n maximum acceleration = %12.4g m/sec^2 at %8.4f sec",amax, tamax);
            printf("\n minimum acceleration = %12.4g m/sec^2 at %8.4f sec",amin, tamin);

            printf("\n\n\n Output Files: \n");
            printf("\n force.out - force input                 (N)  vs time(sec)");
            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)");
	    }

    }
	else
	{
		printf("\n\n\n Output File: \n");


        if(iu==1)
        {
            printf("\n dsrs.out - displacement SRS:  fn(Hz)  Peak Pos Disp(in)      Peak Neg Disp(in)");
            printf("\n vsrs.out - velocity     SRS:  fn(Hz)  Peak Pos Vel (in/sec)  Peak Neg Vel (in/sec)");
            printf("\n asrs.out - acceleration SRS:  fn(Hz)  Peak Pos Acc (G) Peak Neg Acc (G)");
        }
        else
        {
            printf("\n dsrs.out - displacement SRS:  fn(Hz)  Peak Pos Disp(mm)      Peak Neg Disp(mm)");
            printf("\n vsrs.out - velocity     SRS:  fn(Hz)  Peak Pos Vel (mm/sec)  Peak Neg Vel (mm/sec)");
            printf("\n asrs.out - acceleration SRS:  fn(Hz)  Peak Pos Acc (m/sec^2) Peak Neg Acc (m/sec^2)");
        }


//		printf("\n nda_srs.out - acceleration SRS:  fn(Hz)  Peak Pos Acc (m/sec^2) Peak Neg Acc (m/sec^2)");
	}

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

void engine(void)
{
	omega = 2.*pi/(2.*dur);
	omega2=pow(omega,2.);

	for(ijk=1; ijk<=limit; ijk++)
	{

		if(fn>20000.)
		{
			break;
		}

		dt = 1./(32.*fn);

		if(dt>dt2){dt=dt2;}

		nt = 8*long( ( (1./fn) + dur) /dt );

		omegan=tpi*fn;

		k = m*pow(omegan,2.);

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


		yinitial = 0.;
		ypinitial = 0.;

		omegan2=pow(omegan,2.);
		omegad2=pow(omegad,2.);

        domegan=damp*omegan;

        U1=domegan*omega;
        V1=(omega2-omegan2);

        U2=domegan*omegad;
        V2=(-1.+2.*pow(damp,2.))*omegan2;

		a1=(omegan2*f0/k);
		a1/=pow( V1,2.)+pow((2.*damp*omega*omegan),2.);

		a2=a1*omega/omegad;

		a3=a2;


		dmax=-1.0e+12;
		dmin=+1.0e+12;

		vmax=-1.0e+12;
		vmin=+1.0e+12;

		amax=-1.0e+12;
		amin=+1.0e+12;

		tmax=0.;
		tmin=0.;

		tvmax=0.;
		tvmin=0.;

		tamax=0.;
		tamin=0.;


		timecalc();

		printout();


	} // close frequency

}

void timecalc(void)
{
		double cat,sat,doo;

		double tt;

// solve for yi,ypi  at t=dur

        t=dur;

 		omegat=omega*t;
		omegadt=omegad*t;
		omegant=omegan*t;

		ee=exp(-domegan*t);
		cwdt=cos(omegadt);
		swdt=sin(omegadt);

		doo = (domegan/omegad);
		cat = cos(omegat);
		sat = sin(omegat);

		y=0.;

		yone=yinitial*ee*(cwdt+doo*swdt);
		y2=ypinitial*ee*swdt/omegad;

		y3=a1*(-2.*U1*cat-V1*sat);
		y4=a2*ee*(2.*U2*cwdt);
		y5=a3*ee*(omega2 + V2)*swdt;

		dy1=-domegan*yone + omegad*yinitial*ee*(-swdt+doo*cwdt);
		dy2=-domegan*y2 + ypinitial*ee*cwdt/omegad;

		dy3=omega*a1*(+2.*U1*sat-V1*cat);

		dy4=-domegan*y4-omegad*a2*ee*(2.*U2*swdt);
		dy5=-domegan*y5+omegad*a3*ee*(omega2 + V2)*cwdt;

		y=yone+y2+y3+y4+y5;
		v=dy1+dy2+dy3+dy4+dy5;

        yi=y;
        ypi=v;


		for(i=0; i<nt; i++)
		{
			t=dt*i;

			omegat=omega*t;
			omegadt=omegad*t;
			omegant=omegan*t;

			ee=exp(-domegan*t);
			cwdt=cos(omegadt);
			swdt=sin(omegadt);

			doo = (domegan/omegad);
			cat = cos(omegat);
			sat = sin(omegat);


			if(t <= dur)
			{
				f=f0*sat;

				y=0.;

				yone=yinitial*ee*(cwdt+doo*swdt);
				y2=ypinitial*ee*swdt/omegad;

				y3=a1*(-2.*U1*cat-V1*sat);
				y4=a2*ee*(2.*U2*cwdt);
				y5=a3*ee*(omega2 + V2)*swdt;

				dy1=-domegan*yone + omegad*yinitial*ee*(-swdt+doo*cwdt);
				dy2=-domegan*y2 + ypinitial*ee*cwdt/omegad;

				dy3=omega*a1*(+2.*U1*sat-V1*cat);

				dy4=-domegan*y4-omegad*a2*ee*(2.*U2*swdt);
				dy5=-domegan*y5+omegad*a3*ee*(omega2 + V2)*cwdt;

				y=yone+y2+y3+y4+y5;
				v=dy1+dy2+dy3+dy4+dy5;


				acc= (omegan2/k)*f -2.*domegan*v -omegan2*y;

				MAXCHECK

			}
			else
			{
				f=0.;

				tt=t-dur;

				omegat=omega*tt;
				omegadt=omegad*tt;
				omegant=omegan*tt;

				ee=exp(-domegan*tt);
				cwdt=cos(omegadt);
				swdt=sin(omegadt);

				yone=yi *ee*(cwdt+doo*swdt);
				y2=ypi*ee*swdt/omegad;

				dy1=-domegan*yone + omegad*yi*ee*(-swdt+doo*cwdt);
				dy2=-domegan*y2 + ypi*ee*cwdt;

				y=yone+y2;
				v=dy1+dy2;

				acc= -2.*domegan*v -omegan2*y;

				MAXCHECK

			}

			if(iopt==1)
			{
			    if(iu==1)
			    {
                    fprintf(pFile[1],"%14.7e \t %14.7e\n",t,y);
                    fprintf(pFile[2],"%14.7e \t %14.7e\n",t,v);
                    fprintf(pFile[3],"%14.7e \t %14.7e\n",t,acc/386.);
                    fprintf(pFile[4],"%14.7e \t %14.7e\n",t,f);
			    }
			    else
			    {
                    fprintf(pFile[1],"%14.7e \t %14.7e\n",t,1000.*y);
                    fprintf(pFile[2],"%14.7e \t %14.7e\n",t,1000.*v);
                    fprintf(pFile[3],"%14.7e \t %14.7e\n",t,acc);
                    fprintf(pFile[4],"%14.7e \t %14.7e\n",t,f);
			    }
			}


		} // close time
}

void printout(void)
{
		if(iopt==2)
		{
				fprintf(pFile[5],"%14.7e \t %14.7e \t %14.7e\n",fn,dmax*1000.,fabs(dmin)*1000.);
				fprintf(pFile[6],"%14.7e \t %14.7e \t %14.7e\n",fn,vmax*1000.,fabs(vmin)*1000.);
				fprintf(pFile[7],"%14.7e \t %14.7e \t %14.7e\n",fn,amax,fabs(amin));

				if(ioct==1)
				{
					fn*=pow(2.,(1./12.));
				}
				else
				{
					fn+=10.;
				}

		}
}