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


#include <time.h>
#include <sys/types.h>
#include <sys/timeb.h>

#define MINF 0.1

void elapsed(void);



time_t ltime;

time_t ltimen;




#define CORE \
{\
	x[j]=	 a1[j]*xb[j]\
			+a2[j]*xbb[j]\
			+b1[j]*y\
			+b2[j]*yb[j]\
            +b3[j]*ybb[j];\
 \
	if(x[j] > xmax[j])\
	{\
		xmax[j]=x[j];\
	}\
	if(x[j] < xmin[j])\
	{\
		xmin[j]=x[j];\
	}\
 \
	xbb[j]=xb[j];\
	xb[j]= x[j];\
	ybb[j]= yb[j];\
	yb[j]= y;\
 \
}\



#define MAXTRIALS 20000

#define K100 100000





#define NUM 125

#define MAX 125





void inputspec(void);

void intlog(void);



void synth1(void);

void synth2(void);

void synth3(void);





void print_th(void);

void srs(void);

void scale(void);

void gen_time(void);

void max_param(void);

void moreinput(void);

void srs_error(void);

void srs_pre(void);

void rank(void);

void integrate(void);

void zero(void);

void print_srs(void);

void rankfunctions(void);

void weights(void);





FILE *pFile[15];

char filename[15][FILENAME_MAX];





long i,ichoice,idamp,igen,ikk,inn,iseed,istrat,iunit,iwin,j,jk;



long limit;



long nhs[MAXTRIALS][MAX],nspec,nt,ntrials,num,nv;

long rntrials;





float ym[MAXTRIALS],vm[MAXTRIALS],dm[MAXTRIALS],em[MAXTRIALS],im[MAXTRIALS];

float sym[MAXTRIALS],svm[MAXTRIALS],sdm[MAXTRIALS],sem[MAXTRIALS],sim[MAXTRIALS];



long yrank[MAXTRIALS],vrank[MAXTRIALS],drank[MAXTRIALS],erank[MAXTRIALS],irank[MAXTRIALS];

long pyrank[MAXTRIALS],pvrank[MAXTRIALS],pdrank[MAXTRIALS],perank[MAXTRIALS],pirank[MAXTRIALS];



float nrank[MAXTRIALS];





float upper[NUM];



float alpha,beta;



float a1[NUM],a2[NUM];

float b1[NUM],b2[NUM],b3[NUM];





double local_record;



float local_amp[MAX];



float a,amp[MAXTRIALS][MAX],area,az;

float b;



float cosd;

float damp;

float dmax,domegadt;



float dt,dth,dur;

float error;



float f[MAX];

float fa,fb,ff;

float fn[NUM];

float fr[MAX];



float frlast;

float arlast;



float imax,irror;





float octave,omega,omegad;



float omegaf[MAX];



float over_period[NUM];

float onep5_period[NUM];



float r[MAX];



float sind;

float s[MAX],spec[MAX],sr;



float amp_start[MAX];



float t,td[MAXTRIALS][MAX],tt;

float th[K100];

float t1, tb;

float tmax[NUM],tmin[NUM];



float vmax,vth[K100];



float wavelet[NUM][K100];



float xb[NUM], xbb[NUM];

float x[NUM],xmax[NUM],xmin[NUM];

float xxmin[NUM],xxmax[NUM];



float y,ymax,yone,ytwo;

float yb[NUM],ybb[NUM];



float xrmin,irmin;



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

const double tpi=2.*pi;



const double max=32767.;





double exponent;



const long niter=60;



double tstart;



double iw,ew,dw,vw,aw;



char aunit[10],dunit[10],vunit[10];



char ff_wt[20];



char stype[20];


int algorithm;



int main()
{



	printf("\n\n wavelet_synth.cpp,  ver 4.5, December 8, 2011 \n\n by Tom Irvine \n Email: tomirvine@aol.com \n\n");

	printf("\n This program synthesizes a time history using wavelets to satisfy ");

	printf("\n a shock response spectrum (SRS) specification. \n");



	printf("\n The program also optimizes the time history to yield the lowest overall \n error, acceleration, velocity, and displacement. \n");

    printf("\n The optimization is performed via trial-and-error.\n");



	printf("\n Select SRS algorithm:");

	printf("\n 1=Kelly-Richman  2=Smallwood \n");



    scanf("%d",&algorithm);





	zero();



    inputspec();



	intlog();



	moreinput();



	time( &ltime );





	sr=float(10.*fr[num-1]);



	printf("\n\n Recommend sample rate >= %8.4g samples/sec\n",sr);

    printf("\n Enter sample rate \n");

	scanf("%f",&sr);



	if(sr < 4.*fr[num-1])

	{

		printf("\n Error: sample rate is too low \n");

		exit(1);

	}





	dt=float(1./sr);

	dth=float(dt/2.);



	int idur;



	printf("\n\n Select condition  ");

	printf("\n 1=limit time history points   2=set duration \n");



	scanf("%d",&idur);



	if(idur==2)

	{



		printf("\n\n Enter duration (sec) \n");

		printf(" (Recommend %6.3lf or greater) \n",2.0/f[0]);



		scanf("%f",&dur);



	    nt=long(dur/dt);



		if(nt>K100)

		{

			printf("\n\n Warning: duration reduced. \n");



			nt=K100;

		}

	}

	else

	{



		printf("\n\n Enter the number of points \n");



		scanf("%ld",&nt);



		if(nt>K100)

		{

			printf("\n\n Warning: duration reduced. \n");



			nt=K100;

		}



		dur=nt*dt;



	}



	printf("\n dt=%9.4g sec   dur=%8.4f sec  sr=%9.4g sample/sec  nt=%ld \n",dt,dur,sr,nt);





	tstart=float(ltime);





	if(dur < 1.5/f[0])

	{

		printf("\n\n Error: duration is too short.\n\n");

	}

    else

    {



		srs_pre();



		float errorbefore=0.;



		long ijk=0;



		pFile[1]=fopen(filename[1], "w");

		if(nt < K100)

		{



				printf("\n           Peak        Peak         Peak         ");

				printf("\n Trial     Accel       Vel          Disp      T.Error      I.Error");



				fprintf(pFile[1],"\n           Peak        Peak         Peak         ");

				fprintf(pFile[1],"\n Trial     Accel       Vel          Disp      T.Error      I.Error");





				if(iunit==1)

				{

					          printf("\n           (G)         (in/sec)     (in)       \n");

					fprintf(pFile[1],"\n           (G)         (in/sec)     (in)       \n");

				}

				if(iunit==2)

				{

					          printf("\n           (G)         (m/sec)      (mm)       \n");

  					fprintf(pFile[1],"\n           (G)         (m/sec)      (mm)       \n");

				}

				if(iunit==3)

				{

					          printf("\n           (m/sec^2)   (m/sec)      (mm)       \n");

  					fprintf(pFile[1],"\n           (m/sec^2)   (m/sec)      (mm)       \n");



				}





            ichoice = 1;



			for(i=0; i<nspec; i++)

			{

				omegaf[i]=float(tpi*f[i]);



				over_period[i]=float(1.0/f[i]);

			   onep5_period[i]=float(1.5/f[i]);



			   local_amp[i]=0.;

			}

			limit=long(1.2/(dt*f[0]));





			for(inn=0; inn< ntrials; inn++)

			{

				local_record = 1.0e+99;



				for(i=1; i<iseed; i++)

				{

					a=float(rand()/max);

				}



				exponent = 0.5;



				if( rand()/max < 0.4 )

				{

					exponent = ( 0.40 + 0.1*rand()/max );

				}



				if(ichoice==1){synth3();}

				if(ichoice==2){synth3();}

				if(ichoice==3){synth3();}

				if(ichoice==4){synth3();}

				if(ichoice==5){synth3();}

				if(ichoice==6){synth1();}



				ichoice++;



				igen=inn;


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

					upper[i] = float(nhs[igen][i]/(2.*f[i]));

				}


				for(nv=0; nv< niter; nv++)
				{

				    gen_time();

					srs();



					srs_error();





					if( irror < local_record )   // end criteria

					{

						for(i=0; i<nspec; i++)

						{

							local_amp[i]=amp[igen][i];



							local_record = irror;

						}

					}



					if(nv>= 9 && irror >= errorbefore){break;}



					errorbefore=irror;



					scale();



				}



				for(i=0; i<nspec; i++)

				{

					amp[igen][i]=local_amp[i];

				}



				max_param();



				irror=float(local_record);





				          printf(" %ld %10.2f %12.3f %12.3f %11.2f %10.2f  %ld\n",inn,ymax,vmax,dmax,20.*error,20.*irror,nv);

				fprintf(pFile[1]," %ld %10.2f %12.3f %12.3f %11.2f %10.2f  %ld\n",inn,ymax,vmax,dmax,20.*error,20.*irror,nv);



				if(ichoice==7){ichoice=1;}





			    sym[inn]=float(fabs(ymax));

			    svm[inn]=float(fabs(vmax));

				sdm[inn]=float(fabs(dmax));

				sem[inn]=float(fabs(error));

				sim[inn]=float(fabs(irror));



//				printf("\n inn=%ld dminn=%12.4f dmax=%12.4f \n",inn,dm[inn],dmax);



				ijk++;

/*

				if(ijk==1000 && inn < (ntrials-7) )

				{

                    for(jk=0; jk< ntrials; jk++)

					{

						ym[jk]=sym[jk];

						vm[jk]=svm[jk];

						dm[jk]=sdm[jk];

						em[jk]=sem[jk];

						im[jk]=sim[jk];

					}



                    rntrials=inn;



					printf("\n\n intermediate \n");

					rankfunctions();



					ijk=0;

				}

*/

			}

			for(jk=0; jk< ntrials; jk++)

			{

						ym[jk]=sym[jk];

						vm[jk]=svm[jk];

						dm[jk]=sdm[jk];

						em[jk]=sem[jk];

						im[jk]=sim[jk];

			}



			rntrials=ntrials;



			rankfunctions();



		}

		else

		{

			printf("\n\n Error: too many samples.\n\n");

		}

	}

//	pFile[1]=fopen(filename[1], "w");







	elapsed();



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

   getch();

}

void inputspec(void)

{

	printf("\n\n\n The input SRS specification file must have two columns: \n\n   natural frequency(Hz) and SRS accel(G) ");

	printf("\n\n The format is free, but no header lines allowed.\n");

	printf("\n\n Please enter the 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]);

		fclose(pFile[0]);

	}

	else
	{

		printf(" File: %s opened. \n", filename[0]);

	}


	i=0;

	while( fscanf(pFile[0], "%f %f", &fr[i], &r[i]) > 0 )

	{

        frlast=fr[i];

		arlast= r[i];



	    i++;

    }

    num=i;



    fclose(pFile[0]);



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

	strcpy(filename[2],"wavelet");

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

	strcpy(filename[4],"srs");

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

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



	strcpy(filename[7],"velox");

	strcpy(filename[8],"disp");





}

void intlog(void)

{



//  Calculate slopes between input points





       for(i=0; i<= num-2; i++)

	   {

          a=float(log(r[i+1])-log(r[i]));

          b=float(log(fr[i+1])-log(fr[i]));

          s[i]=a/b;



//		  printf("%lf\n",s[i]);

       }



//       fprintf(pFile[1],"%15.6e %15.6e\n",f[0],psd[0] );



//  Interpolate





	   f[0]= fr[0];

	   fa=fr[0];

       fb=fr[0];



	   spec[0]=r[0];

	   amp_start[0]=float(spec[0]/16.);



	   i=1;





       long ioct;



	   printf("\n\n Enter octave spacing. \n 1= 1/3   2= 1/6   3= 1/12 \n");

	   scanf("%ld",&ioct);



	   octave=float(1./3.);



	   if(ioct==2)

	   {

		  octave=float(1./6.);

	   }

	   if(ioct==3)

	   {

		  octave=float(1./12.);

	   }





	   while(1)

	   {

		  ff=float(pow(2.,octave)*fb);

		  fb=ff;



		  f[i]=ff;



		  if(ff>fr[num-1]){break;}



		  if( ff >= fr[0])

          {



				for( j=1; j<num; j++)

				{

					if(ff < fr[j])

					{

						f[i]=ff;



						az=float(log10(r[j-1]));

						az+=float(s[j-1]*(log10(ff)-log10(fr[j-1])));

						spec[i]=float(pow(10.,az));

						amp_start[i]=float(spec[i]/16.);

						break;

					}

					if(ff == fr[j])

					{

						f[i]=ff;

						spec[i]=r[j];

						amp_start[i]=float(spec[i]/16.);

						break;

					}

				}

		  }

		  i++;

       }

	   nspec=i;



	   if(frlast > f[nspec-1])

	   {

          nspec++;



		     f[nspec-1]=frlast;

		  spec[nspec-1]=arlast;

          amp_start[nspec-1]=float(arlast/16.);

	   }



	   if(nspec > NUM)

	   {

		   printf("\n Warning: number of specification points reduced.");

		   nspec=NUM;

	   }

/*

	  for( j=0; j<nspec; j++)

	  {

		  printf(" amp_start[%ld]=%8.4g \n",j,amp_start[j]);

	  }

*/



}

void synth1(void)  // random
{

	strcpy(stype,"synth1");

	int nflag=0;

	while(nflag==0)
	{

		nflag=1;


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

			amp[inn][i]= amp_start[i];


			if( (rand()/max) < 0.5 )
			{
				amp[inn][i]*=-1.;

			}



			td[inn][i]=float((rand()/max)*(dur - onep5_period[i]));

		}

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


			nhs[inn][i]= -1 + 2*long(60.*rand()/max);


			if( nhs[inn][i] < 3 )
			{
				nhs[inn][i]=3;

			}



			while(1)
			{

				if( (float(nhs[inn][i])/(2.*f[i]))+td[inn][i] >= dur )
				{
					nhs[inn][i]-=2;

				}
				else
				{
					break;
				}
			}
		}

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

			if( float(nhs[inn][i]) < 3 )
			{

				printf("s1: error: i=%ld  nhs= %9.4g  td=%9.4g   dur=%9.4g \n",i,nhs[inn][i],td[inn][i],dur);
				exit(1);

				nflag=0;

			}
		}
	}


}

void synth2(void)  // max nhs  (currently unused)
{



}

void synth3(void)  // reverse
{

	strcpy(stype,"synth3");


	int nflag=0;


	while(nflag==0)
	{

		nflag=1;



		alpha=float(1.1*rand()/max);

		beta=float(rand()/max);



		float gamma=float(rand()/max);



		long nlimit = (2*long(65.*beta))-1;

		if(nlimit<5){nlimit=5;}



		if(inn< long(0.10*ntrials)){alpha=float(2./3.);}



		td[inn][nspec-1]=0.;

		nhs[inn][nspec-1]= nlimit;

		amp[inn][nspec-1]= amp_start[nspec-1];



		for(i=(nspec-2); i>=0; i--)
		{

			amp[inn][i]= amp_start[i];



			if(gamma<0.8)
			{

				amp[inn][i]*=float(pow(-1,i));

			}



			td[inn][i]=td[inn][i+1]+alpha*onep5_period[i];



			while( (onep5_period[i])+td[inn][i] >= dur )
			{

				td[inn][i]=float((rand()/max)*(dur - onep5_period[i]));



			}

		}



//	printf(" inn=%ld\n",inn);



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


			nhs[inn][i]= nlimit;

			while(1)
			{

				if( (float(nhs[inn][i])/(2.*f[i]))+td[inn][i] >= dur )
				{

					nhs[inn][i]-=2;

				}

				else
				{

					break;
				}

			}



//		printf("%ld  %12.4g  %12.4g  %ld  %12.4g \n",i,f[i],amp[inn][i],nhs[inn][i],td[inn][i]);

		}



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

			if( float(nhs[inn][i]) < 3 )
			{

				printf("s3: error: i=%ld  nhs= %9.4g  td=%9.4g   dur=%9.4g \n",i,nhs[inn][i],td[inn][i],dur);

				exit(1);

				nflag=0;

			}

		}

	}

}

void print_th(void)
{

	char ff[20],temp[20];



	strcpy(ff,filename[3]);

	sprintf (temp,"_%ld", iwin );

	strcat (ff, temp );

	strcat (ff, ".out" );



	pFile[3]=fopen(ff, "w");





	printf("\n Acceleration file:  %s   - time(sec) & accel(%s) ",ff,aunit);





	for(j=0; j<nt; j++)

	{

		a=0;

		t=dt*j;





		for(i=0; i<nspec; i++)

		{

			tt=t-td[igen][i];



			if( t>=td[igen][i] && tt <= nhs[igen][i]/(2.*f[i])  )

			{

				a+=float(amp[igen][i]*sin(omegaf[i]*tt/float(nhs[igen][i]) )*sin(omegaf[i]*tt));

			}

		}

		fprintf(pFile[3],"%14.6f \t %12.4f \n",t,a);

		th[j]=a;

	}

	fclose(pFile[3]);

}

void srs()

{



/******************************** initialize arrays ******************/

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

   {

      xmax[i]=0.;

      xmin[i]=0.;

	  xxmin[i]=0.;

	  xxmax[i]=0.;



	  x[i]=0.;

	  xb[i]=0.;

      xbb[i]=0.;

	  yb[i]=0.;

	  ybb[i]=0.;

   }



/******************************** initialize constants ***************/



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



/*********************** initialize filter coefficients ***************/





//   printf(" nspec= %ld \n",nspec);





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



		i=0;

		while(1)

		{

			y=th[i];





			for(j=0; j<nspec; ++j)

			{

				CORE

			}

			i++;



			if(i>=nt){break;}



		}



		for(j=0; j<nspec; ++j)

		{

			if( fabs(xmax[j]) <= 1.0e-20 || fabs(xmin[j]) <= 1.0e-20 )

			{

				printf(" Error SRS  j=%ld  xmax=%8.4g  xmin=%8.4g \n",j,xmax[j],xmin[j]);

				printf("            f[j]=%8.4g  \n",f[j]);



				exit(1);

			}

        }



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



   float tlast=t;

   long ijk;





   y=0.;



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

   {

	    t=tlast+dt*ijk;



		for(j=0; j<=nspec; ++j)

		{

			if( dt*ijk > over_period[j] ){break;}



			CORE

		}

   }



   for(j=0; j<=nspec; ++j)

   {

        xxmax[j]=xmax[j];

   }



   for(j=0; j<=nspec; ++j)

   {

		xmin[j]=float(fabs(xmin[j]));



		if(xxmax[j]<float(fabs(xmin[j])))

		{

			xxmax[j]=float(fabs(xmin[j]));

		}

   }



   for(j=0; j<=nspec; ++j)

   {

		xxmin[j]=float(fabs(xmin[j]));



		if(xmax[j]<float(fabs(xxmin[j])))

		{

			xxmin[j]=float(fabs(xmax[j]));

		}

//		printf(" j=%ld xmax=%lf xmin=%lf xxmin=%lf\n",j,xmax[j],xmin[j],xxmin[j]);

   }



}

void scale(void)

{

	float ave,ss;



	for(i=0; i<nspec; i++)

	{

        ave = float(( fabs(xmax[i]) + fabs(xmin[i]) )/2.);



		ss= float(pow( ( spec[i]/ave ), exponent));



		if(ss>=1.0e-20 && ss <=1.0e+20)

		{

		}

		else

		{

			printf("\n scale error: i=%ld  ss=%8.4g  ave=%8.4g  spec=%8.4g\n",i,ss,ave,spec[i]);



			exit(1);

		}



		amp[inn][i]*=ss;



//		printf("%ld amp=%10.3g spec=%10.3g xxmin=%10.3g ss=%10.3g\n",i,amp[inn][i],spec[i],xxmin[i],ss);



//		getch();

	}





}

void gen_time(void)
{

		float ta,rms;

//		printf(" ref1 \n");


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

			if( float(nhs[igen][i]) < 3 )
			{

				printf("gt error: i=%ld  nhs= %9.4g \n",i,nhs[igen][i]);


				exit(1);

			}

		}

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

			if(  omegaf[i] <= MINF || omegaf[i] > 1.0e+20 )
			{

				printf("gt error: i=%ld  omegaf= %9.4g \n",i,omegaf[i]);

				exit(1);

			}

		}

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

			if(  fabs(amp[igen][i])<1.0e-20 )

			{

				printf("gt error: i=%ld  amp= %9.4g \n",i,amp[igen][i]);



				exit(1);

			}

		}

		for(i=0; i<nspec; i++)

		{

			if(  fabs(upper[i])<1.0e-20 )

			{

				printf("error: i=%ld  upper= %9.4g \n",i,upper[i]);



				exit(1);

			}

		}







		if(nv==0)

		{

			for(j=0; j<nt; j++)

			{



				t=dt*j;



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



				        if(  omegaf[i] <= MINF)
						{

							printf("ref 1: omegaf error \n");

							exit(1);

						}


					wavelet[i][j]=0.;



					tt=t-td[igen][i];


					ta= omegaf[i]*tt;



//					printf(" i=%ld  j=%ld  t=%8.4g  \n",i,j,t);



					if( t>=td[igen][i] && tt <= upper[i] )

					{

						wavelet[i][j]=float(sin(ta/float(nhs[igen][i]) )*sin(ta));



//						printf(" w=%9.4g  nhs=%ld  ta=%8.4g  omegaf[i]=%8.4g\n",wavelet[i][j],nhs[igen][i],ta,omegaf[i]);

						if(  omegaf[i] <= MINF)
						{

							printf("ref 2: omegaf error \n");

							exit(1);

						}



						if( fabs(wavelet[i][j]) > 1.0e+10)
						{

							printf(" wavelet error: i=%ld  j=%ld %9.4g\n",i,j,wavelet[i][j]);

							exit(1);

						}
					}
				}

			}

		}



//		printf(" ref2 \n");





		rms=0.;



		for(j=0; j<nt; j++)

		{

				a=0;



				for(i=0; i<nspec; i++)

				{

						a+=amp[igen][i]*wavelet[i][j];

				}



				th[j]=a;



				rms+=float(pow(a,2.));

		}



		rms=float(sqrt(rms/double(nt)));







//		printf(" nt=%ld nspec=%ld  rms = %12.4g \n",nt,nspec,rms);



		if(rms<1.0e-20)

		{

			printf("\n rms error \n");

			exit(1);

		}



//		printf(" ref3 \n");



}

void max_param(void)

{





	ymax=0.;

    vmax=0.;

	dmax=0.;



	area=0.;



	for(j=0; j<nt; j++)

	{

		vth[j]=0.;

	}





	for(j=0; j<nt; j++)

	{

		if( float(fabs(th[j])>ymax)){ymax=float(fabs(th[j]));}

	}



	for(i=0; i<=nt; ++i)

	{

      if(i==0 || i==nt)

      {

	      area+=th[i]*dth;

      }

      else

      {

	      area+=th[i]*dt;

      }

      if( float(fabs(area)>vmax)){vmax=float(fabs(area));}



      vth[i]=area;



	}



	area=0.;



	for(i=0; i<=nt; ++i)

	{

      if(i==0 || i==nt)

      {

	      area+=vth[i]*dth;

      }

      else

      {

	      area+=vth[i]*dt;

      }

      if( float(fabs(area))>dmax){dmax=float(fabs(area));}



	}



	if(iunit==1)

	{

		vmax*=float(386.);

		dmax*=float(386.);

	}

	if(iunit==2)

	{

		vmax*=float(9.81);

		dmax*=float(9810.);

	}

	if(iunit==3)

	{

		vmax*=float(1.);

		dmax*=float(1000.);

	}



}

void srs_error(void)

{



	float emax,emin;



	error=0.;

	irror=0.;





    for(i=0; i<nspec; i++)

	{



	  if(spec[i] < 1.0e-20)

	  {

		  printf("\n spec error. %ld %12.6e \n",i,spec[i]);

	  }



      emax=float(fabs( log10(xmax[i]/spec[i]) ));

      emin=float(fabs( log10(xmin[i]/spec[i]) ));



      error+=emax;

      error+=emin;



	  if( emax > irror){irror= emax;}

	  if( emin > irror){irror= emin;}

	}





}

void moreinput(void)

{





	printf("\n Enter damping format for SRS. \n 1= damping ratio   2= Q \n");

	scanf("%ld",&idamp);



	  if(idamp == 1 )

	  {

	     printf("\n Enter SRS damping ratio (typically 0.05) ");

         scanf("%f",&damp);

	  }

	  else

	  {

	     printf("\n Enter SRS amplification factor Q (typically 10) ");

         scanf("%f",&damp);

	     damp = float(1./(2.*damp));

	  }



	  if(damp >= 1.0)

	  {

		  printf("\n Error:  damping value must be < 1 \n");

		  exit(1);

	  }





	  printf("\n\n Enter the number of trials. \n");

	  scanf("%ld",&ntrials);



	  if(ntrials> MAXTRIALS)

	  {

		  ntrials=MAXTRIALS;



          printf("\n Warning: number of trials reduced to %ld \n",ntrials);

	  }





	  while(1)

	  {

		printf("\n\n Enter units ");

		printf("\n  1=English:  G,       in/sec, in   ");

		printf("\n  2=metric:   G,       m/sec,  mm   ");

		printf("\n  3=metric:   m/sec^2, m/sec,  mm  \n");



		scanf("%ld",&iunit);



        if( iunit==1 || iunit==2 || iunit==3 ){break;}



	  }





	  strcpy(dunit,"mm");

	  strcpy(vunit,"m/sec");

	  strcpy(aunit,"G");



	  if(iunit==1)

	  {

		  strcpy(dunit,"inch");

		  strcpy(vunit,"in/sec");

	  }

	  if(iunit==3)

	  {

		  strcpy(aunit,"m/sec^2");

	  }





//	  printf("\n\n Enter strategy \n  1=Random  2=Reverse Sine Sweep  3=Combination\n");

//	  scanf("%ld",&ichoice);





	  printf("\n\n Enter arbitrary integer (1 to 1000) \n");

	  scanf("%ld",&iseed);



	  for(i=1; i<iseed; i++)

	  {

		a=float(rand()/max);

	  }



}

void rank(void)

{



    long itemp;



	float temp;

/*

	for(i=0; i<ntrials; i++)

	{

	    printf("ref 2  i=%ld dmi=%12.6f \n",i,dm[i]);

	}

*/

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

		yrank[i]=i;

		vrank[i]=i;

		drank[i]=i;

		erank[i]=i;

		irank[i]=i;



		pyrank[i]=0;

		pvrank[i]=0;

		pdrank[i]=0;

		perank[i]=0;

		pirank[i]=0;

	}



	for(i=0; i< (rntrials-1); i++)
	{

		for(j=i+1; j<rntrials; j++)
		{

 //           fprintf(pFile[6],"i=%ld j=%ld ymi=%12.3f ymj=%12.3f\n",i,j,ym[i],ym[j]);

//			        printf("i=%ld j=%ld dmi=%12.6f dmj=%12.6f \n",i,j,dm[i],dm[j]);



			if(ym[i]<ym[j])

			{



				temp=ym[i];

                ym[i]=ym[j];

				ym[j]=temp;



				itemp=yrank[i];

                yrank[i]=yrank[j];

				yrank[j]=itemp;

			}



 //           printf("i=%ld j=%ld vmi=%12.4f vmj=%12.4f \n",i,j,vm[i],vm[j]);

//			        printf("i=%ld j=%ld dmi=%12.6f dmj=%12.6f \n",i,j,dm[i],dm[j]);



			if(vm[i]<vm[j])

			{

				temp=vm[i];

                vm[i]=vm[j];

				vm[j]=temp;



				itemp=vrank[i];

                vrank[i]=vrank[j];

				vrank[j]=itemp;

			}



//            printf("i=%ld j=%ld dmi=%12.6f dmj=%12.6f \n",i,j,dm[i],dm[j]);



			if(dm[i]<dm[j])

			{

				temp=dm[i];

                dm[i]=dm[j];

				dm[j]=temp;



			    itemp=drank[i];

                drank[i]=drank[j];

				drank[j]=itemp;

			}





			if(em[i]<em[j])

			{

				temp=em[i];

                em[i]=em[j];

				em[j]=temp;



			    itemp=erank[i];

                erank[i]=erank[j];

				erank[j]=itemp;

			}



			if(im[i]<im[j])

			{

				temp=im[i];

                im[i]=im[j];

				im[j]=temp;



			    itemp=irank[i];

                irank[i]=irank[j];

				irank[j]=itemp;

			}



		}

    }



	pFile[6]=fopen(filename[6], "w");



	fprintf(pFile[6],"\n\n accel \n");

	for(i=0; i<rntrials; i++)

	{

		fprintf(pFile[6],"%ld %12.3f \n",yrank[i],ym[i]);

	}



	fprintf(pFile[6],"\n\n vel \n");

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

		fprintf(pFile[6],"%ld %12.3f \n",vrank[i],vm[i]);

	}



	fprintf(pFile[6],"\n\n disp \n");

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

		fprintf(pFile[6],"%ld %14.4g \n",drank[i],dm[i]);

	}



	fprintf(pFile[6],"\n\n error \n");

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

		fprintf(pFile[6],"%ld %12.3f \n",erank[i],em[i]);

	}

	fprintf(pFile[6],"\n\n irror \n");

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

		fprintf(pFile[6],"%ld %12.3f \n",irank[i],im[i]);

	}



	fclose(pFile[6]);



	long nmin=10000;





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

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

            if(yrank[i]==j)
			{

				pyrank[j]=i;


                break;

            }
        }
	}


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

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

            if(vrank[i]==j)
			{

				pvrank[j]=i;



                break;

            }

        }

	}



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

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

            if(drank[i]==j)
			{

				pdrank[j]=i;



                break;

            }

        }

	}



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

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

            if(erank[i]==j)
			{

				perank[j]=i;



                break;

            }

        }

	}



	for(i=0; i<rntrials; i++)

	{

		for(j=0; j<rntrials; j++)

		{

            if(irank[i]==j)

			{

				pirank[j]=i;



                break;

            }

        }

	}

}

void weights()

{



	pFile[5]=fopen(filename[5], "w");



    float nmax=0.;





	for(i=0; i<rntrials; i++)

	{



//  weights



//		nrank[i]=(1*pirank[i]+1*perank[i])+(2*pdrank[i]+1*pvrank[i]+1*pyrank[i]);





		nrank[i]=float((iw*pirank[i]+ew*perank[i])+(dw*pdrank[i]+vw*pvrank[i]+aw*pyrank[i]));





		if( nrank[i]>nmax)

		{

			nmax=nrank[i];

			iwin=i;

		}

		fprintf(pFile[5],"%ld \t %lf \t %ld \t %ld \t %ld \t %ld \t %ld\n",i,nrank[i],pyrank[i],pvrank[i],pdrank[i],perank[i],pirank[i]);

	}

	fclose(pFile[5]);





	char ff[20],temp[20];



	strcpy(ff,filename[2]);

	sprintf (temp,"_%ld", iwin );

	strcat (ff, temp );

	strcat (ff, ".out" );



	pFile[2]=fopen(ff, "w");



	strcpy(ff_wt,ff);









	fprintf(pFile[2]," No.   Freq(Hz)   Polarity    NHS     Delay(msec)    Amp(%s)\n",aunit);







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

		if(amp[iwin][i] >=0 )
		{

		   fprintf(pFile[2],"%ld \t %8.1f \t + \t %ld \t %8.2f \t %8.2f \n",i,f[i],nhs[iwin][i],1000.*td[iwin][i],fabs(amp[iwin][i]));

		}
		else
		{

		   fprintf(pFile[2],"%ld \t %8.1f \t - \t %ld \t %8.2f \t %8.2f \n",i,f[i],nhs[iwin][i],1000.*td[iwin][i],fabs(amp[iwin][i]));

		}

	}



	fclose(pFile[2]);



}

void srs_pre(void)

{



	printf("\n");



	for(j=0;j<nspec;++j)

	{





		if(algorithm==1)

		{



			omega=float(tpi*f[j]);

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

			cosd=float(cos(omegad*dt));

			sind=float(sin(omegad*dt));

			domegadt=float(damp*omega*dt);



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

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

			b1[j]=float(2.*domegadt);

			b2[j]=float(omega*dt*exp(-domegadt));

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

			b3[j]=0.;



		}

		else

		{

			double C,E,K,S,Sp;





			omega=float(tpi*f[j]);

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



			E=float(exp(-damp*omega*dt));

			K=float(omegad*dt);

			C=float(E*cos(K));

			S=float(E*sin(K));



			Sp=S/K;



			a1[j]=float(2*C);

			a2[j]=float(-pow(E,2));



			b1[j]=float(1.-Sp);

			b2[j]=float(2.*(Sp-C));

			b3[j]=float(pow(E,2)-Sp);



		}



//		printf("%ld om=%8.4g a1=%8.3g a2=%8.3g b1=%8.3g b2=%8.3g b3=%8.3g\n",j,omega,a1[j],a2[j],b1[j],b2[j],b3[j]);

   }



}

void zero(void)

{

	for(i=0; i<MAXTRIALS; i++)

	{

	     ym[i]=0.;

		 vm[i]=0.;

		 dm[i]=0.;

		 em[i]=0.;

    }

	for(i=0; i< K100; i++)

	{

		 th[i]=0.;

		vth[i]=0.;

	}

}

void print_srs(void) // 555

{

	char ff[20],temp[20];



	strcpy(ff,filename[4]);

	sprintf (temp,"_%ld", iwin );

	strcat (ff, temp );

	strcat (ff, ".out" );



	pFile[4]=fopen(ff, "w");



	printf("          SRS file:  %s  - fn(Hz) pos(%s) neg(%s) \n",ff,aunit,aunit);

    printf("     Wavelet Table:  %s \n",ff_wt);



	for(i=0; i<nspec; i++)

	{

				fprintf(pFile[4],"%12.3f \t %12.5e \t %12.5e\n",f[i],xmax[i],xmin[i]);

	}

	fclose(pFile[4]);



}

void rankfunctions(void)
{

			rank();   // 333



		    iw=3.0;

			ew=0.5;

			dw=5.;

			vw=1.;

			aw=1.;



			strcpy(filename[5],"disp_rank");

			weights();

		    igen=iwin;





			printf("\n\n Optimum case for displacement = %ld \n",iwin);

		    printf("\n   Peak Accel = %12.3f %s    ",ym[pyrank[iwin]],aunit);

            printf("\n   Peak Velox = %12.3f %s    ",vm[pvrank[iwin]],vunit);

		    printf("\n   Peak Disp  = %12.3f %s    ",dm[pdrank[iwin]],dunit);

	        printf("\n   Max Error  = %12.3f dB \n\n",20.*im[pirank[iwin]]);



			gen_time();

			print_th();

			integrate();

			srs();

			print_srs();





		    iw=3.0;

			ew=0.5;

			dw=1.;

			vw=1.;

			aw=5.;



			strcpy(filename[5],"accel_rank");

		    weights();

			igen=iwin;



			printf("\n\n Optimum case for acceleration = %ld \n",iwin);

		    printf("\n   Peak Accel = %12.3f %s    ",ym[pyrank[iwin]],aunit);

            printf("\n   Peak Velox = %12.3f %s    ",vm[pvrank[iwin]],vunit);

		    printf("\n   Peak Disp  = %12.3f %s    ",dm[pdrank[iwin]],dunit);

	        printf("\n   Max Error  = %12.3f dB \n\n",20.*im[pirank[iwin]]);





			gen_time();

			print_th();

			integrate();

			srs();

			print_srs();



}

void integrate(void)

{

	float lambda;

	float apb[MAX],amb[MAX];



	float alpha_hat[MAX];

	float beta_hat[MAX];



	float a,v,d;



	float da,va;



	float ta;



	char ff[20],temp[20];



	strcpy(ff,filename[7]);

	sprintf (temp,"_%ld", iwin );

	strcat (ff, temp );

	strcat (ff, ".out" );



	pFile[7]=fopen(ff, "w");



	if(pFile[7]==NULL)

	{

		printf("\n Error: %s \n",ff);

		exit(1);

	}

	printf("\n     Velocity file:  %s   - time(sec) & vel(%s) ",ff,vunit);







	strcpy(ff,filename[8]);

	sprintf (temp,"_%ld", iwin );

	strcat (ff, temp );

	strcat (ff, ".out" );



	pFile[8]=fopen(ff, "w");




	if(pFile[8]==NULL)
	{

		printf("\n Error: %s \n",ff);

		exit(1);

	}

	printf("\n Displacement file:  %s   - time(sec) & disp(%s) \n",ff,dunit);





	if(iunit==1)
	{

		lambda=float(386.);

	}
	if(iunit==2)
	{

		lambda=float(9.81);

	}
	if(iunit==3)
	{

		lambda=1.;

	}





	float upper[MAX];



	for(i=0; i<nspec; i++)

	{

		upper[i] = float(nhs[igen][i]/(2.*f[i]));

	}





	for(i=0; i<nspec; i++)

	{

		alpha_hat[i] = omegaf[i]/float(nhs[iwin][i]);

		beta_hat[i]  = omegaf[i];



		apb[i]=alpha_hat[i]+beta_hat[i];

		amb[i]=alpha_hat[i]-beta_hat[i];



		if(fabs(amb[i]) <= 1.0e-20)

		{

		    printf("\n warning:  i=%ld  nhs=%ld  alpha=%10.4g  beta=%10.4g \n",i,(nhs[iwin][i]),alpha_hat[i],beta_hat[i]);



		    exit(1);

		}



	}





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

		a=0;

		v=0.;

		d=0.;

		t=dt*j;



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



			tt=t-td[igen][i];



			ta= omegaf[i]*tt;



			if( t>=td[igen][i] && tt <= upper[i] )

			{

				a+=float(amp[iwin][i]*sin(alpha_hat[i]*tt )*sin(beta_hat[i]*tt));



				va=float(-(1./apb[i])*sin(apb[i]*tt)+(1./amb[i])*sin(amb[i]*tt));



				v+=float(va*amp[iwin][i]);



				da=  float((1./pow(apb[i],2.))*(cos(apb[i]*tt)-1));

				da-= float((1./pow(amb[i],2.))*(cos(amb[i]*tt)-1));



				d+= float(da*amp[iwin][i]);

			}

		}

		fprintf(pFile[7],"%14.6e \t %14.6e \n",t,v*lambda/2.);

		fprintf(pFile[8],"%14.6e \t %14.6e \n",t,d*lambda/2.);



	}

	fclose(pFile[7]);

	fclose(pFile[8]);



}

void elapsed()
{



	time( &ltimen );



	double tnow=double(ltimen);



//    printf("\n %10.4g %10.4g ", ltimen,tstart);



    long irem=long(tnow-tstart);



	if(  (tnow-tstart) > 0.)
	{

  	  long ihr= long( (irem                      ) / 3600.);

      long imin=long( (irem - ihr*3600           ) / 60.  );

      long isec=long( (irem - ihr*3600 - imin*60 ) / 1.   );



//	  printf("\n k= %ld\n",*k);



      printf("\n\n Elapsed time (hr:min:sec)\n");



      printf("      = %ld:%ld:%ld \n", ihr,imin,isec );

	}

}