/*
	program qsrs.cpp

	by Tom Irvine   Email: tomirvine@aol.com
*/


#include <iostream>
#include <fstream>
#include <cmath>
#include <cstring>
#include <string>

#include <stdlib.h>
#include <stdio.h>

using namespace std;


#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];\
		tmax[j]=t;\
	}\
	if(x[j] < xmin[j])\
	{\
		xmin[j]=x[j];\
		tmin[j]=t;\
	}\
 \
	xbb[j]=xb[j];\
	xb[j]= x[j];\
	ybb[j]=yb[j];\
	yb[j]= y;\
 \
}



void line_test(void);

void linear_interpolation(const double &scale);
void files(void);
void read_data_check(void);
void initialize_filter_coefficients(const int & algor, const double & dam, const double fn[]); //  alias method
void select_algorithm(int & algor, double & dam);
void select_interpolation(const double &scale);
void initialize_frequencies(const double &ssr);
void initialize_arrays(void);
void output_files(void);
void primary(const double &scale);
void residual(const double &dt,const double fn[]);
void time_parameters();

void line_test_read(const double &scale);

#define numBytes 300

const long MAX=3000000;
const long NUM=20000;


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

char string[numBytes];

double a1[NUM],a2[NUM];
double b1[NUM],b2[NUM],b3[NUM];


double dt;

double fn[NUM];

double scale;

double sr;
double t, t1, tb;
double tlast;

double tmax[NUM],tmin[NUM];

double xb[NUM], xbb[NUM];

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

double y,yone,ytwo;

double yb[NUM],ybb[NUM];



double tt[MAX],aa[MAX];

long number_lines;

double peak_f=0.;
double peak_m=0.;

int ire;
int mflag=0;
int nmax=0;

long idamp,ioct;
long last;

int linear;
int kflag;

unsigned long i,j,jnum;

ifstream inFile;

char filename[12][FILENAME_MAX];
ofstream outFile[12];


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

	if(argc > 1 )
	{

			cout << "\n\n  argc = " << argc << endl;
			cout << "\n argv[0] = " << argv[0] << endl;
			cout << "\n argv[1] = " << argv[1] << endl;

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

			mflag=1;


	}

	cout << "\n qsrs.cpp, ver 7.5,  May 29, 2010 " << endl;
	cout << "  \n  by Tom Irvine " ;
	cout << "  \n  Email:  tomirvine@aol.com \n";

	cout << "\n  Shock response spectrum of an arbitrary acceleration " ;
	cout << "\n  base input time history supplied by the user.\n" << endl;

    int algorithm;
    double damp;

    initialize_arrays();

	files();

	line_test();

	time_parameters();

	select_interpolation(scale);

	select_algorithm(algorithm,damp);

    initialize_frequencies(sr);

	initialize_filter_coefficients(algorithm,damp,fn);



	primary(scale);

//	cout<<"\n ire= "<<ire<<endl;

    if(ire==1)
    {
       residual(dt,fn);
    }
	output_files();

	cout<< "\n Press any key to exit.\n" << endl;
	cin.get();
	cin.get();

	exit(1);

}
void linear_interpolation(const double &scale)
{

//	cout << "\n\n Linear Interpolation \n";

	double aaa;
	double ttt;
	double c1,c2;

	double x,l;

	long i,j,jfirst;

	long last,nt;

     if( inFile.is_open() )
     {
        inFile.clear();
        inFile.seekg(0,std::ios::beg);
        inFile.close();
     }
	 inFile.open(filename[0]);

     while( !inFile.is_open() )
     {
				cout << "\n R1:  Failed to open file: " << filename[0];
				cout << "\n\n Please enter the input filename: \n";

                cin >> filename[0];
				inFile.open(filename[0]);
     }
     cout << " File: " << filename[0] << " opened. \n";


	 cout << "\n Enter output time history filename for interpolation. " << endl;
     cin >> filename[1];

     outFile[1].open(filename[1]);
     outFile[1].precision(6);

	cout << "\n\n Enter Output Sampling Rate (samples/sec): \n" << endl;

    cin >> sr;

	dt=1./sr;

	i=0;

    char string[numBytes];

	while(1)
	{
        inFile.getline(string,numBytes);
        if(inFile.eof()){break;}

		kflag=0;


		if(kflag != 0)
		{
		}
		else
		{
			if(strrchr(string,',' ) != NULL )
			{
				sscanf(string,"%lf, %lf", &tt[i],&aa[i] );
			}
			else
			{
				sscanf(string,"%lf %lf", &tt[i],&aa[i] );
			}

			tt[i]*=scale;
		}

		i++;

		if(i>=MAX){break;}


	}

	last=i-1;

    inFile.clear();
    inFile.seekg(0,std::ios::beg);
	inFile.close();


	nt=long( (tt[last]-tt[0])/dt );

    cout << "\n t[" << last << "]=" << tt[last] <<"   t[0]=" << tt[0] <<endl;

    cout << "\n last = " << last << endl;

    cout << "\n nt= " << nt << endl;


	if(nt<1)
	{
		cout <<  "\n nt error \n" << endl;

		exit(1);
	}

    cout << "\n dt= " << dt << endl;

	outFile[1] << tt[0] << "\t" << aa[0] << endl;

	cout << endl;

    jfirst=1;


	for(i=1; i<=nt; i++)
	{
		ttt = dt*i+tt[0];

//		cout<<" ttt= "<<ttt<<endl;

		for(j=jfirst; j<=last; j++)
		{

			if( (j-1) >= 0 && ttt==tt[j-1])
			{
                outFile[1] << ttt << "\t" << aa[j-1] << endl;

				break;
			}
			else
			{
				if( j >= 0 && ttt==tt[j])
				{
					outFile[1] << ttt << "\t" << aa[j] << endl;
					break;
				}
				else
				{
					if( ttt> tt[j-1] && ttt < tt[j])
					{
						x=ttt-tt[j-1];
						l=tt[j]-tt[j-1];

                        c2=x/l;
						c1=1.-c2;

						aaa=c1*aa[j-1] + c2*aa[j] ;

					    outFile[1] << ttt << "\t" << aaa << endl;
						break;
					}
				}
			}
		}

		jfirst=j-1;

		if(jfirst <0 ){jfirst = 0;}

    }
    outFile[1].close();



}
void files()
{

	if(mflag==0)
	{

		cout << "\n Input must be time and accel(G) " << endl;
		cout << "\n Please enter the input filename. \n" << endl;

        cin >> filename[0];
	}

    if( inFile.is_open() )
    {
        inFile.clear();
        inFile.seekg(0,std::ios::beg);
        inFile.close();
    }
	inFile.open(filename[0]);

	while(!inFile.is_open() )
	{
		cout << "\n R2: Failed to open file: " << filename[0] << endl;

		cout << "\n Please enter the input filename: " << endl;

		cin >> filename[0];

        inFile.open(filename[0]);
	}

	cout << " File: " << filename[0] << " opened " << endl;

	strcpy(filename[2], "srs.out");

    outFile[2].open(filename[2]);

	int itu;

	scale=1.;

	cout << "\n Select input file time unit \n 1= sec   2= msec " << endl;

//	scanf("%d",&itu);

    cin >> itu;

	if(itu==2)
	{
		scale=0.001;
	}

    cout << "\n Include residual? ";
    cout << "\n   1=yes  2=no \n";

    cin >> ire;

}
void initialize_filter_coefficients(const int & algor, const double & dam, const double fn[])
{
    double omega;
    double omegad;

	if(algor==1)
	{

        double cosd;
        double sind;
        double domegadt;


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

			omega=2.*pi*fn[j];

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

			cosd=cos(omegad*dt);

			sind=sin(omegad*dt);

			domegadt=dam*omega*dt;

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

			b1[j]=2.*domegadt;
			b2[j]=omega*dt*exp(-domegadt);
			b2[j]*=( (omega/omegad)*(1.-2.*pow(dam,2))*sind -2.*dam*cosd );

			b3[j]=0.;

		}
	}
	else
	{

        double C,E,K,S,Sp;

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

			omega=2.*pi*fn[j];

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

			E=exp(-dam*omega*dt);
			K=omegad*dt;
			C=E*cos(K);
			S=E*sin(K);

			Sp=S/K;

			a1[j]=2*C;
			a2[j]=-pow(E,2);

			b1[j]=1.-Sp;
			b2[j]=2.*(Sp-C);
			b3[j]=pow(E,2)-Sp;

		}
   }
}
void select_algorithm(int & algor, double & dam)
{

	cout<< "\n Enter damping format. \n 1= damping ratio   2= Q \n";

    cin >> idamp;


	if(idamp == 1 )
	{
	     cout<<"\n Enter damping ratio (typically 0.05) ";
		 cin >> dam;
	}
	else
	{
	     cout<<"\n Enter amplification factor Q (typically 10) ";
         cin >> dam;


	     dam = 1./(2.*dam);
	}

	if(dam >= 1.0)
	{
		  cout << "\n Error:  damping value must be < 1 \n";

		  exit(1);
	}



	cout << "\nSelect alogrithm:";

	cout << "\n 1=Kelly-Richman  2=Smallwood \n";

    cin >> algor;

	outFile[2] << "\n  SHOCK RESPONSE SPECTRUM OUTPUT \n";
    outFile[2] << "\n  Input file = " << filename[0] << endl;
    outFile[2] << "dam= " << dam << endl;
    outFile[2] << "\n  dt= "<< dt << " sr= " << sr << " samples/sec" << endl;

}
void select_interpolation(const double &scale)
{
	cout << "\n\n Perform linear interpolation on time history? ";
	cout << "\n 1=yes  2=no \n";

	cin >> linear;

	if(linear==1)
	{
		linear_interpolation(scale);
    }
}

void initialize_frequencies(const double &ssr)
{
    double bw;

	while(1)
	{
		cout << "\n\n Enter starting frequency (Hz) \n";

        cin >> fn[0];

		if(fn[0] > ssr/4.)
		{
			cout << "\n\n Starting frequency must be < " << ssr/4. << endl;
		}
		else
		{
			if(fn[0] <= 1.0e-20 )
			{
				cout << "\n Starting frequency must be > 0. " <<endl;
			}
			else
			{
				break;
			}
		}
	}

   while(1)
   {
		cout << "\n\n Enter frequency spacing option \n";

		cout << "\n 1= 1/12 octave   2= constant bandwidth \n";

        cin >> ioct;

		if(ioct==2)
		{
            cout << "\n\n Enter bandwidth (Hz) \n";

			cin >> bw;
		}

		if(ioct==1 || ioct==2 ){break;}

   }

   for(j=1; j< NUM; j++)
   {

		if(ioct == 1)
		{
			fn[j]=fn[0]*pow(2.,j*(1./12.));
		}

		if(ioct == 2)
		{
			fn[j]=fn[j-1]+bw;
		}

		if( fn[j] > ssr/8.)break;

   }
   jnum=j;

}
void initialize_arrays(void)
{

	for(i=0; i< NUM; ++i)
	{
		xmax[i]=0.;
		xmin[i]=0.;

		tmax[i]=0.;
		tmin[i]=0.;

		x[i]=0.;
		xb[i]=0.;
		xbb[i]=0.;

		yb[i]=0.;
		ybb[i]=0.;
	}

}

void output_files()
{
	/**************************** print response spectra *****************/

   strcpy(filename[1], "srs.grp");

//  outFile[1]=fopen(filename[1], "w");
   outFile[1].open(filename[1]);

   strcpy(filename[3], "srs2.grp");

//   outFile[3]=fopen(filename[3], "w");
    outFile[3].open(filename[3]);

   outFile[2] << "\n total samples= " << i << endl;

   outFile[2] << "\n duration= " << i*dt << " sec" << endl;

   outFile[2] << "\n primary start time= " << t1 <<" sec \n primary end time  = " << tlast << "sec\n";

   outFile[2] << "\n     fn(Hz)      pos srs(G)      neg srs(G)    pos time(sec)    neg time(sec) \n";


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

        outFile[2] << fn[j] <<"\t"<< xmax[j] <<"\t"<< xmin[j] <<"\t"<< tmax[j] <<"\t"<< tmin[j] <<endl;


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

        outFile[3] <<fn[j] <<"\t"<< xmax[j] <<"\t"<< xmin[j] <<endl;

		if(xmax[j]<fabs(xmin[j]))
		{
			xmax[j]=fabs(xmin[j]);
		}

		outFile[1]  <<fn[j] <<"\t"<< xmax[j] <<endl;


		if(xmax[j]>peak_m)
		{
			peak_m=xmax[j];
			peak_f=fn[j];
		}
   }

   outFile[1].close();
   outFile[2].close();
   outFile[3].close();


   cout << "\n\n Absolute peak = " << peak_m << " G at " << peak_f << " Hz \n";
   cout << "\n\n Output files: " << endl;
   cout << "\n  srs.out - text output \n" << endl;

   cout << "\n srs.grp  - natural frequency (Hz) " << endl;
   cout << "\n             and peak absolute acceleration (G) " << endl;

   cout << "\n srs2.grp - natural frequency (Hz)  " << endl;
   cout << "\n             and peak positive acceleration (G)  ";
   cout << "\n             and peak negative acceleration (G)  \n";

   cout << "\n Please call the *.grp files into your own graphics program. " << endl;
   cout << "\n Note that the data files are tab delimited. " << endl;


}
void primary(const double &scale)
{
	/********************** primary response calculation *******************/

	cout << "\n Calculating primary response..... \n\n";



    if( inFile.is_open() )
    {
        inFile.clear();
        inFile.seekg(0,std::ios::beg);
        inFile.close();
    }


	if(linear==2)
	{
		inFile.open(filename[0]);
	}
	else
	{
		inFile.open(filename[1]);
	}

    if( !inFile.is_open() )
    {
         cout<<"\n linear = " << linear << endl;

        if(linear==2)
        {
	       cout << "\n R4: Failed to open file: " << filename[0];
        }
        else
        {
 	       cout << "\n R5: Failed to open file: " << filename[1];
        }

    }
    else
    {
   	    if(linear==2)
	    {
           cout << " File: " << filename[0] << " opened. \n";
        }
        else
        {
           cout << " File: " << filename[1] << " opened. \n";
        }
    }

	i=0;

    char string[numBytes];

    cout<<"\n nmax= "<<nmax<<endl;

    for(i=0;i<nmax;i++)
    {
        inFile.getline(string,numBytes);
    }

    while(1)
    {
        inFile.getline(string,numBytes);

        if(inFile.eof())
        {
//            cout<<"break"<<endl;
            break;
        }


		kflag=0;

		if(kflag != 0)
		{
		}
		else
		{
			if(strrchr(string,',' ) != NULL )
			{
				sscanf(string,"%lf, %lf", &t,&y );
			}
			else
			{
				sscanf(string,"%lf %lf", &t,&y );
			}

			t*=scale;
		}


		for(j=0; j<=jnum; ++j)
		{
			CORE
		}
		i++;

	}
//    inFile.clear();
//    inFile.seekg(0,std::ios::beg);
	inFile.close();

    cout<<"\n end primary "<<endl;

}
void residual(const double &dt,const double fn[])
{
	/********************** residual response calculation *******************/


   cout<<"\n Calculating residual response..... \n";

   tlast=t;

   long ijk;

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

   y=0.;


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

	    t=tlast+dt*ijk;

		for(j=0; j<=jnum; ++j)
		{
			if( dt*ijk > 1.2/fn[j] ){break;}

			CORE

		}
   }

}
void time_parameters()
{

	long j=0;
	long ijk=0;

	double t1=0.;
	double t,a;

	double tb=-1.0e+90;

	double dtmin=1.0e+90;
	double dtmax=-dtmin;

    number_lines=0;

	char string[numBytes];


     if( inFile.is_open() )
     {
        inFile.clear();
        inFile.seekg(0,std::ios::beg);
        inFile.close();
     }
	 inFile.open(filename[0]);

     while( !inFile.is_open() )
     {
				cout << "\n R1:  Failed to open file: " << filename[0];
				cout << "\n\n Please enter the input filename: \n";

                cin >> filename[0];
				inFile.open(filename[0]);
     }


    for(i=0;i<nmax;i++)
    {
        inFile.getline(string,numBytes);
    }

	while(1)
	{

        inFile.getline(string,numBytes);


        if(inFile.eof())
        {
            cout<<"\n break \n";
                        break;
        }


		kflag=0;

		if(kflag !=0)
		{
			ijk++;

			if(ijk>20)
			{
				cout<<"\n excessive number of line_test lines \n";

				exit(1);
			}
		}
		else
		{
			i=number_lines;

			if(strrchr(string,',' ) != NULL )
			{
				sscanf(string,"%lf, %lf", &t,&a );
			}
			else
			{
				sscanf(string,"%lf %lf", &t,&a );
			}

            t*=scale;


            if(number_lines==0)
			{
				t1=t;
			}
			if(t<tb && number_lines>=1)
			{
				cout<<"\n Error: time values must ascend. \n";

				cout << "\n D: " <<string;

                cout << "\n t= "<<t<< " tb= " << tb <<endl;

				cout<<"\n Press any key to exit  \n";
				cin.get();
			}

			if(number_lines>=1)
			{
				if(dtmin>(t-tb))
				{
 //                   cout<< "i="<<i<<"t="<<t <<"tb=" <<tb<<endl;

					dtmin=(t-tb);
				}

				if(dtmax<(t-tb))
				{
					dtmax=(t-tb);
				}
			}

			tb=t;


			number_lines++;

			if(number_lines==MAX)
			{
				cout<< "\n\n ** warning:  maximum line number reached. ** \n";
				break;
			}
		}
		j++;

	}
    inFile.clear();
    inFile.seekg(0,std::ios::beg);
    inFile.close();

	dt=(t-t1)/double(number_lines-1);

	sr=1./dt;

	cout<< " Time span: "<< t1 <<" to "<< t <<" sec " << endl;

	cout<< "\n dtmin= " << dtmin <<" sec";
	cout<< "\n dt   = " << dt    <<" sec";
	cout<< "\n dtmax= " << dtmax <<" sec \n";

	cout<< "\n srmin= " << 1./dtmax <<" samples/sec";
	cout<< "\n sr   = " << sr       <<" samples/sec";
	cout<< "\n srmax= " << 1./dtmin <<" samples/sec";

    last=number_lines;
    cout << "\n\n number of input points = " << last << endl;

	if( dtmin < 0.99 *dtmax)
	{

	   int nsr;

	   cout << "\n *** Time Step Warning ***   \n";
       cout << "\n dtmax = " << dtmax <<"   dtmin = " <<dtmin <<endl;
	   cout << "\n The time step must be constant.  \n";
       cout << "\n Enter new sample rate?  1=yes  2=no  \n";

       cin >> nsr;

	   if(nsr==1)
       {
			cout << "\n\n Assume constant sample rate. Enter sample rate (samples/sec) " << endl;
			cin >> sr;

			dt=1./sr;

			cout << "\n\n sr   = " << sr << " samples/sec";
			cout << "\n Nyquist Freq = " << sr/2. << " Hz" ;
			cout << "\n dt   = " << dt <<" sec	\n"<<endl;

	   }
	   else
	   {
			cout<<"\n Press any key to exit \n";
			cin.get();
			cin.get();
			exit(1);

	   }
	}

}
void line_test(void)
{
	nmax=0;
	long j=1;
    char string[300];

	while(1)
	{
        inFile.getline(string,numBytes);
        if(inFile.eof()){break;}

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

        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)
        {nmax=j;}

		j++;


	}
	//	         rewind(inFile);

    inFile.clear();
    inFile.seekg(0,std::ios::beg);
    inFile.close();

    cout<< "\n nmax = " << nmax <<endl;
}