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

#define MINF 0.001
#define LIMIT 10000

double level[LIMIT];
double f[LIMIT],a[LIMIT],v[LIMIT],d[LIMIT];
double db, ra,rd,rv,grms,drms,vrms;
double oa, ov, od;
double aa,ff;

double df;

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

int iflag;
int i, ic, imethod;
int m, n, nn;

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

void fileoper(void);
void calc(double *level);

int mflag=0;

void main( int argc, char *argv[] )	
{
	if(argc > 1 )
	{
		printf("\n\n  argc = %d ", argc );
		printf("\n argv[0] %s ", argv[0] );
		printf("\n argv[1] %s ", argv[1] );

		strcpy (filename[0], argv[1] );

		mflag=1;

	}
	   
			printf("\n  psdtrap, ver 2.1,   October 31, 2008  \n");
			printf("  by Tom Irvine                \n");
			printf("  Email:  tomirvine@aol.com  \n\n");
  
			printf(" Integration of an acceleration power spectral \n");
			printf(" density curve to determine overall level. \n\n");
			printf(" The power spectral density must have a constant  \n");
			printf(" frequency resolution, delta f.  \n\n");

            for(i=0; i<LIMIT; i++)
			{
		       f[i]=0.;
               a[i]=0.;
			   v[i]=0.;
			   d[i]=0.;
            
			}
	        oa=0.;
	     
			if(mflag==0)
			{
       
				printf( " The file must contain two columns: \n");
				printf( " Freq(Hz)  Level(G^2/Hz)    \n");

				printf( " \n An error message will occur if the file does not exist.\n");
				printf( " Input filename \n");
			  
				scanf("%s",filename[0]);
			}

			 
			 pFile[0]=fopen(filename[0], "rb");

			 if(pFile[0] == NULL)
			 {
				 printf( "\n Input file error \n");

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

			 i=0;
			 while( fscanf(pFile[0], "%lf %lf", &ff, &aa) > 0  )
			 {
	
				 	 f[i]=ff;
					 a[i]=aa;

					 i++;

                     if(i>=LIMIT)break;  
			 }
			 fclose( pFile[0]);

			 m=i;
			 printf("\n m=%d ", m);



			 df = ( f[m-1] - f[0] )/( unsigned long (m-1) );
			 printf("\n\n df= %10.3e Hz\n", df);


			 ra=0.;
			 rv=0.;
			 rd=0.;
	
			 if(f[0]> MINF)
			 {
				v[0]=a[0]/pow((tpi*f[0]),2.);
				d[0]=a[0]/pow((tpi*f[0]),4.);
             }
			 
			 ra+=a[0]*df/2.;
			 rv+=v[0]*df/2.;
			 rd+=d[0]*df/2.;
			 

	         for(i=1; i<m; i++)
			 {
				 
					 v[i]=a[i]/pow((tpi*f[i]),2.);
				     
					 d[i]=a[i]/pow((tpi*f[i]),4.);

					 ra+=a[i]*df;

					 rv+=v[i]*df;

					 rd+=d[i]*df;
			 }


			 grms = sqrt( ra );

			 vrms = 386.*sqrt( rv );

			 drms = 386.*sqrt( rd );
		  	  
             printf("\n\n Overall Acceleration = %15.4g GRMS ",grms);               
             printf("\n\n Overall Velocity     = %15.4g (inch/sec) RMS \n",vrms);
			 printf("                      = %15.4g (mm/sec) RMS ",vrms*25.4); 			 
             
			 printf("\n\n Overall Displacement = %15.4g inch RMS \n",drms); 
                 printf("                      = %15.4g mm RMS ",drms*25.4);  
			 
			printf("\n\n Program complete. \n Press any key to exit.");
			getch();
}