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


#define MAXCHECK \
if(acc>amax){amax=acc; tamax=t;}\
if(acc<amin){amin=acc; tamin=t;}\
if(rd>rd_max){rd_max=rd; trmax=t;}\
if(rd<rd_min){rd_min=rd; trmin=t;}\




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

void transition(void);


const double mm_per_inch=25.4;

double rd_initial;

double rv_initial;

double damp;
double y;

double ain;
double amax,amin;

double dt,dt2,t;
double fn;

double ee, cwdt, swdt;

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

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

double fff;

double vT,zT;

double tamin,tamax;
double trmin,trmax;


double rd,rd1,rd2,rd3;

double rv,rv1,rv2,rv3;


double rd_max,rd_min;

double abefore;

double rscale;


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


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


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

void main()
{

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

	printf("\n ahsine.cpp  version 1.7  September 19, 2009\n");
	printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

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

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

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


	strcpy(filename[1],"accel.in");
	strcpy(filename[2],"accel.out");
	strcpy(filename[3],"asrs.out");
	strcpy(filename[4],"rd.out");

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




	printf("\n\n");
	for(i=1; i<=5; 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/12 oct  2= constant 10 Hz \n");

		scanf("%ld",&ioct);

	}


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



	printf("\n\n Enter acceleration unit. \n 1= G   2= m/sec^2\n");

	scanf("%ld",&iunit);


	if(iunit==1)
	{
		printf("\n\n Enter amplitude of half-sine base input shock (G) \n");
	}
	else
	{
		printf("\n\n Enter amplitude of half-sine base input shock (m/sec^2) \n");	

	}
		
	rscale = 386.;
	

	scanf("%lf",&amp);

	if(iunit==2)
	{
		amp/=9.81;
	}

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


		printf("\n\n\n maximum acceleration = %12.4g G at %8.4f sec",amax, tamax);
		printf("\n minimum acceleration = %12.4g G at %8.4f sec",amin, tamin);
		
		printf("\n\n maximum acceleration = %12.4g m/sec^2 ",amax*9.81);
		printf("\n minimum acceleration = %12.4g m/sec^2 ",amin*9.81);
	
	
		printf("\n\n\n\n maximum relative displacement = %12.4g inch at %8.4f sec",rd_max*rscale, trmax);
		printf("\n minimum relative displacement = %12.4g inch at %8.4f sec",rd_min*rscale, trmin);

	
		printf("\n\n maximum relative displacement = %12.4g mm ",rd_max*rscale*mm_per_inch);
		printf("\n minimum relative displacement = %12.4g mm ",rd_min*rscale*mm_per_inch);

		

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


		if(iunit==1)
		{

	    printf("\n accel.in  - base input acceleration:  time(sec)  accel(G)");
	    printf("\n accel.out - response   acceleration:  time(sec)  accel(G)");

	    printf("\n    rd.out - relative displacement  :  time(sec)  disp(inch)");

		}

		else

		{

	    printf("\n accel.in  - base input acceleration:  time(sec)  accel(m/sec^2)");

	    printf("\n accel.out - response   acceleration:  time(sec)  accel(m/sec^2)");

	    printf("\n    rd.out - relative displacement  :  time(sec)  disp(mm)");

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


		if(iunit==1)
		{
		printf("\n   asrs.out - acceleration  SRS:  fn(Hz)  Pos Acc(G)      Neg Acc(G)");

		printf("\n rd_srs.out - relative disp SRS:  fn(Hz)  Pos Disp(inch)  Neg Disp(inch) \n");
		}

		else
		{

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

		printf("\n rd_srs.out - relative disp SRS:  fn(Hz)  Pos Disp(mm)    Neg Disp(mm) \n");		

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

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


		fff = ( (1./fn) + dur) /dt;


		nt = 12*long( fff );

		omegan=tpi*fn;

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

		rd_initial = 0.;
		rv_initial = 0.;

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


		domegan=damp*omegan;


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


		rd_max=-1.0e+12;
		rd_min=+1.0e+12;

		tamax=0.;
		tamin=0.;


		timecalc();

		printout();


	} // close frequency
	
}

void timecalc(void)
{
		double cat,sat;

		double a1,a2;

		double b1,b2;

		double c1,c2;

		double den;



		a1=rd_initial;

		a2=(rv_initial+domegan*rd_initial)/omegad;


        b1=2.*damp*omega*omegan;

		b2=(omega2-omegan2);



		c1=2.*domegan*omegad;

		c2=omega2-omegan2*(1.-2*pow(damp,2.));



        den=pow((omega2-omegan2),2.)+pow((2.*damp*omega*omegan),2.);



        transition();


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

			if(t <= dur)
			{



				omegat=omega*t;

				omegadt=omegad*t;

				omegant=omegan*t;

				domegant=domegan*t;



				ee=exp(-domegant);

				cwdt=cos(omegadt);

				swdt=sin(omegadt);



				cat = cos(omegat);

				sat = sin(omegat);



				ain= amp*sin(omegat);

				

				rd1=ee         *(a1*cwdt + a2*swdt);

				rd2=(amp/den)  *(b1*cat  + b2*sat );

				rd3=-((amp/den)*ee*omega/omegad)*(c1*cwdt + c2*swdt);

				

				rd=rd1+rd2+rd3;





				rv1=-domegan*rd1;

                rv1+=omegad*ee*(-a1*swdt + a2*cwdt);



				rv2=omega*(amp/den)  *(-b1*sat  + b2*cat );				



				rv3=-domegan*rd3;

                rv3+=-((amp/den)*ee*omega)*(-c1*swdt + c2*cwdt);



				rv=rv1+rv2+rv3;



//				zT=rd;

//				vT=rv;
			}
			else
			{

				ain=0.;



				omegat=omega*(t-dur);

				omegadt=omegad*(t-dur);

				omegant=omegan*(t-dur);

				domegant=domegan*(t-dur);



				ee=exp(-domegant);

				cwdt=cos(omegadt);

				swdt=sin(omegadt);



				a1=zT;

				a2=(vT+domegan*zT)/omegad;



                rd = ee *(a1*cwdt + a2*swdt);



				rv= -domegan*rd;

				rv+= omegad*ee *(-a1*swdt + a2*cwdt);


			}



			acc= -omegan2*rd  - 2.*domegan*rv;

 //           printf(" t=%8.4g acc=%8.4g  rd=%8.4g  rv=%8.4g \n",t,acc,rd,rv); 


			MAXCHECK

			if(iopt==1)
			{
				if(iunit==1)
				{

					fprintf(pFile[1],"%14.7e \t %14.7e\n",t,ain);	
					fprintf(pFile[2],"%14.7e \t %14.7e\n",t,acc);

					fprintf(pFile[4],"%14.7e \t %14.7e\n",t,rd*rscale);
				}
				else
				{
					fprintf(pFile[1],"%14.7e \t %14.7e\n",t,ain*9.81);	
					fprintf(pFile[2],"%14.7e \t %14.7e\n",t,acc*9.81);

					fprintf(pFile[4],"%14.7e \t %14.7e\n",t,rd*rscale*mm_per_inch);
				}
			}



			if( i > 0.92*nt)

			{

				if( (abefore*acc)<=0 )

				{

					break;

				}

			}

			abefore=acc;




		} // close time
}

void printout(void)
{

		if(iopt==2)
		{
				if(iunit==1)
				{

					fprintf(pFile[3],"%14.7e \t %14.7e \t %14.7e\n",fn,amax,fabs(amin));
					fprintf(pFile[5],"%14.7e \t %14.7e \t %14.7e\n",fn,rd_max*rscale,fabs(rd_min)*rscale);
				}
				else
				{
					fprintf(pFile[3],"%14.7e \t %14.7e \t %14.7e\n",fn,amax*9.81,fabs(amin*9.81));
					fprintf(pFile[5],"%14.7e \t %14.7e \t %14.7e\n",fn,rd_max*rscale*mm_per_inch,fabs(rd_min)*rscale*mm_per_inch);
				}


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

				}

				else
				{

					fn+=10.;

				}


		}
}

void transition(void)

{

	double cat,sat;

	double a1,a2;

	double b1,b2;

	double c1,c2;

	double den;



	a1=rd_initial;

	a2=(rv_initial+domegan*rd_initial)/omegad;



	b1=2.*damp*omega*omegan;

	b2=(omega2-omegan2);



	c1=2.*domegan*omegad;

	c2=omega2-omegan2*(1.-2*pow(damp,2.));



	den=pow((omega2-omegan2),2.)+pow((2.*damp*omega*omegan),2.);



	t=dur;



	omegat=omega*t;

	omegadt=omegad*t;

	omegant=omegan*t;

	domegant=domegan*t;



	ee=exp(-domegant);

	cwdt=cos(omegadt);

	swdt=sin(omegadt);



	cat = cos(omegat);

	sat = sin(omegat);



	ain= amp*sin(omegat);

				

	rd1=ee         *(a1*cwdt + a2*swdt);

	rd2=(amp/den)  *(b1*cat  + b2*sat );

	rd3=-((amp/den)*ee*omega/omegad)*(c1*cwdt + c2*swdt);

				

	rd=rd1+rd2+rd3;



	rv1=-domegan*rd1;

	rv1+=omegad*ee*(-a1*swdt + a2*cwdt);



	rv2=omega*(amp/den)  *(-b1*sat  + b2*cat );				



	rv3=-domegan*rd3;

	rv3+=-((amp/den)*ee*omega)*(-c1*swdt + c2*cwdt);



	rv=rv1+rv2+rv3;



	zT=rd;

	vT=rv;

}