//
// Butterworth filter, infinite impulse response
//
// cascade with refiltering option for phase correction
//

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

#define MAX 10000000

//************ complex number handling ********************
typedef struct 
{
    double r,i;
	
}doublecomplex;

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

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

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

int ik;
int iphase;

int kflag;

unsigned long i;
unsigned long ns;

unsigned long iband;
unsigned long iflag;
unsigned long j;
unsigned long k;
unsigned long l;





double nmaxp=-1000000;

double nminp= 1000000;
      
double a[4][4],b[4][5];
double alpha[20];

double dt;
double f;
double fl;
double fh;
double freq;

double maxp;

double minp;



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

double sr;
double ttime;
double ttime1;

float  y[MAX];
float yt[MAX];

doublecomplex s[20];





long numBytes=300;

char string[300];

void files(void);
void read1(void);
void read2(void);
void butterworth(void);
void normf(void);
void poles(void);
void alphac(void);
void lco(void);
void hco(void);
void dtrans(void);
void stab(void);
void apply(void);
void applymethod(void);
void stage1(void);
void stage2(void);
void printco(void);

void line_test(void);



int mflag=0;


int 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[1], argv[1] );



			mflag=1;

	}



   iflag=1;

   printf("\n  Program filter.cpp  ver 2.7 \n  July 23, 2010 \n  by Tom Irvine \n"); 
   printf("\n Butterworth 6th order with refiltering ");
   printf("\n for option phase correction. \n\n");

//
//   iband=1 for LP filter  
//         2 for HP filter
//

//****** initialize *******

   fl=0.;
   fh=0.;
   iband=0;

//  l = order

   l=6;

//****** open files *******

   files();

//****** read input data **************

   read1();

   read2();

//****** calculate coefficients *******

   if(iflag != 999 && iband !=3)
   {
       butterworth();
   }

   if(iflag < 900 )
   {
      if(iband == 1 || iband ==2)
      {
			applymethod();
      }

	  if(iband == 3)
      {
		 f=fh;
		 freq=f;

		 printf("\n Step 1");
		 iband=2;
         butterworth();
		 applymethod();

		 f=fl;
		 freq=f;

		 printf("\n\n Step 2");
		 iband=1;
		 butterworth();
		 applymethod();
      }


		double mu=0.;
		double ms=0.;
		double rms=0.;
		double std=0.;



		for(i=0; i<ns; i++)
		{
			fprintf(pFile[2]," %14.7e %14.4e \n", (dt*i)+ttime1, y[i]);

			mu+=y[i];

			ms+=pow(y[i],2.);


			if(y[i]>nmaxp)
			{

				nmaxp=y[i];
			}
			if(y[i]<nminp)
			{
				nminp=y[i];

			}

		} 

		mu/=double(ns);
		ms/=double(ns);

		rms=sqrt(ms); 

		std=sqrt(ms - pow(mu,2));


		printf("\n\n Output ttime History Statistics \n");


		printf("\n maximum = %8.4g   ",nmaxp);

		printf("\n minimum = %8.4g \n",nminp);


		printf("\n    mean = %12.4g ",mu);
		printf("\n std dev = %12.4g ",std);
		printf("\n     rms = %12.4g \n",rms);

		printf("\n\n ttime history output file: %s \n",filename[2]);
		printf("\n Please call the output file into your own plotting program.\n");
   }
   else
   {
       printf("\n\n  Abnormal termination.  No output file generated. ");
   }
   fclose( pFile[2] );

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


}
void applymethod(void)
{
	if(iphase==1)
	{
			apply();
			apply();  
	}
	else
	{
			apply();	
	}
}

void read2(void)
{

 /*  fscanf(pFile[0],"%i %lf %lf", &iband, &fl, &fh ); */       

      printf("\n Enter filter type");
	  printf("\n 1=lowpass  2=highpass  3=bandpass \n");
      scanf("\n %ld", &iband);

	  if((iband != 1) && (iband != 2) && (iband != 3))
      {
         printf("\n  iband error: iband = %i ",iband);
      }


	  f=0.;
	  fl=0.;
	  fh=0.;

      if(iband == 1)
      {
          printf("\n Enter lowpass frequency (Hz) ");
		  scanf("\n %lf", &fl );
          f=fl;
	  } 
      
      if(iband == 2)
      {
          printf("\n Enter highpass frequency (Hz) ");
		  scanf("\n %lf", &fh );      
          f=fh;
	  }
      if(iband == 3)
      {
		  printf("\n\n This bandpass filter is implemented as a"
			     "\n highpass filter and lowpass filter in series.\n");
    
		  printf("\n Enter lower frequency (Hz) ");
		  scanf("\n %lf", &fh ); 
		  printf("\n Enter upper frequency (Hz) ");
		  scanf("\n %lf", &fl );
          f=fh;
	  }
	  printf("\n\n Refiltering for phase correction?");
	  printf("\n 1=yes (appropriate for transient signal)");
	  printf("\n 2=no  (appropriate for steady-state signal) \n");
	  scanf("\n %d", &iphase );

      if(iphase != 1 && iphase !=2 ){iphase=1;}

	  freq=f;


	if( f > 0.5*sr)
	{

		printf("\n error: cutoff frequency must be < Nyquist frequency");    

		iflag=999;

	}


}
void read1(void)
{

	double dtmin=1.0e+99;
	double dtmax=-dtmin;
    double tb;

	printf("\n\n reading time history ");

	i=0;

	printf("\n");



	while( fgets(string,numBytes,pFile[1] )>0 )
	{

		kflag=0;



		line_test();

//		printf(" %ld %s  \n",i,string);



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

			if(i>0)
			{
				dt=ttime-tb;

				if(dt<dtmin)
				{
					dtmin=dt;
				}
				if(dt>dtmax)
				{
					dtmax=dt;
				}
			}

			tb=ttime;

			i++;



			if(i==MAX)
			{

				printf("\n Warning:  maximum number of data points reached. \n");

				break;

			}

		}

	}

	ns=i;

	double duration=(ttime-ttime1);

	fclose( pFile[1] );

	double ave=0.;
	double ms=0.;
	double rms=0.;


	maxp=-1.0+90;

	minp=-maxp;


	for(i=0; i <ns; i++)
	{ 
		ave+=y[i];

		ms+=pow(y[i],2);



		if(y[i]>maxp)

		{

			maxp=y[i];

		}

		if(y[i]<minp)

		{

			minp=y[i];

		}


	}

	ave/=double(ns);

	rms=sqrt(ms/double(ns));



	double std;



	std=sqrt(pow(rms,2)-pow(ave,2));


	printf("\n Input File Statistics \n");

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

	printf("\n maximum = %8.4g   ",maxp);

    printf("\n minimum = %8.4g \n",minp);


	printf("\n average = %8.4g   ",ave);

	printf("\n std dev = %8.4g   ",std);
	printf("\n     rms = %8.4g \n",rms);




	printf("\n start ttime=%9.5g ",ttime1);
	printf("\n   end ttime=%9.5g ",ttime);
	printf("\n   duration=%9.5g \n",duration);


	if( ns > 1 )
    {

	    dt=duration/float(ns-1);

        sr=1./dt;

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

		printf("\n srmin= %11.6g samples/sec	   ",1./dtmax);
		printf("\n sr   = %11.6g samples/sec	   ",sr);
		printf("\n srmax= %11.6g samples/sec	 \n",1./dtmin);

		printf("\n ttime step \n");
	
		printf("\n dtmin= %11.6g sec	   ",dtmin);
		printf("\n dt   = %11.6g sec	   ",dt);
		printf("\n dtmax= %11.6g sec	 \n",dtmax);

		printf("\n Nyquist Freq = %lf Hz \n",sr/2.);

    }
	else
    {
		 printf("\n  dt error ");

	}

	if( dtmin < 0.99 *dtmax)
	{

	   int nsr;

	   printf("\n *** ttime Step Warning ***   \n");

       printf("\n dtmax = %8.4g   dtmin = %8.4g \n",dtmax,dtmin);

	   printf("\n The ttime step must be constant. \n");


       printf("\n Enter new sample rate?  1=yes  2=no \n");

	   scanf("%d",&nsr);



	   if(nsr==1)

       {

			printf("\n\n Assume constant sample rate. Enter sample rate (samples/sec) \n");



			scanf("%lf",&sr);

			dt=1./sr;



			printf("\n\n sr   = %11.6g samples/sec	   ",sr);

			printf("\n Nyquist Freq = %lf Hz \n",sr/2.);

			printf("\n dt   = %11.6g sec	   \n",dt);

	   }

	   else

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

			getch();



			exit(1);

	   }
	   

	}

}

void files(void)
{



	if(mflag==0)

	{


		printf("\n Input file must have two columns: ttime(sec) and amplitude(units) ");
		printf("\n The format is free, but no header lines are allowed.\n");
		printf("\n Please enter the input filename. \n");



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

	}




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



	while(pFile[1] == NULL )
	{

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

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

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

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

	}

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


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

   printf( "\n\n Enter the output filename: \n");	              
   scanf("%s",filename[2]);
   pFile[2]=fopen(filename[2], "w");


   printf("\n");


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

   strcpy(filename[3], "s_trans.out");
   pFile[3]=fopen(filename[3], "w");

   strcpy(filename[4], "coeff.out");
   pFile[4]=fopen(filename[4], "w");

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

   for(i=2; i<=4; i++)
   {
      if(pFile[i] == NULL )
      {
         printf("%ld  Failed to open file: %s \n", i,filename[i]);
         fclose(pFile[i]);
      }
      else
      {
         printf("%ld  File: %s opened. \n", i,filename[i]);
      }
   }
}

void apply(void)
{
	float temp;

    if(iphase==1)
	{

		printf("\n begin reversal ");	

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

//
//  cascade stage 1
//

    printf("\n stage 1");
    ik=1;
	stage1();
	
//
//  cascade stage 2
//

    printf("   stage 2");
    ik=2;
	stage2();
	
//
//  cascade stage 3
//
	 printf("  stage 3");
	 ik=3;
	 stage1();

	 for(i=0; i<ns; i++)
	 {
	      y[i]=yt[i];
	 }
}

void stage1(void)
{
	for(i=0; i< MAX; i++)
    {
		yt[i]=float(0.);
    }
               
    yt[0]=float( b[ik][0]*y[0] );       
    yt[1]=float( b[ik][0]*y[1]+b[ik][1]*y[0]-a[ik][1]*yt[0] );
	  
    for(i=2; i< ns; i++)
    {
         yt[i]=float( b[ik][0]*y[i] + b[ik][1]* y[i-1] +b[ik][2]* y[i-2]
                             - a[ik][1]*yt[i-1] -a[ik][2]*yt[i-2] );         
	}
}

void stage2(void)
{
	for(i=0; i< MAX; i++)
    {
		y[i]=float(0.);
    }
               
    y[0]=float( b[ik][0]*yt[0] );       
    y[1]=float( b[ik][0]*yt[1]+b[ik][1]*yt[0]-a[ik][1]*y[0] );
	  
    for(i=2; i< ns; i++)
    {
         y[i]=float( b[ik][0]*yt[i] + b[ik][1]*yt[i-1] +b[ik][2]*yt[i-2]
                             - a[ik][1]* y[i-1] -a[ik][2]* y[i-2] );         
	}
}

void butterworth(void)
{

// 6th order
    	
	double alpha[20];
	
	for(j=0; j<=3; j++)
    {
		for( k=0; k<=3; k++)
        {
			a[j][k]=0.;
		}
	}

	for(j=0; j<=3; j++)
    {
		for( k=0; k<=4; k++)
        {
			b[j][k]=0.;
		}
	}

	for(j=0; j<=19; j++)
    {
		alpha[j]=0.;
		s[j].r=0.;
		s[j].i=0.;
    }

//*** normalize the frequency ***

	normf();

//*** solve for the poles *******

	poles();

//*** solve for alpha values ****

    alphac();

//*** solve for filter coefficients **

    if( iband == 1 )
    {
		lco();
    }
	else
	{
		hco();
    }

//*** plot digital transfer function **

	dtrans();

//*** check stability ****************

	stab();

}

void normf(void)
{
	   double targ=0.;

       targ=pi*f*dt;
   //    printf("\n  arg = %lf radians", targ );
       om=tan(targ);  
  //     printf("\n   om = %lf ", om );   
}

void poles(void)
{
        doublecomplex h;
		doublecomplex theta;
		doublecomplex temp;

		double a1, a2, a3;
		double arg=0.;
        double xr, xi;

//		printf("\n  calculate poles ");
		
		for( k=0; k<= (2*l-1); k++)
        {
			arg = (2.*(k+1) +l-1)*pi/(2.*l);
			
			s[k].r= cos(arg);
            s[k].i= sin(arg);

//			printf("\n %i %12.6g %12.6g ", k, s[k].r, s[k].i );
        }

        for( i=0; i<=200; i++)
        {
           arg = i/40.;
		   
		   h.r =     - s[0].r;
		   h.i = arg - s[0].i;

           for( j=1; j<=(l-1); j++)		   
           {
              theta.r =     - s[j].r;
			  theta.i = arg - s[j].i;
			   
			  COMPLEXMUL(h.r,h.i,theta.r,theta.i,temp.r,temp.i)

			  h.r = temp.r;
			  h.i = temp.i; 
		   }
          
		   COMPLEXRECIP(h.r, h.i, xr, xi)

           h.r=xr;
		   h.i=xi;

		   a1 = freq*arg;
		   
		   COMPLEXMAG( h.r, h.i, a2 )
           
		   a3 = pow(a2,2);

		   fprintf(pFile[3]," %lf %lf %lf \n", a1, a2, a3);
      
        }
        fclose( pFile[3]);
}

void alphac(void)
{
		printf("\n");
		for(j=0; j <= (l-1); j++)
        {
			alpha[j] = 2.*s[j].r;
 //           printf("\n alpha[%ld]= %12.5g ",j,alpha[j]);
		}
}

void lco(void)
{
		double den;
		double om2;

		printf("\n");

		om2=pow(om,2);

		for(k=1; k<=l/2; k++)
        {
		    den = om2-alpha[k-1]*om+1.;
		
			a[k][0]=0.;
			a[k][1]=2.*(om2 -1.)/den;
            a[k][2]=( om2 +alpha[k-1]*om+ 1.)/den;

			  b[k][0]=om2/den;
              b[k][1]=2.*b[k][0];
              b[k][2]=b[k][0];
		
		printco();
		}
	
}	 


void hco(void)
{
	    double den;
		double om2;

		printf("\n");

		om2=pow(om,2);

		for(k=1; k<=l/2; k++)
        {
			den = 1. -alpha[k-1]*om +om2;
		
			a[k][0]=0.;
			a[k][1]=2.*(-1.+ om2)/den;
            a[k][2]=( 1.+alpha[k-1]*om+ om2)/den;

			  b[k][0]= 1./den;
              b[k][1]=-2.*b[k][0];
              b[k][2]=    b[k][0];
		}
		
}

void printco()
{
		printf("\n filter coefficients");

		printf("\n a[%i][1]=%10.5g  a[%i][2]=%10.5g",k,a[k][1],k,a[k][2]);
		printf("\n b[%i][0]=%10.5g  b[%i][1]=%10.5g  b[%i][2]=%10.5g",
			                                        k,b[k][0],k,b[k][1],k,b[k][2]);
		printf("\n");

		fprintf(pFile[4],"\n %i ",k);
		fprintf(pFile[4],"\n a[%i][1]=%10.5g  a[%i][2]=%10.5g",k,a[k][1],k,a[k][2]);
		fprintf(pFile[4],"\n b[%i][0]=%10.5g  b[%i][1]=%10.5g  b[%i][2]=%10.5g",
			                                        k,b[k][0],k,b[k][1],k,b[k][2]);
		fprintf(pFile[4],"\n");

		fclose( pFile[4] );
}

void dtrans(void)
{
         double fd;
		 double omd;
		 double zden;
		 double zz;
		 double znzd;

         doublecomplex z;
		 doublecomplex num[4];
		 doublecomplex den[4];
		 doublecomplex st;

//  (sr/2.) = Nyquist frequency

         long nnn = long(sr);

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

	         for(k=0; k<= 3; k++)
             {
			     num[k].r=0.;
			     num[k].i=0.;
			     den[k].r=0.;
			     den[k].i=0.;
             }

             fd=j*.5*dt;
			 omd=2.*pi*fd;
		
			 z.r=cos(omd);
			 z.i=-sin(omd);

             znzd=1.;

			 for(k=1; k<= 3; k++)
             {
				 st.r=0.;
				 st.i=0.;

				 COMPLEXMUL(z.r,z.i,z.r,z.i,st.r,st.i)
                
                 num[k].r=b[k][0] +b[k][1]*z.r  +b[k][2]*st.r;
				 num[k].i=        +b[k][1]*z.i  +b[k][2]*st.i;
				 
                 den[k].r=1.      +a[k][1]*z.r  +a[k][2]*st.r;
				 den[k].i=        +a[k][1]*z.i  +a[k][2]*st.i;
							       
                 zden=0.;

				 COMPLEXMAG(den[k].r,den[k].i,zden)
				 
				 if(zden < 0. ){zden=-zden;}

			     if( zden < 1.0e-14)
				 {
                     printf("\n error ");
                     printf("\n %14.5g ",zz);
                     iflag=908;                 
                 }		
            
			     zz=0.;
			     st.r=0.;
			     st.i=0.;

			     COMPLEXDIV(num[k].r,num[k].i,den[k].r,den[k].i,st.r,st.i)
                 COMPLEXMAG(st.r,st.i,zz)

	             znzd*=zz;
             } 
			 fprintf(pFile[5]," %12.4g %12.4g %12.4g \n",fd/dt,znzd,pow(znzd,2));
		 }
		 fclose( pFile[5] );
            

}
void stab(void)
{
		double a1;
		double d1, d2, d3;
        double dlit;

		double at1, at2;
		double als=0.5e-06;
		double h2;

		als*=6.;

		printf("\n stability reference threshold= %14.7e ", als);
    
	    for( i=1; i<= l/2; i++)
        {
			at1= -a[i][1];
			at2= -a[i][2];

		//		printf("\n\n stability coordinates: (%12.7g, %14.7g) ",at1,at2);
        
		    h2=at2;
 
            a1=h2-1.;
            d3=at1-a1;
         
            a1=1.-h2;
            d2=a1-at1;
            d1=at2+1.;
		
          //  printf("\n d1=%14.5g  d2=%14.5g  d3=%14.5g",d1,d2,d3);

            dlit=d1;

            if(dlit > d2){dlit=d2;}
            if(dlit > d3){dlit=d3;}

			printf("\n\n stage %ld     dlit= %14.5g ",i, dlit);

            if(dlit > als)
			{
				printf("\n good stability"); 
            }
			if( (dlit < als) && (dlit > 0.))
		    {
                printf("\n marginally unstable ");
            }
		    if(dlit < 0.)
		    {  
			    printf("\n unstable ");
		        iflag=905;
		    }
		}
	    printf("\n");
}

void line_test()

{



//		printf("j=%ld.  %s\n",j,string);



		if( strrchr(string,'A' ) != NULL ){kflag=1;}

		if( strrchr(string,'B' ) != NULL ){kflag=2;}

		if( strrchr(string,'C' ) != NULL ){kflag=3;}

		if( strrchr(string,'D' ) != NULL ){kflag=4;}

		if( strrchr(string,'E' ) != NULL && strstr(string,"E+")==NULL && strstr(string,"E-")==NULL ){kflag=5;}

		if( strrchr(string,'F' ) != NULL ){kflag=6;}

		if( strrchr(string,'G' ) != NULL ){kflag=7;}

		if( strrchr(string,'H' ) != NULL ){kflag=8;}

		if( strrchr(string,'I' ) != NULL ){kflag=9;}

		if( strrchr(string,'J' ) != NULL ){kflag=10;}

		if( strrchr(string,'K' ) != NULL ){kflag=11;}

		if( strrchr(string,'L' ) != NULL ){kflag=12;}

		if( strrchr(string,'M' ) != NULL ){kflag=13;}

		if( strrchr(string,'N' ) != NULL ){kflag=14;}

		if( strrchr(string,'O' ) != NULL ){kflag=15;}

		if( strrchr(string,'P' ) != NULL ){kflag=16;}

		if( strrchr(string,'Q' ) != NULL ){kflag=17;}

		if( strrchr(string,'R' ) != NULL ){kflag=18;}

		if( strrchr(string,'S' ) != NULL ){kflag=19;}

		if( strrchr(string,'T' ) != NULL ){kflag=20;}

		if( strrchr(string,'U' ) != NULL ){kflag=21;}

		if( strrchr(string,'V' ) != NULL ){kflag=22;}

		if( strrchr(string,'W' ) != NULL ){kflag=23;}

		if( strrchr(string,'X' ) != NULL ){kflag=24;}

		if( strrchr(string,'Y' ) != NULL ){kflag=25;}

		if( strrchr(string,'Z' ) != NULL ){kflag=26;}



//



		if( strrchr(string,'a' ) != NULL ){kflag=27;}

		if( strrchr(string,'b' ) != NULL ){kflag=28;}

		if( strrchr(string,'c' ) != NULL ){kflag=29;}

		if( strrchr(string,'d' ) != NULL ){kflag=30;}

		if( strrchr(string,'e' ) != NULL  && strstr(string,"e+")==NULL && strstr(string,"e-")==NULL ){kflag=31;}

		if( strrchr(string,'f' ) != NULL ){kflag=32;}

		if( strrchr(string,'g' ) != NULL ){kflag=33;}

		if( strrchr(string,'h' ) != NULL ){kflag=34;}

		if( strrchr(string,'i' ) != NULL ){kflag=35;}

		if( strrchr(string,'j' ) != NULL ){kflag=36;}

		if( strrchr(string,'k' ) != NULL ){kflag=37;}

		if( strrchr(string,'l' ) != NULL ){kflag=38;}

		if( strrchr(string,'m' ) != NULL ){kflag=39;}

		if( strrchr(string,'n' ) != NULL ){kflag=40;}

		if( strrchr(string,'o' ) != NULL ){kflag=41;}

		if( strrchr(string,'p' ) != NULL ){kflag=42;}

		if( strrchr(string,'q' ) != NULL ){kflag=43;}

		if( strrchr(string,'r' ) != NULL ){kflag=44;}

		if( strrchr(string,'s' ) != NULL ){kflag=45;}

		if( strrchr(string,'t' ) != NULL ){kflag=46;}

		if( strrchr(string,'u' ) != NULL ){kflag=47;}

		if( strrchr(string,'v' ) != NULL ){kflag=48;}

		if( strrchr(string,'w' ) != NULL ){kflag=49;}

		if( strrchr(string,'x' ) != NULL ){kflag=50;}

		if( strrchr(string,'y' ) != NULL ){kflag=51;}

		if( strrchr(string,'z' ) != NULL ){kflag=52;}



//



		if( strrchr(string,'!' ) != NULL ){kflag=53;}

		if( strrchr(string,'@' ) != NULL ){kflag=54;}

		if( strrchr(string,'#' ) != NULL ){kflag=55;}

		if( strrchr(string,'$' ) != NULL ){kflag=56;}

		if( strrchr(string,'%' ) != NULL ){kflag=57;}

		if( strrchr(string,'^' ) != NULL ){kflag=58;}

		if( strrchr(string,'&' ) != NULL ){kflag=59;}

		if( strrchr(string,'*' ) != NULL ){kflag=60;}

		if( strrchr(string,'/' ) != NULL ){kflag=61;}

		if( strrchr(string,'>' ) != NULL ){kflag=62;}

		if( strrchr(string,'<' ) != NULL ){kflag=63;}

		if( strrchr(string,'?' ) != NULL ){kflag=64;}

		if( strrchr(string,'=' ) != NULL ){kflag=65;}

		if( strrchr(string,';' ) != NULL ){kflag=66;}



		if( 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){ kflag=67;}

/*

		if(kflag==67)

		{



		if( strrchr(string,'1')==NULL)printf(" no 1 \n");

		if( strrchr(string,'2')==NULL)printf(" no 2 \n");

		if( strrchr(string,'3')==NULL)printf(" no 3 \n");

		if( strrchr(string,'4')==NULL)printf(" no 4 \n");

		if( strrchr(string,'5')==NULL)printf(" no 5 \n");

		if( strrchr(string,'6')==NULL)printf(" no 6 \n");

		if( strrchr(string,'7')==NULL)printf(" no 7 \n");

		if( strrchr(string,'8')==NULL)printf(" no 8 \n");

		if( strrchr(string,'9')==NULL)printf(" no 9 \n");

		if( strrchr(string,'0')==NULL)printf(" no 0 \n");



		}

*/		

}