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


//************ complex number handling ********************
typedef struct 
{
    float r,i;
	
}floatcomplex;



#define  COMPLEXMAG(x,y,z) z=sqrt( pow( x,2) + pow(y,2) );
#define COMPLEXMAG2(x,y,z) z=    ( pow( x,2) + pow(y,2) );

#define COMPLEXRECIP(ar, ai, xr, xi)              \
{                                                 \
           xr=  ar/( pow( ar,2) + pow(ai,2) );    \
           xi= -ai/( pow( ar,2) + pow(ai,2) );    \
}                                                 \

#define COMPLEXMUL(ar, ai, br, bi,cr,ci)          \
{                                                 \
           cr=  ar*br -ai*bi;                     \
           ci=  ar*bi +ai*br;                     \
}                                                 \

#define COMPLEXDIV(nr, ni, dr, di,cr,ci)          \
{                                                 \
           cr =  nr*dr +ni*di;                    \
           cr/= ( pow( dr,2) + pow(di,2) );       \
           ci =  ni*dr -nr*di;                    \
           ci/= ( pow( dr,2) + pow(di,2) );       \
}                                                 \

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

//   262144 
//   131072 
//    65536 
//    32768
//    16384

#define FFTMAX   65536  
#define FFTMAXP1 32768

#define NUM 3000

void inputpsd(void);

void white(void);
void gauss(void);
void ifft(void);
void scalefft(void);
void fft(void);

void intlog(void);

void slope_calc(void);
void slope_calc_v(void);
void slope_calc_d(void);

void bitrev(void);
void ctfourier(void);

void forcesym(void);
void ifft_bitrev(void);
void ifft_ctfourier(void); 

void scaleth(void);
void area(void);

void start(void);

void print_output(void);

void header(void);


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

        float velox[FFTMAX];
		float disp[FFTMAX];

		float vsum=0.;
		float vave=0.;

		float dsum=0.;
		float dave=0.;


		float amax=0.;
		float amin=0.;

		float vmax=0.;
		float vmin=0.;

		float dmax=0.;
		float dmin=0.;

long i,ich,ig,j,jj,jjj,k;
long kt=0; 
long m,n,nf,nn,nout,np,npsd,nstages,nsegment,ns;

long iseed;

const float pi=float(atan2(0.,-1.));
const float tpi= float(2.*pi);
const float e= float(1./sqrt(2.*pi));
const float max= float(32767.);

floatcomplex y[FFTMAX];
floatcomplex ww[FFTMAX];



float a[FFTMAX],ae,aq,ave;
float b[NUM];
float d,dam,delta,df;

float f,ff,fft_freq[100000],freq[NUM],fmax,f1,f2;
float grms,grmsout;

float mag,ms;

float psd[NUM];
float vpsd[NUM];
float dpsd[NUM];

float ratio;

float sdg;
float slope[NUM],slopev[NUM],sloped[NUM];
float ss,std,stddev,sum[10000];
float t,theta,tt;
float x;
float yt;
float var;
float z;

long iflag=0;


double dt;
double sr;
double dur;
double tmax,tmax_old;
double ttt,tlast;
double vrms,drms;


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


void main()
{
		vsum=0.;
		dsum=0.;


	    printf("\n\n psdgen.cpp    ver 3.5   January 15, 2009");
		printf("\n\n By Tom Irvine,  Email: tomirvine@aol.com \n");


	    printf("\n\n This program synthesizes a time history to satisfy");
		printf("\n a specified power spectral density function. ");

		printf("\n\n The input file must have two columns: ");
		printf("\n    Frequency(Hz) & PSD(G^2/Hz) ");

		printf("\n\n The input file may contain up to 100 coordinates. ");

		
	    inputpsd();
		intlog();

		long nop = FFTMAX;

//		printf("\n The program is limited to %ld output points. \n",nop);

		printf("\n Enter the duration (seconds) \n");
		scanf("%lf",&tmax);

		printf("\n Enter the sample rate (samples/sec) \n");
		scanf("%lf",&sr);

	    tmax_old = tmax;

		if(tmax*sr > FFTMAX  )
		{
//			printf("\n **Warning: too many samples. \n Lowering duration. \n");
			tmax= float(FFTMAX)/ sr;
		}

        if( sr< 2.*fmax)
		{
			printf("\n Error:  the sample rate must be >= %14.5g \n",2.*fmax);
			exit(1);
		}

		
		printf("\n\n Enter arbitrary integer (1 to 1000) \n");
		scanf("%ld",&iseed);


		nout=long(tmax*sr);

		dur = float(2.1*tmax);
		
        nf=long(fmax/df);

		dt= float(1./sr);

		long kjk;


		start();


        for(kjk=0; kjk < 400; kjk++)
        {
			if(iflag==1){break;}
	 
			printf("\n white ");
			white();
			printf("\n fft ");
			fft();
			printf("\n scalefft ");
			scalefft();		
			printf("\n ifft ");
			ifft();
			printf("\n scaleth ");
			scaleth();	

		    print_output();

		}
		

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

		printf("\n\n max accel = %10.3g G",amax);
		printf("\n min accel = %10.3g G \n",amin);

		
		printf("\n max vel = %10.3g in/sec",vmax);
		printf("\n min vel = %10.3g in/sec \n",vmin);

		
		printf("\n max disp = %10.3g in",dmax);
		printf("\n min disp = %10.3g in \n",dmin);

		printf("\n Output files \n");
		printf("\n %s \n %s \n %s \n",filename[2],filename[3],filename[4]);


//	   printf("\n\n Please call the output files into your own plotting program.");
       printf("\n Output file units: accel(G), vel(in/sec), disp(in) \n");

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


}
void slope_calc()
{
	//  calculate slopes

    
	      for( i=0; i < npsd-1; i++)
		  {
			  if(freq[i] > 1.0e-12)
			  {

				  if( freq[i+1]/freq[i] <= 0. )
				  {
					  printf("\n frequency error: %12.4g %12.4g \n",freq[i+1],freq[i]);
					  break;
				  }
				  if( psd[i+1]/psd[i] <= 0. )
				  {
					  printf("\n amplitude error: %12.4g %12.4g \n",psd[i+1],psd[i]);
					  break;
				  }

					slope[i]=float( log( psd[i+1]/psd[i] )/log( freq[i+1]/freq[i] ) );
			  }
			  else
			  {
				  slope[i] = 0.;
			  }
//	          fprintf(pFile[3],"\n freq[i] = %14.7e  slope = %14.7e\n",freq[i],slope[i]);
		  }

		float ra=0.;

		grms=0.;

		for( i=0; i < npsd-1; i++)
		{
			if(slope[i] < -1.0001 ||  slope[i] > -0.9999 )
			{
				ra+= float( ( psd[i+1]*freq[i+1]- psd[i]*freq[i])/( slope[i]+1.) );
			}
			else
			{
				if(freq[i] < 1.0e-12)
				{ 
					freq[i]= float(1.0e-12);
				}

				ra+= float( psd[i]*freq[i]*log( freq[i+1]/freq[i]) );
			}
		}
		grms=float( sqrt(ra) );



}

void inputpsd()
{
	
		
		printf("\n Enter the name of the input file. \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]);
			exit(1);
		}
		else
		{
			printf(" File: %s opened. \n", filename[i]);
		}

        npsd=0;   

		printf("\n\n");

        double fff,aaa;

		long k=0;

		header();

		for(i=0; i<NUM; i++)
		{
			if(fscanf(pFile[0],"%lf %lf",&fff,&aaa) <=0 ){break;}
	      
			if(fff < 0.)
			{
				printf("\n Error: frequencies must be >0;");
				exit(1);
			}
			if(aaa < 0.)
			{
				printf("\n Error: psd must be > 0;");
				exit(1);
			}
			if(aaa == 0.)
			{
				aaa = 1.0e-30;
			}
            if( fff > 0 )
			{
                freq[k]=fff;
				psd[k]=aaa;
				fmax=freq[k];

				printf("%14.5g  %14.5g \n",freq[k],psd[k]);


				if(k>1 && freq[k] < freq[k-1])
				{

					printf("\n k=%ld  freq[k]=%14.5g  freq[k-1]=%14.5g \n",k,freq[k],freq[k-1]);

					printf("\n Error: frequencies must be in ascending order");
					exit(1);
				}

				k++;
				npsd++; 
            }

		}
		if(npsd<=1)
		{
			printf("\n Error:  at least two coordinates are required.");
			exit(1);
		}

		printf("\n npsd=%ld \n",npsd);


        printf("\n ******* Output Filenames ******** \n");

		printf("\n\n Enter the acceleration time history filename. \n");
	    scanf("%s",&filename[2]);
	    pFile[2]=fopen(filename[2],"w");


		printf("\n\n Enter the velocity time history filename. \n");
	    scanf("%s",&filename[3]);
	    pFile[3]=fopen(filename[3],"w");


		printf("\n\n Enter the displacement time history filename. \n");
	    scanf("%s",&filename[4]);
	    pFile[4]=fopen(filename[4],"w");

		strcpy(filename[5],"dv.out");
	    pFile[5]=fopen(filename[5],"w");

}
void white(void)
{
       sdg=1.;

       np=ns;

	   if(np > FFTMAX){np=FFTMAX;}

       printf("\n white 1 np= %ld",np);

	   if( np > 2 )
	   {
//  The gauss function determines the probability for 
//  each Nsigma value.


          printf("\n gauss ");
          gauss();

	      ave = 0.;
	      ms = 0.;

		  long kj;
		  
		  if(iseed < 1){iseed=1;}

	      for(i=0; i<np; i++)
	      {  
			 for(kj=0;kj<=iseed;kj++)
			 {
				x=rand()/max;
			 }

             a[i]=0.;
		            
		     for(j=1; j<= 599; j++)
		     {

//   x is a probability value.  0 < x < 1 
			  
                if( x >= sum[j] &&  x <= sum[j+1] )
			    {

//	 a = Nsigma value which corresponds to probability x.

                    a[i]=b[j];
                    break;
                }
             }
//			 printf(" %14.7e \n",a[i]);
          

// 	  50% of the amplitude points are multiplied by -1.

		     yt=rand()/max;

		     if(yt <= 0.5){a[i]*=-1;}

		     ave+=a[i];
		     ms+=float(pow(a[i],2.));
	      } 
		  printf("\n white 2  np= %ld",np);
	   
	      ave/=float(np);
	      ms/=float(np);

	      var=ms-pow(ave,2.);
	      stddev=sqrt(fabs(var));

//	print the output file

//		  printf("\n sratio = %14.7e\n",(sdg/stddev) );
	   
//	      printf("\n\n Enter the name of the white noise time history output file. \n");
//	      scanf("%s",&filename[1]);

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

	      for(i=0; i<np; i++)
	      {
			 t=dt*i;
		     a[i]=(a[i]-ave)*(sdg/stddev);
//	         fprintf(pFile[1],"%14.7e %14.7e \n",t,a[i]);

//			  printf("%14.7e %14.7e \n",t,a[i]);
          }

//	      printf("\n\n The output file is:  a.out ");
//	      printf("\n\n This file has two columns:  time(sec) and amplitude.");
//	      printf("\n Please call this file into your own plotting program.\n\n");
	   }
	   fclose(pFile[1]);

}
      
void gauss(void)
{
       strcpy(filename[1],"gauss.out");
	   pFile[1]=fopen(filename[1],"w");

/*
	b[j] is a standard deviation.  The maximum b[j] is eight sigma.
     	
	sum[j] is the area under the Gaussian curve from 0 to b[j]. 
*/
	   float ab;
   	       
       for( j=1;j<=800;j++)
	   {
            b[j]=j/100.;
            ab=-b[j];
            ss=0.;
            delta=(b[j]-ab)/1000.;

            for( i=0; i<=1000; i++)
			{ 
               z=(delta*i)+ab;
               ss+=float( exp( ( -pow(z,2.) )/2.) );
            }

            sum[j]=ss*e*delta;
            aq=float( 1.-sum[j] );

//            fprintf(pFile[1]," %14.7e %14.7e %14.7e \n",b[j],sum[j],aq);

       }
//       fclose(pFile[1]);            
}		
		
//******************************************************


void fft(void)
{

    for(i=0; i< FFTMAX; i++)
    {
		 y[i].r=float(0.);
	     y[i].i=float(0.); 
	}
	for(i=0; i< FFTMAXP1; i++)
    {
		 ww[i].r=float(0.);
	     ww[i].i=float(0.); 
	}

	for(i=0; i<np; i++)
	{
		 y[i].r = a[i];

//		 printf(" %14.7e \n", a[i]);
	}


//****** Perform bit reversal *******

	   printf("\n bitrev ");
	   bitrev();

//******* Cooley-Tukey Method ******
	   printf("\n begin ctm");

	   for(i=0; i< FFTMAXP1; i++)
       {
		 ww[i].r=float(0.);
	     ww[i].i=float(0.); 
	   }
       printf("\n ctfourier ");
       ctfourier();  


}
void ctfourier(void)
{
	long i1;
	
	long iy=1;

    long l;
	long lg;
	long lk;
	long lwd2;

	long nworld;
	long wd2;
	long world=4;

	float arg;
    float am;

	floatcomplex asum;
	floatcomplex a;
	floatcomplex b;
	floatcomplex c;

    m=nsegment;

	printf("\n m= %ld  nstages= %ld \n\n",m, nstages );

    for(j=0; j<= m-2; j+=2)
    {
	    jj=j+1;
		a.r=  y[j].r;
	    a.i=  y[j].i;
		b.r= y[jj].r;
		b.i= y[jj].i;
		y[j].r = a.r + b.r;
		y[j].i = a.i + b.i;

		y[jj].r = a.r - b.r;
		y[jj].i = a.i - b.i;
	}
	world=4;

//     printf("\n begin world loop ");
	
    for(i=2; i<=nstages; i++)
	{ 

       wd2=world/2;

	   nworld=int(m/world);    
	  
//     printf("\n i=%ld / %ld ",i,nstages );
//     printf("\n  (wd2-1)=%ld ", (wd2-1) );

       for(lk=0; lk<= (wd2-1); lk++)
	   {
	        arg=tpi*float(lk)/world;
		 
		    ww[lk].r =float(cos(arg));
		    ww[lk].i =float(-sin(arg));

	   }  
	   lwd2=0;

//       printf("\n  nworld=%ld ",nworld); 
	    
	    for( jjj=1; jjj<= nworld; jjj++)
		{  
		       i1=1+lwd2;
	      
			   for(lg=1; lg <= wd2; lg++)      
			   {       
		           l=i1;
		           lwd2=l+int( pow(2,iy) );
		     
				   a.r=y[l-1].r;
				   a.i=y[l-1].i;
		           b.r=y[lwd2-1].r;
		           b.i=y[lwd2-1].i;   
					 
		           COMPLEXMUL(ww[lg-1].r,ww[lg-1].i,b.r,b.i,c.r,c.i)

				   y[l-1].r  = a.r + c.r;
			       y[l-1].i  = a.i + c.i;

				   y[lwd2-1].r  = a.r - c.r;
			       y[lwd2-1].i  = a.i - c.i;
					 
				   i1++;
		      }
	     }
	     iy++;
	     world*=2;

//                printf("\n iy=%ld world=%ld", iy, world);
	   }

//*** am is the number of points in the segment
	
	   am=float(m);
    
       asum.r=float(0.);
	   asum.i=float(0.);

//	   printf("\n\n Enter the name of the white noise FFT output file. \n");
//	   scanf("%s",&filename[1]);

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

       for(ig=0; ig< m; ig++)      
       {
	         y[ig].r/=am;
			 y[ig].i/=am;

			 fft_freq[ig]=ig*df;

             ff=ig*df;

//			 fprintf( pFile[1],"%15.6e %15.6e %15.6e \n", ff, y[ig].r, y[ig].i);
			 
			 
//			 COMPLEXMAG(y[ig].r,y[ig].i,z)


	         if(ig<(m/2))  
			 {
				 asum.r += y[ig].r;
				 asum.i += y[ig].i; 
			 }
       }
	   fclose(pFile[1]);
	   
}

void bitrev(void)
{
	  float dd;
 
	  long ia, ib, ik, ilk;
      long kk;
	  long nm1;
	  long md2, mm1;
	  
	   m=nsegment;
       mm1=m-1;
       dd=float( nstages );
       nm1=long(dd);
       md2=m/2;

       printf("\n m=%ld nm1= %ld  md2=%ld  ns=%ld ", m, nm1, md2, ns);
	   
       for( i=1; i<= nm1; i++)   
       { 
		  for(ik=0; ik<m; ik++)
		  {
		    ww[ik].r=y[ik].r;
		  }
	  
	      kk=0;
	      ia=0;

	      ilk=long( pow(2,(i-1)) );

	      for(j=1; j<=ilk; j++)
		  {
			  ib=ia+md2-1;
	     
	         for(k=ia; k<=ib; k++)
			 {
		        if(k==ia)
				{
					kk=k;
				}
		
		        y[k].r = ww[kk].r;
		        y[k].i = float(0.);
		
				y[k+md2].r = ww[kk+1].r;
				y[k+md2].i =float(0.);
		
				kk=kk+2;
	         }
	         ia += md2*2;

	      }
	      md2/=2;   
	  }
 }
 void scalefft(void)
 {
	 float scale,spec,z;

//     printf("\n\n Enter the name of the scaled FFT output file. \n");
//	 scanf("%s",&filename[1]);

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


	 printf("\n ig = %ld \n",ig);
 

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

//		 z=sqrt( pow( y[i].r,2) + pow(y[i].i,2) );

		 z=0.;
		 z=float(  atan2(y[i].i,y[i].r) );

		 scale=float( 1.0e+16 );

		 for(j=0; j<= npsd; j++)
		 {
			 spec=0.;

			 if( fft_freq[i] >= freq[j] && fft_freq[i] < freq[j+1] )
			 {
				
				 spec=float( psd[j]*( pow( (fft_freq[i]/freq[j] ),slope[j]) ) );

				 scale = float( sqrt(spec)/z );
/*
				 if( fft_freq[i]>1000)
				 {
					 fprintf(pFile[3],"%ld %14.4e %14.4e %14.4e %14.4e %14.4e %14.4e %14.4e\n",j,fft_freq[i],freq[j],freq[j+1],slope[j],spec,scale,z);
				 }
*/
				 if( fabs(spec) > 1000000.)
				 {
					 printf("\n spec = %12.4g \n",spec);
					 exit(1);
					 break;
				 }

				 break;
			 }
		 }

         mag= float( sqrt(spec*df) );

		 y[i].r=float(mag*cos(z));
		 y[i].i=float(mag*sin(z));

		 float ymag = float( ( pow( y[i].r,2) + pow(y[i].i,2) ) );
/*
		 if(df*i > 1000. )
		 {
			fprintf(pFile[3],"%14.5e %14.4e %14.4e %14.4e %14.4e\n",df*i,spec,y[i].r,y[i].i,ymag);
         }
*/

	 }
  
//   fclose(pFile[1]);

 }
//**********************************************************
//**********************************************************

 /*

      ifft.cpp, ver 1.0, 11 Aug 1998

      by 
      Tom Irvine
   

*********************************************************************    

  This program calculates the inverse Fourier transform of an
  input time history.  A inverse Fast Fourier transform 
  implementation is used based on the Cooley-Tukey method.


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

//   262144 
//   131072 
//    65536 
//    32768
//    16384

void ifft()
{
	for(i=0; i< FFTMAXP1; i++)
    {
		 ww[i].r=float(0.);
         ww[i].i=float(0.); 
	}

//*** Fourier transform *****
//*** force symmetry ******

	   forcesym();

//****** Perform bit reversal *******

	   ifft_bitrev();

//******* Cooley-Tukey Method ******

	   for(i=0; i< FFTMAXP1; i++)
       {
		 ww[i].r=float(0.);
         ww[i].i=float(0.); 
 	   }
//       printf("\n call function ctfourier ");
       ifft_ctfourier();  
    
}
void ifft_ctfourier(void)
{
	long i1;
	
	long iy=1;
	
    long l;
	long lg;
	long lk;
	long lwd2;

	long nworld;
	long wd2;
	long world=4;

	float arg;

    float am;

	floatcomplex asum;
	floatcomplex a;
	floatcomplex b;
	floatcomplex c;

//	printf("\n  in ctfourier1 ");

    m=nsegment;

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

    for(j=0; j<= m-2; j+=2)
    {
        jj=j+1;
		a.r=  y[j].r;
        a.i=  y[j].i;
		b.r= y[jj].r;
		b.i= y[jj].i;
		 y[j].r = a.r + b.r;
		 y[j].i = a.i + b.i;

		y[jj].r = a.r - b.r;
		y[jj].i = a.i - b.i;
	}
	world=4;

//	printf("\n begin world loop ");
	
    for(i=2; i<=nstages; i++)
	{ 

       wd2=world/2;

	   nworld=int(m/world);	   
          
//	   printf("\n i=%ld / %ld ",i,nstages );
//       printf("\n  (wd2-1)=%ld ", (wd2-1) );

       for(lk=0; lk<= (wd2-1); lk++)
	   {
               arg=tpi*float(lk)/world;
                 
		    ww[lk].r =float(cos(arg));
        	ww[lk].i =float(+sin(arg));

	    }  
        lwd2=0;

 //       printf("\n  nworld=%ld ",nworld); 
        for( jjj=1; jjj<= nworld; jjj++)
		{  
		       i1=1+lwd2;
              
			   for(lg=1; lg <= wd2; lg++)	   
			   {	   
                     l=i1;
                     lwd2=l+int( pow(2,iy) );
                     
					 a.r=y[l-1].r;
					 a.i=y[l-1].i;
                     b.r=y[lwd2-1].r;
                     b.i=y[lwd2-1].i;   
					 
                       COMPLEXMUL(ww[lg-1].r,ww[lg-1].i,b.r,b.i,c.r,c.i)

				        y[l-1].r  = a.r + c.r;
                        y[l-1].i  = a.i + c.i;

					 y[lwd2-1].r  = a.r - c.r;
	              	 y[lwd2-1].i  = a.i - c.i;
					 
					 i1++;
                }
	      }
          iy++;
          world*=2;

//		  printf("\n iy=%ld world=%ld", iy, world);
	   }

//*** am is the number of points in the segment
	
	   am=float(m);
    
//	   printf("\n am=%f ",am);

       asum.r=float(0.);
	   asum.i=float(0.);

//	   printf("\n m=%d ",m);
       for(ig=0; ig< m; ig++)      
       {
//             y[ig].r/=am;
//			 y[ig].i/=am;

//			 freq=ig*df;

//			 fprintf( pFile[2],"%15.6e %15.6e %15.6e \n", freq, y[ig].r, y[ig].i);
//			 fprintf( pFile[3],"%15.6e %15.6e  \n", freq, y[ig].r);

             if(ig<(m/2))  
			 {
				 asum.r += y[ig].r;
				 asum.i += y[ig].i; 
			 }
       }

//     printf("\n\n sum real=%f   sum imag=%f",asum.r,asum.i);
//       printf("\n\n complex time history written to:  a.out \n");
//	   printf("  \n    real time history written to:  a.r \n");
//	   fclose( pFile[2] );
	   
}
void forcesym(void)
{
//      printf("\n\n ns = %ld ",ns);

      for(i=1; i < (ns/2); i++)
      {
              y[ns-i].r =  y[i].r;
              y[ns-i].i = -y[i].i;
      }


//	  printf("\n\n Enter name of forced symmetry file. \n");
//	  scanf("%s",filename[1]);

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


//	  for(i=0; i<ns; i++)
//	  {
//		  fprintf(pFile[1],"%14.7e %14.7e %14.7e \n",i*df,y[i].r,y[i].i);
//	  }
//	  fclose(pFile[1]);
}


void ifft_bitrev(void)
{
	   float dd;
 
	   long ia, ib, ik, ilk;
       long kk;
	   long nm1;
	   long md2, mm1;
	  
	   m=nsegment;
       mm1=m-1;
       dd=float( nstages );
       nm1=long(dd);
       md2=m/2;

  //    printf("\n nm1= %ld  md2=%ld  ns=%ld ", nm1, md2, ns);
	   
       for( i=1; i<= nm1; i++)   
       { 

		  for(ik=0; ik<m; ik++)
		  {
    		    ww[ik].r=y[ik].r;
				ww[ik].i=y[ik].i;
		  }
          
          kk=0;
          ia=0;

          ilk=long( pow(2,(i-1)) );

          for(j=1; j<=ilk; j++)
		  {
			  ib=ia+md2-1;
             
             for(k=ia; k<=ib; k++)
			 {
                if(k==ia)
				{
					kk=k;
				}
                
                y[k].r = ww[kk].r;
                y[k].i = ww[kk].i;
                
				y[k+md2].r = ww[kk+1].r;
				y[k+md2].i = ww[kk+1].i;
                
				kk=kk+2;
             }
             ia+= md2*2;
          
          }
          md2/=2;   
		  
        }
 }

void scaleth(void)
{

	ave=0.;
	float ms=0.;

	grmsout=0.;


	for(i=0; i<nout; i++)
	{
		ave+=y[i].r;
		ms+=pow(y[i].r,2);
	}
	ave/=float(nout);
	ms/=float(nout);

	grmsout=float(sqrt(ms));

	ratio=grms/grmsout;
	
	printf("\n grms = %14.7e   grmsout = %14.7e\n",grms,grmsout);
	printf("\n ratio = %14.7e \n",ratio);

	for(i=0; i<nout; i++)
	{
		y[i].r*=ratio;
	}
	
}
void start(void)
{
	long ielk = long( (log(dur/dt)/log(2.)))-1;

	ielk++;

	if( pow(2.,ielk) <= FFTMAX)
	{
		nsegment=long( pow(2.,ielk) );
	}
	else
	{
		nsegment= FFTMAX;
	}
	np=nsegment;
	ns=nsegment;
	nstages=long( (log(ns)/log(2))+0.5 );


	dur=ns/sr;
	df=float(1./dur);
	fmax= float(ns*df/2.);
		
	printf("\n df   = %14.4g Hz",df);
	printf("\n fmax = %14.4g Hz",fmax);

	printf("\n dur  = %14.4g sec\n",dur);

    printf("\n\n ns = %ld \n",ns);


	slope_calc();
	slope_calc_v();
	slope_calc_d();

	printf("\n Input PSD \n\n");
	printf("   acceleration = %10.4e grms\n",grms);
	printf("   velocity     = %10.4e inch/sec rms\n",vrms);
	printf("   displacement = %10.4e inch rms\n",drms);

        
}
void print_output(void)
{
	
	    double vv=0.;
		double dd=0.;


		vave=0.;
		dave=0.;

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

				ttt=kt*dt;

				kt++;


				if( ttt > tmax_old )
				{
//					printf("\n r1:  ttt=%12.4g  tmax_old=%12.4g \n",ttt,tmax_old);
					iflag=1;
					break;
				}
		
				fprintf(pFile[2],"%15.9e \t %14.7e\n",ttt,y[j].r);
		
				if(amax < y[j].r){ amax = y[j].r; }
				if(amin > y[j].r){ amin = y[j].r; }

				if(j==0)
				{
					velox[j]=float( ( y[j].r )*dt*386./2. +vsum );
				}
				else
				{
					velox[j]=float( ( y[j].r )*dt*386. +vsum);
				}

				vsum=velox[j];
				vave+=velox[j];
				vv+=pow(velox[j],2.);

		
		}


			vave/=float(j);
			vv/=float(j);

			double vstd = sqrt( vv - pow(vave,2.)); 



			float vratio = float(vrms/vstd);
//			printf("\n vratio = %12.4g \n",vratio);

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

 				ttt=j*dt + tlast;
			 
				if( ttt > tmax_old )
				{
//				    printf("\n r2:  ttt=%12.4g  tmax_old=%12.4g \n",ttt,tmax_old);
					iflag=1;
					break;
				}
			
				velox[j]-=vave;

//				velox[j]*=vratio;

				fprintf(pFile[3],"%15.9e \t %14.7e\n",ttt,velox[j]);

			
				if(vmax < velox[j]){ vmax = velox[j]; }
				if(vmin > velox[j]){ vmin = velox[j]; }


				if(j==0)
				{
					disp[j]=float(( velox[j] )*dt/2. +dsum);
				}
				else
				{
					disp[j]=float(( velox[j] )*dt +dsum);
				}
				dsum=disp[j];
				dave+=disp[j];
				dd+=pow(disp[j],2.);

//				printf(" %ld  %12.4g  %12.4g  \n",j,velox[j],disp[j]);
//                getch();
			}

			dave/=float(j);
			dd/=float(j);

			double dstd = sqrt( dd - pow(dave,2.)); 


			for(j=0; j< nout; j++)
			{
              
				ttt=j*dt + tlast;
				if( ttt > tmax_old )
				{
//				    printf("\n r3:  ttt=%12.4g  tmax_old=%12.4g \n",ttt,tmax_old);
					iflag=1;
					break;
				}
				disp[j]-=dave;

//			    disp[j]*=drms/dstd;

				fprintf(pFile[4],"%15.9e \t %14.7e           \n",ttt,disp[j]);
				fprintf(pFile[5],"%15.9e \t %14.7e \t %14.7e \n",ttt,disp[j],velox[j]);

				if(dmax < disp[j]){ dmax = disp[j]; }
				if(dmin > disp[j]){ dmin = disp[j]; }

			}

			tlast=ttt+dt;

}
void intlog(void)
{
	printf("\n intlog \n");

	float fq,f_temp[NUM],psd_temp[NUM];

	double N;

	double conv = pow(386.,2.);

/*
	freq[i]
	psd[i]
	      
	fmax=freq[i];
	npsd; 
*/
	for(i=0;i<NUM;i++)
	{
		f_temp[i]=freq[i];
		psd_temp[i]=psd[i];
	}

    long nnn = 2000;

	if( nnn > NUM ){nnn=long(NUM/2.);}

	df = float( ( fmax-freq[0])/double(nnn) );


    freq[0]=f_temp[0];
    psd[0]=psd_temp[0];

    fq = f_temp[0];


	vpsd[0]=float(conv*psd[0]/( pow( (tpi*fq),2.)));
	dpsd[0]=float( vpsd[0]/( pow( (tpi*fq),2.)) );


	long num=1;


	while(1)
	{
		fq +=df;

		if(fq > fmax ){break;}

		for(i=0;i<(npsd-1);i++)
		{
			if(fq==f_temp[i])
			{
                 psd[num]=psd_temp[i];
				vpsd[num]=float( conv*psd[num]/( pow( (tpi*fq),2.)) );
				dpsd[num]=float( vpsd[num]/( pow( (tpi*fq),2.)) );

				freq[num]=fq;

				num++;

				break;
			}
			if( fq > f_temp[i] && fq < f_temp[i+1])
			{
				N = log10(psd_temp[i+1]/psd_temp[i])/log10(f_temp[i+1]/f_temp[i]);

				 psd[num] = float( psd_temp[i]*pow( (fq/f_temp[i]), N) );
				vpsd[num]=float( conv*psd[num]/( pow( (tpi*fq),2.)) );
				dpsd[num]=float( vpsd[num]/( pow( (tpi*fq),2.)) );

				freq[num]=fq;

		        num++;

				break;
			}
		}
	}
	npsd=num;

	printf("\n  df=%12.4g  fmax=%12.4g  npsd=%ld \n",df,fmax,npsd);
    
}
void slope_calc_v()
{
	//  calculate slopes

    
	      for( i=0; i < npsd-1; i++)
		  {
			  if(freq[i] > 1.0e-12)
			  {
					slopev[i]=float(log( vpsd[i+1]/vpsd[i] )/log( freq[i+1]/freq[i] ));
			  }
			  else
			  {
				  slopev[i] = 0.;
			  }
//	          printf("\n freq[i] = %14.7e  slope = %14.7e\n",freq[i],slope[i]);
		  }

		float ra=0.;

		vrms=0.;

		for( i=0; i < npsd-1; i++)
		{
			if(slopev[i] < -1.0001 ||  slopev[i] > -0.9999 )
			{
				ra+= float( ( vpsd[i+1]*freq[i+1]- vpsd[i]*freq[i])/( slopev[i]+1.) );
			}
			else
			{
				if(freq[i] < 1.0e-12){ freq[i]= float(1.0e-12);}

				ra+= float( vpsd[i]*freq[i]*log( freq[i+1]/freq[i]) );
			}
		}
		vrms=sqrt(ra);

//		printf("\n input vrms = %14.7e\n",vrms);

}
void slope_calc_d()
{
	//  calculate slopes

    
	      for( i=0; i < npsd-1; i++)
		  {
			  if(freq[i] > 1.0e-12)
			  {
					sloped[i]=float( log( dpsd[i+1]/dpsd[i] )/log( freq[i+1]/freq[i] ) );
			  }
			  else
			  {
				  sloped[i] = 0.;
			  }
//	          printf("\n freq[i] = %14.7e  slope = %14.7e\n",freq[i],slope[i]);
		  }

		float ra=0.;

		drms=0.;

		for( i=0; i < npsd-1; i++)
		{
			if(sloped[i] < -1.0001 ||  sloped[i] > -0.9999 )
			{
				ra+= float( ( dpsd[i+1]*freq[i+1]- dpsd[i]*freq[i])/( sloped[i]+1.) );
			}
			else
			{
				if(freq[i] < 1.0e-12){ freq[i]= float(1.0e-12); }

				ra+= float( dpsd[i]*freq[i]*log( freq[i+1]/freq[i]) );
			}
		}
		drms=sqrt(ra);

//		printf("\n input drms = %14.7e\n",drms);

}
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,'0' ) == NULL && strrchr(string,'1' ) == NULL && strrchr(string,'2' ) == NULL &&  strrchr(string,'3' ) == NULL && strrchr(string,'4' ) == NULL && strrchr(string,'5' ) == NULL && strrchr(string,'6' ) == NULL && strrchr(string,'7' ) == NULL && strrchr(string,'8' ) == NULL && strrchr(string,'9' ) == NULL && strrchr(string,'0' ) == 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,'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);
	}
}