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

void factors(void);
void apply(void);
void files(void);
void check(void);

const double PI=atan2(0.,-1.);
const double TPI=2.*atan2(0.,-1); 

double uago,uvgo,udgo;

double dam,pv; 

double 	alphaAm;
double 	alphaVm;
double 	alphaDm;

double 	alphaAsig;
double 	alphaVsig;
double 	alphaDsig;

double n_med,n_sig;
double alphaVmed;
double a_med_1,a_med_2;
double a_sig_1,a_sig_2;


double ta,tb,tc,td,te,tf;


int iunits;

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

void main()
{
	printf("\n\n eds.cpp, ver 1.0,  Aug 3, 2002 \n");
	printf("\n by Tom Irvine, tomirvine@aol.com ");
	printf("\n\n This program calculates an elastic design spectrum. ");
	printf("\n Reference: Chopra, Dynamics of Structures \n");
    
	printf("\n Select units ");
	printf("\n 1= English  2= metric \n");
	scanf("%d",&iunits);



	if(iunits == 1)
	{
		printf("\n\n Enter the peak ground acceleration (G) \n");
        scanf("%lf",&uago);

		uago*=386.;

		printf("\n\n Enter the peak ground velocity (in/sec) \n");
        scanf("%lf",&uvgo);	

		printf("\n\n Enter the peak ground displacement (in) \n");
        scanf("%lf",&udgo);	
	}
	else
	{
		printf("\n\n Enter the peak ground acceleration (m/sec^2) \n");
        scanf("%lf",&uago);

		printf("\n\n Enter the peak ground velocity (m/sec) \n");
        scanf("%lf",&uvgo);	

		printf("\n\n Enter the peak ground displacement (mm) \n");
        scanf("%lf",&udgo);
		
		udgo/=1000.;
	}


	printf("\n\n Enter the damping ratio \n");
	scanf("%lf",&dam);

	
    tc = TPI*uvgo/uago;
    td = TPI*udgo/uvgo;


	printf("\n\n Cross-over period 1 = %8.3lf sec  ",tc);
	  printf("\n Cross-over period 2 = %8.3lf sec\n",td);


	ta = 1./33;
	tb = 1./8.;
	te = 10.;
	tf = 33.;


    check();
	files();
	factors();
	apply();

	printf("\n Please call the output file into your own graphics program. \n");
    printf("\n The output file has four columns. \n");

	printf("\n period   ground     median    1 sigma "); 
	printf("\n          velocity   velocity  velocity "); 

	if(iunits==1)
	{
 	printf("\n  (sec)   (in/sec)   (in/sec)   (in/sec) \n"); 
	}
	else
	{
 	printf("\n  (sec)   (m/sec)   (m/sec)   (m/sec) \n");
	}

    printf("\n Press any key to exit. \n");
	getch();
}
void factors(void)
{
	dam = log(dam*100.);

	alphaAm=3.21 - 0.68*dam;
	alphaVm=2.31 - 0.41*dam;
	alphaDm=1.82 - 0.27*dam;

	alphaAsig=4.38 - 1.04*dam;
	alphaVsig=3.38 - 0.67*dam;
	alphaDsig=2.73 - 0.45*dam;

}
void apply(void)
{
	double period;
    double oct = pow(2.,(1./12.));

	double vmed,vsig;

    period = 0.01;



	while(1)
	{


        if( period <= ta )
		{
			pv = uago *period/TPI;

			vmed=pv;
			vsig=pv;
		}
		if( period > ta && period <=tb )
		{
			pv = uago *period/TPI;

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

			n_sig = log(  alphaVsig*(uago *tb/TPI) /(uago *ta/TPI) )/log(tb/ta);
			n_med = log( 	alphaVm*(uago *tb/TPI) /(uago *ta/TPI) )/log(tb/ta);

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

            vsig= (uago *ta/TPI) * pow( (period/ta), n_sig );
			vmed= (uago *ta/TPI) * pow( (period/ta), n_med );

		}
		if( period > tb && period <=tc )
		{
			pv = uago *period/TPI;
			vmed=pv*alphaVm;
			vsig=pv*alphaVsig;
		}
		if( period > tc && period <=td )
		{
			pv = uvgo;
			vmed=pv*alphaVm;
			vsig=pv*alphaVsig;
		}
		if( period > td && period <=te )
		{
			pv= udgo * TPI/period;
			vmed=pv*alphaVm;
			vsig=pv*alphaVsig;
		}
		if( period > te && period <=tf )
		{
			pv= udgo * TPI/period;


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

			n_sig = log(  (udgo *TPI/tf) /(alphaVsig*udgo *TPI/te) )/log(tf/te);
			n_med = log(  (udgo *TPI/tf) /(alphaVm  *udgo *TPI/te) )/log(tf/te);

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

            vsig= (alphaVsig*udgo *TPI/te) * pow( (period/te), n_sig );
			vmed= (alphaVm  *udgo *TPI/te) * pow( (period/te), n_med );

		}
        if( period > tf )
		{
			pv= udgo * TPI/period;
			vmed=pv;
			vsig=pv;
		}


		fprintf(pFile[0]," %10.4g  %10.4g  %10.4g  %10.4g\n",period,pv,vmed,vsig);

		period*=oct;

		if(period > 52. ){break;}
	}
	fclose(pFile[0]);

}
void files(void)
{
	printf("\n\n Enter the output filename. \n");
	scanf("%s",filename[0]);
	pFile[0]=fopen(filename[0],"w");

}
void check(void)
{
	if( tc < tb )
	{
		tb=tc;
	}

	if(td > te )
	{
		te=tb;
	}

	if(   td < tc   ||   td > tf   ||   tc < ta   )
	{
		printf("\n Input error. \n");
		
		printf("\n Press any key to exit. \n");
	    getch();
	    exit(1);
	}

}