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


#define INPUTLABELS \
if(iunits==1)\
{\
if(idim==1)\
{\
printf(" Freq(Hz)  Level(G^2/Hz)    \n");\
}\
if(idim==2)\
{\
printf(" Freq(Hz)  Level((inches/sec)^2/Hz)    \n");\
}\
if(idim==3)\
{\
printf(" Freq(Hz)  Level(inches^2/Hz)    \n");\
}\
}\
if(iunits==2)\
{\
if(idim==1)\
{\
printf(" Freq(Hz)  Level((m/sec^2)^2/Hz)    \n");\
}\
if(idim==2)\
{\
printf(" Freq(Hz)  Level((mm/sec)^2/Hz)    \n");\
}\
if(idim==3)\
{\
printf(" Freq(Hz)  Level(mm^2/Hz)    \n");\
}\
}\


#define MINF 0.001




void fileoper(void);

void calc(double *level);

void read_data(void);

void input_setup(void);

void input_check(void);

void acceleration_vd(void);

void header(void);



long numBytes=300;

char string[300];


double level[5000];
double f[5000],a[5000],s[5000],v[5000],d[5000];
double db, ra,grms;
double oa, ov, od;
double aa,ff;

int iflag;
int i, ic, id, ifile, iunits, idim;
int m, n, nn;

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



FILE *pFile[8];

char filename[8][FILENAME_MAX];



int main( int argc, char *argv[] )
{

		printf("                               \n");
		printf(" psdint.cpp ver 3.3, April 17, 2010   \n");
		printf("							   \n");
		printf(" by Tom Irvine                 \n");
		printf("    Email:  tomirvine@aol.com  \n");
		printf("							   \n");
		printf(" This program accepts an acceleration, velocity, or displacement \n");
		printf(" power spectral density function.  Given one of these functions, the \n");
		printf(" program then integrates or differentiates as appropriate to find \n");
		printf(" the other two functions.\n");

		printf("\n The program also calculates the overall level of each of the three \n");
		printf(" functions. \n\n");

		printf("\n Note that the initial input frequency must be > 0. ");
		printf("\n Please edit the input file accordingly. \n\n");



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

			iunits=1;

			id=1;

			idim=1;

			ifile=1;

			printf("\n");

		}





		ic=1;



		while(ic==1)
        {

		

			input_setup();

		    fileoper();



			input_check();



			acceleration_vd();

//
// Acceleration
//
		    calc(a);
		    oa=grms;

			
			if(iunits==1)
			{
				printf("\n\n Acceleration Power Spectral Density\n");
				printf("\n\tFreq\tAPSD\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t(G^2/Hz)\tlog-log\t\t(dB/octave) \n");
				printf("\t%6.1lf\t%11.4e \n", f[0],a[0]);    
				fprintf(pFile[4],"\n\n Acceleration Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tAPSD\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t(G^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],a[0]);

			    fprintf(pFile[1],"\t%6.1lf\t%11.4e \n", f[0],a[0]);
				
				for( i=0; i < m-1; i++)
				{
					db=( 10.*log10(2.) )*s[i];
					printf("\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],a[i+1],s[i],db);   
					fprintf(pFile[4],"\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],a[i+1],s[i],db);   
			        fprintf(pFile[1],"\t%11.4e\t%11.4e \n",f[i+1],a[i+1]);		

				}

			 }
			 else
			 {
				printf("\n\n Acceleration Power Spectral Density\n");
				printf("\n\tFreq\tAPSD\t\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t((m/sec^2)^2/Hz)\tlog-log\t\t(dB/octave) \n");
				printf("\t%6.1lf\t%11.4e \n", f[0],a[0]);    
				fprintf(pFile[4],"\n\n Acceleration Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tAPSD\t\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t((m/sec^2)^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],a[0]);

			    fprintf(pFile[1],"\t%6.1lf\t%11.4e \n", f[0],a[0]);

				for( i=0; i < m-1; i++)
				{
					db=( 10.*log10(2.) )*s[i];
					printf("\t%6.1lf\t%11.4e\t\t%11.4e\t%11.4e  \n", f[i+1],a[i+1],s[i],db);   
					fprintf(pFile[4],"\t%6.1lf\t%11.4e\t\t%11.4e\t%11.4e  \n", f[i+1],a[i+1],s[i],db); 

			        fprintf(pFile[1],"\t%11.4e\t%11.4e \n",f[i+1],a[i+1]);		
				}

			 }
	         
//
// Velocity
//
			 if(iunits==1)
			 {
				printf("\n\n Velocity Power Spectral Density\n");
				printf("\n\tFreq\tVPSD\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t((IN/S)^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\n\n Velocity Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tVPSD\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t((IN/S)^2/Hz)\tlog-log\t\t(dB/octave) \n");

			 
				if( idim !=2){ v[0]=(pow( 386.,2.) )*a[0]/(pow( ( tpi*f[0] ),2.)) ;}
				printf("\t%6.1lf\t%11.4e \n", f[0],v[0]);   
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],v[0]);    

				if(idim !=2)
				{
					for(i=0; i<m; i++)
					{
						v[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),2.)); 	 
					}
				}
			 }
			 else
			 {
				printf("\n\n Velocity Power Spectral Density\n");
				printf("\n\tFreq\tVPSD\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t((mm/sec)^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\n\n Velocity Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tVPSD\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t((mm/sec)^2/Hz)\tlog-log\t\t(dB/octave) \n");

			 
				if( idim !=2){ v[0]=(pow( 1000.,2.) )*a[0]/(pow( ( tpi*f[0] ),2.)) ;}
				printf("\t%6.1lf\t%11.4e \n", f[0],v[0]);   
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],v[0]);    

				if(idim !=2)
				{
					for(i=0; i<m; i++)
					{
						v[i]=(pow( 1000.,2.) )*a[i]/(pow( ( tpi*f[i] ),2.)); 	 
					}
				}
			 }


	         calc(v);
             ov=grms;

		//	 printf("\n istart=%ld m-1=%ld",istart, m-1);

			 for( i=0; i < m-1; i++)
			 {
                  db=( 10.*log10(2.) )*s[i];
                  printf("\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],v[i+1],s[i],db);    
				  fprintf(pFile[4],"\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],v[i+1],s[i],db); 
			 }

//
// Displacement
//
			 if(iunits==1)
			 {
				printf("\n\n Displacement Power Spectral Density\n");
				printf("\n\tFreq\tDPSD\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t(IN^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\n\n Displacement Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tDPSD\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t(IN^2/Hz)\tlog-log\t\t(dB/octave) \n");

				if( idim != 3){ d[0]=(pow( 386.,2.) )*a[0]/(pow( ( tpi*f[0] ),4.)); }
				printf("\t%6.1lf\t%11.4e \n", f[0],d[0]); 
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],d[0]);   
			 
				for(i=0; i<m; i++)
				{
					if( idim != 3)
					{
						if( f[i] >= MINF)
						{
							d[i]=(pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.));
						}
					}
				}
			 }
			 else
			 {
				printf("\n\n Displacement Power Spectral Density\n");
				printf("\n\tFreq\tDPSD\t\tSlope\t\tSlope  \n");
				printf("\t(Hz)\t(mm^2/Hz)\tlog-log\t\t(dB/octave) \n");
				fprintf(pFile[4],"\n\n Displacement Power Spectral Density\n");
				fprintf(pFile[4],"\n\tFreq\tDPSD\t\tSlope\t\tSlope  \n");
				fprintf(pFile[4],"\t(Hz)\t(mm^2/Hz)\tlog-log\t\t(dB/octave) \n");

				if( idim != 3){ d[0]=(pow( 1000.,2.) )*a[0]/(pow( ( tpi*f[0] ),4.)); }
				printf("\t%6.1lf\t%11.4e \n", f[0],d[0]); 
				fprintf(pFile[4],"\t%6.1lf\t%11.4e \n", f[0],d[0]);   
			 
				for(i=0; i<m; i++)
				{
					if( idim != 3)
					{
						if( f[i] >= MINF)
						{
							d[i]=(pow( 1000.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.));
						}
					}
				}
			 }


	         calc(d);
             od=grms;

			 for( i=0; i < m-1; i++)
			 {
                  db=( 10.*log10(2.) )*s[i];
                  printf("\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],d[i+1],s[i],db); 
				  fprintf(pFile[4],"\t%6.1lf\t%11.4e\t%11.4e\t%11.4e  \n", f[i+1],d[i+1],s[i],db);    
			 }
//
//
//
			 for(i=0; i < m; i++)
			 {
				 fprintf( pFile[2],"\t%11.4e\t%11.4e \n",f[i],v[i]);
				 fprintf( pFile[3],"\t%11.4e\t%11.4e \n",f[i],d[i]);
             }
			 fclose( pFile[1] );
			 fclose( pFile[2] );
			 fclose( pFile[3] );
//
//
            char alabel1[20],alabel2[20],vlabel[20],dlabel[20];

			if(iunits==1)
			{
				strcpy(alabel1,"GRMS");
				strcpy(alabel2,"G 3-sigma");
				strcpy(vlabel,"in/sec");
				strcpy(dlabel,"inch");
			}
			else
			{
				strcpy(alabel1,"(m/sec^2) RMS");
				strcpy(alabel2,"(m/sec^2) 3-sigma");
				strcpy(vlabel,"mm/sec");
				strcpy(dlabel,"mm");
			}

				printf("\n\n Overall Acceleration = %12.4g %s ",oa,alabel1);
			      printf("\n                      = %12.4g %s \n\n",3.*oa,alabel2);

				printf(" Overall Velocity     = %12.4g (%s) RMS \n",ov,vlabel);       
				printf("                      = %12.4g (%s 3-sigma) \n\n",ov*3,vlabel); 
					
				printf(" Overall Displacement = %12.4g (%s RMS) \n",od,dlabel);
				printf("                      = %12.4g (%s 3-sigma) \n",od*3,dlabel); 

                if(iunits==1)
			    {
                    od*=25.4;         
                    strcpy(dlabel,"mm");         
                    printf("\n\n ** metric ** \n\n");
   				    printf(" Overall Displacement = %12.4g (%s RMS) \n",od,dlabel);
				    printf("                      = %12.4g (%s 3-sigma) \n",od*3,dlabel);                     
                 }	

			 printf("\n  Output files \n");
			 printf("\n  Complete Text: %s ",filename[4]); 
             printf("\n  Acceleration:  %s ",filename[1]);
             printf("\n  Velocity    :  %s ",filename[2]);
             printf("\n  Displacement:  %s \n",filename[3]);

             printf("\n  Another calculation?  1=yes 2=no     \n");
             scanf("%d",&ic);    



			 ifile=0;

	}
}
	   
void calc(double *level)
{
		  ra=0.;
          grms=0.;

		  for(i=0; i < 5000; i++)
		  {
			  s[i]=0.;
		  }
		  

//  calculate slopes

          nn=m-1;

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

			  s[i]=log( level[i+1]/level[i] )/log( f[i+1]/f[i] );

              if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			  {	   
                 ra+= ( level[i+1] * f[i+1]- level[i]*f[i])/( s[i]+1.);
			  }
		      else
			  {
                 ra+= level[i]*f[i]*log( f[i+1]/f[i]);
			  }
		  }
          grms=sqrt(ra);
}

void fileoper(void)
{
       strcpy(filename[1], "a.psd");
       pFile[1]=fopen(filename[1], "w");

	   strcpy(filename[2], "v.psd");
       pFile[2]=fopen(filename[2], "w");

	   strcpy(filename[3], "d.psd");
       pFile[3]=fopen(filename[3], "w");

	   strcpy(filename[4], "psdint.txt");
       pFile[4]=fopen(filename[4], "w");
}

void read_data(void)

{

        if(ifile==1)

		{

			idim=1;

		}



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

		

		while(pFile[0] == NULL )

		{

			printf("\n Failed to open file: %s \n", filename[0]);

				

			printf("\n Please enter the input filename: \n");



			scanf("%s",filename[0]);

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

		}

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



		header();



		i=0;



//		while( fscanf(pFile[0], "%lf %lf", &ff, &aa) > 0 && iflag < 999)



		while( fgets(string,numBytes,pFile[0] )>0 && iflag < 999)

		{



			if(strrchr(string,',' ) != NULL )
			{

				sscanf(string,"%lf, %lf", &ff,&aa );

			}

			else
			{

				sscanf(string,"%lf %lf", &ff,&aa );

			}



//				printf("\n %lf %lf \n",ff,aa); //***



			if(aa<=0)
			{



				printf("\n %lf %lf \n",ff,aa);



				printf("\n Input error:  all amplitudes must be > 0 \n");



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



				getch();



				exit(1);



			}

			if( ff > MINF )

			{

				f[i]=ff;



				if(idim==1){a[i]=aa;}

				if(idim==2){v[i]=aa;}

				if(idim==3){d[i]=aa;}



//				printf("\n %lf %lf %d \n",ff,a[i],idim); //***



				i++;

			}

			if(i>=5000)break;  

		}

		fclose( pFile[0]);



		m=i;

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

		printf(" \n");

}

void input_setup(void)
{

		iflag=0;



		if(iunits !=1 )
		{

			iunits=99;



			while( iunits !=1 && iunits!=2 && ifile !=1 )
			{

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

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

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


				scanf("%d",&iunits);

				printf("\n");

			}

		}





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

		{

			f[i]=0.;

			a[i]=0.;

			s[i]=0.;

			v[i]=0.;

			d[i]=0.;

		}

		oa=0.;

		ov=0.;

		od=0.;



		int id=99;

   

		while( id !=1 && id !=2 && id !=3 && ifile !=1)

		{

			printf( " \n\n Select input data dimension. \n" );

			printf( " 1= acceleration power spectral density       \n");

			printf( " 2= velocity power spectral density       \n");

			printf( " 3= displacement power spectral density       \n");



			scanf("%d",&id);

		}



		idim = id;



		int im=99;





		if(ifile == 1 )

		{

			INPUTLABELS

			read_data();				

        }

		else
		{

			while( im !=1 && im !=2 )
			{

				printf( " \n\n Select data input method \n" );

				printf( " 1=keyboard  2=file       \n");

	     

				scanf("%d",&im);



			}

			const int imethod = im;          



			iflag=0;

       

			if(imethod==1)

			{	    

				printf( " How Many Breakpoints (minimum=2)? \n");

          		scanf("%d",&n);



				if(n < 2)
				{

						printf( "Input error \n");

						iflag=950;

				}

				else
				{

						printf(" Enter Points (free-format) \n");

                  	  

						INPUTLABELS



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

							if(idim==1){scanf("%lf %lf",&f[i],&a[i]);}

							if(idim==2){scanf("%lf %lf",&f[i],&v[i]);}

							if(idim==3){scanf("%lf %lf",&f[i],&d[i]);}

						}
				}

				m=n;

			}

			else
			{

				printf( " The file must contain two columns: \n");



				printf( " \n An error message will occur if the file does not exist.\n");

				printf( " Input filename \n");

			 

			    scanf("%s",filename[0]);



				INPUTLABELS

				read_data();

			}

		}

   		  

}

void input_check(void)

{

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

			 {

//				printf("\n check loop ");



                if((i >0) && (f[i] < f[i-1]))

				{

                   printf(" Error:  frequencies must ascend. \n");

				}

                if(f[i] < MINF)

				{

                   printf(" Error: frequencies must be > %7.3f \n",MINF);

                   printf(" Input frequency = %lf \n", f[i]);



				   iflag=981;

				} 

                

				

				if(idim==1)

				{

				   if(a[i] <= 0.)

				   {

                      printf("Error: all input numbers must be > 0.\n");

                      printf(" %lf %lf \n",f[i],a[i]);

				      iflag = 981;

				   }

//			      fprintf( pFile[1],"\t%11.4e\t%11.4e \n",f[i],a[i]);

                }

	            if(idim==2)

				{

				   if(v[i] <= 0.)

				   {

                      printf("Error: all input numbers must be > 0.\n");

                      printf(" %lf %lf \n",f[i],v[i]);

				      iflag = 981;

				   }

//			      fprintf( pFile[1],"\t%11.4e\t%11.4e \n",f[i],v[i]);

                }

				if(idim==3)

				{

				   if(d[i] <= 0.)

				   {

                      printf("Error: all input numbers must be > 0.\n");

                      printf(" %lf %lf \n",f[i],d[i]);

				      iflag = 981;



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

					  getch();

					  exit(1);

				   }

//			      fprintf( pFile[1],"\t%11.4e\t%11.4e \n",f[i],d[i]);

                }



			 } 

		

}

void acceleration_vd(void)

{

	

			 if(idim==2)

			 {

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

				 {

                    if( f[i] >= MINF)

					{

						if(iunits==1)

						{

							a[i]=v[i]*(pow( ( tpi*f[i] ),2.))/(pow( 386.,2.) );

						}

						else

						{

							a[i]=v[i]*(pow( ( tpi*f[i] ),2.))/(pow( 1000.,2.) );

						}

					}       

				 }

			 }

			 if(idim==3)

			 {

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

				 {

                    if( f[i] >= MINF)

					{

						if(iunits==1)

						{

							a[i]=d[i]*(pow( ( tpi*f[i] ),4.))/(pow( 386.,2.) );

						}

						else

						{

							a[i]=d[i]*(pow( ( tpi*f[i] ),4.))/(pow( 1000.,2.) );

						}

					}       

				 }

			 }

		     

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

	}

}