/*

      spl.cpp

      by 
      Tom Irvine
      Email:  tomirvine@aol.com

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

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

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

  Two input files are required.  One output file is generated.

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

  Input Files:


     The input file is a free-format ASCII text file with two columns:

     Time(sec) and Pressure(units)

     
     The delta time between adjacent data points must be consistent
     throughout the file, otherwise the Fourier transform will 
     be in error.  


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

  Output File:


     This is an ASCII text file with two columns:  

         Freq (Hz)   Power Spectral Density (eu^2/Hz)
      
     This file can be called into a graphics program, such as
     Excel, Axum, Matlab, EasyPlot, etc.

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

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

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

#define MAX   131072
#define MAXP1 65536

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

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

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

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

floatcomplex y[MAX];
floatcomplex ww[MAXP1];
  
float hold[MAXP1];
  
double ae, ave, aw[200],bw[200],cw[200];
double df, dt, dur;
double freq,fband[200],fc[200];
double grms;
double kurt; 

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

double reference;
double sr, std, sum[100];
double time, time1, t, tt;

double mag[MAXP1];

unsigned long i, j, k, n, ns;
unsigned long ich;

unsigned long m, nave;

unsigned long nstages;
unsigned long nsegment;
unsigned long nhs;

int iflag, imr, ihw, iwn, iscale;


//**** function prototypes ****

void files(void);
void read2(void);
void read3(void);
void bitrev(void);
void fourier(void);
void meanremove(void);
void hanning(void);
void ctfourier(void);
void psd(void);
void holddat(void);
void advise(void);
void weights(void);

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

		mflag=1;

	}

	iflag=0;

	printf("\n\n\n\n");
	printf("\n spl   ver 1.6   February 3, 2009\n");
	printf("\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com \n");

	printf("\n This program calculates a sound pressure level from a pressure time history.");
	printf("\n The medium is assumed to be air.\n");
    
	for(i=0; i< MAX; i++)
    {
		 y[i].r=float(0.);
         y[i].i=float(0.); 
	}
	for(i=0; i< MAXP1; i++)
    {
		 ww[i].r=float(0.);
         ww[i].i=float(0.);
		 
		 mag[i]=0.;
		 hold[i]=float(0.);
	}

//***** open files ****************
    
	files();

	advise();
    rewind(pFile[1]);
    
//***** read segment length *******
	
//	fscanf(pFile[0],"%ld", &nsegment );
    

	iflag=1;
	while( 1>0)
	{
//	   printf("\n Enter samples per segment value from third column. \n ");
//	   scanf("%ld",&nsegment );

       if(nsegment> MAXP1){nsegment=MAXP1;}
	   n=nsegment;
       

       nhs = nsegment/2;

	   ich=0;
	   for(i=1; i<=16; i++)
       {
          if(nsegment == unsigned long (pow(2,i))	 )
	      { 
		      ich=1;
              nstages=i;	   
	      }
	   } 
       if( ich == 0)
       {
	     	printf("\n Error: illegal n value.  n must be 2^i, where i is an integer.");    
	 	    if(n > MAXP1 )
		    { 
			   printf("\n maximum allowed value is %ld ", MAXP1);
		    }
		    iflag=999;
	   }
	   if( iflag == 1){break;}
	}
	
   
//***  read mean removal and window parameters

//  fscanf(pFile[0],"%d", &imr );
//  fscanf(pFile[0],"%d", &ihw );

//    printf("\n\n Mean removal? ");
//    printf("\n 1=yes 2=no \n");

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

	imr=1;

	printf("\n\n Enter window type ");
    printf("\n 1=rectangular  2=Hanning \n");

	scanf("%d", &ihw );

    printf("\n");

	printf("\n\n Apply Weighting Network? ");
    printf("\n     1=No (physics & engineering) ");
    printf("\n     2=Yes (human hearing) \n");

	scanf("%d", &iwn );

    if(iwn == 2)
	{
	   printf("\n\n Select Scale ");
       printf("\n     1=A-weight   2=B-weight   3=C-weight \n");

	   scanf("%d", &iscale );
	}

    printf("\n\n processing data....\n");


//  fclose( pFile[0] );

    int irun=0;
	
	nave=0;
	while( irun==0)
	{
//*** read time history

	   if( nave == 0)
	   {
	      read2();
		  holddat();
	   }
	   else
	   {
		  read3();
		  holddat();

		  if(iflag==905){break;}
	   }
    
	   if( iflag < 900 )
       {

//*** perform mean removal if selected.

	      if(imr==1){ meanremove();}

//*** perform hanning window if selected.

	      if(ihw==2){ hanning(); }

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

	      bitrev();

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

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

	}

	fclose( pFile[1] );

    psd();

	printf("\n\n\n Output file:  %s \n",filename[2]);
	printf("\n Dimensions are: freq(Hz) & SPL (dB) ");

    printf("\n\n zero dB Reference = 20 micro Pascal \n");

	printf("\n Please call the output file into your own plotting program. \n");

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

	exit(1);
}
void ctfourier(void)
{
	int i1;

	unsigned long ig;
	
	int iy=1;
	int jj;
	int jjj;
    int l;
	int lg;
	int lk;
	int lwd2;

	unsigned long m;

	int nworld;
	int wd2;
	int world=4;

	double arg;

	double z=0.;

    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=tp*double(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(nsegment);
	   m=nsegment;
    
//     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);

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


             COMPLEXMAG2( y[ig].r,y[ig].i,z)
			 mag[ig]+=z;

//			 printf("\n mag[ig] = %10.3e ", mag[ig] );
       }

//     printf("\n\n sum real=%lf   sum imag=%lf",asum.r,asum.i);
//       printf("\n\n output data written to:  a.out \n");
//	   fclose( pFile[2] );
	   
}
void read2(void)
{
    double sum, sum2, sum4, tn;

	printf("\n\n reading time history ");
	
	if( ( fscanf(pFile[1]," %lf %f ",&time1, &y[0].r) ) <= 0 )
    {
		printf("\n fscanf error ");
    }
	sum=     y[0].r;
    sum2=pow(y[0].r,2);
    sum4=pow(y[0].r,4);

    i=1;
	while ( fscanf(pFile[1]," %lf %f ",&tn, &y[i].r) > 0 )
    {
     	 sum+=     y[i].r;
         sum2+=pow(y[i].r,2);
         sum4+=pow(y[i].r,4);

		 i++;
		 if(i == nsegment)
		 { 
			 time=tn;
			 break;
		 }

	}                               
	ns=i;

//	printf("\n\n  %ld samples read out of %ld specified\n",ns, nsegment);

	if(ns < nsegment)
    {
		iflag=905;
	}

	if(iflag != 905)
    {
	if( ns > 1)
    {
	    dt=(time-time1)/float(ns-1);
        sr=1./dt;

		dur = (time-time1);
        df=1./dur;
	
        ave=sum/i;
        grms=sqrt(sum2/i);
        std = sqrt( pow(grms,2) - pow(ave,2) );
        kurt=sum4/(i*pow(std,4));

        printf("\n samples= %ld ", i);
        printf("\n dt = %10.5g sec  sr= %12.6e samples/sec", dt, sr);
        printf("\n\n time span= %12.6g sec to %12.6g sec", time1, time);
   
        printf("\n\n  ave= %9.4g     std dev = %9.4g    rms= %9.4g ",
	                            ave, std, grms);
        printf("\n kurtosis= %9.4g", kurt );
        printf("\n");

		printf("\n delta f= %14.4g ",df);

	}	
	else
    {
		 printf("\n  dt error ");
	}
	}
//	fclose( pFile[1] );
}

void read3(void)
{
    double sum, sum2, sum4, tn;

	sum=0.;
	sum2=0.;
	sum4=0.;
	
	printf("\n\n reading time history ");

	for(i=0; i< MAX; i++)
    {
         y[i].i=float(0.); 
	}
		
	for(i=0; i < nhs; i++)
	{
		y[i].r = hold[i];
	}
	for( i=nhs; i < nsegment; i++)
	{
		y[i].r = float(0.);
	}
	
    i=0;
	while ( fscanf(pFile[1]," %lf %f ",&tn, &y[i+nhs].r) > 0 )
    {
  
		 i++;
		 if(i+nhs == nsegment)
		 { 
//			 time=tn;
			 break;
		 }

	}                               
	ns=i;

//	printf("\n\n  %ld samples read out of %ld specified\n",ns, nsegment);

	if(ns < nhs)
    {
		iflag=905;
//		printf("\n  error:  insufficient samples in data ");
	    printf("\n   ns=%ld  nhs=%ld ",ns,nhs);
	}

	if(iflag != 905)
    {
	if( ns > 1)
    {
//	    dt=(time-time1)/float(ns-1);
//        sr=1./dt;

//		dur = (time-time1);
//        df=1./dur;

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

		     sum+=    y[i].r;
            sum2+=pow(y[i].r,2);
            sum4+=pow(y[i].r,4);
		}
	
        ave=sum/i;
        grms=sqrt(sum2/i);
        std = sqrt( pow(grms,2) - pow(ave,2) );
        kurt=sum4/(i*pow(std,4));

 //       printf("\n samples= %ld ", i);
 //       printf("\n dt = %10.5g sec  sr= %12.6e samples/sec", dt, sr);
 //       printf("\n\n time span= %12.6g sec to %12.6g sec", time1, time);
   
        printf("\n\n  ave= %9.4g     std dev = %9.4g    rms= %9.4g ",
	                            ave, std, grms);
        printf("\n kurtosis= %9.4g", kurt );
        printf("\n");

//		printf("\n delta f= %14.4g ",df);

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

}

void holddat(void)
{
	for(i=0; i<nsegment; i++)
	{
		hold[i]=y[i+nhs].r;
	}
}


void psd(void)
{
//    printf("\n nave=%ld");

	double rms=0.;

// mag has already undergone complex conjugate multiplication

	m=nsegment;

//  The zero case is treated separately.  Neither two-sided to one-sided nor
//  rms conversion is needed at this frequency.
	
   	mag[0]*=(1./df);
    mag[0]*=(1./double(nave));

	freq=0.;

//	fprintf(pFile[2],"%10.4e %10.4e \n",freq,mag[0]);
  
    for( i=1; i < m/2 ; i++)
	{
//		printf("\n *** mag[i]= %10.3e", mag[i] );
		
	    freq=i*df;

	 	mag[i]*=(1./df);

		mag[i]*=(1./double(nave));

		mag[i]*=4.;        // two-sided to one-sided conversion.
		mag[i]*=0.5;	   // convert from peak^2 to rms^2

        rms+=mag[i];

//		fprintf(pFile[2],"%10.4e %10.4e \n",freq,mag[i]);
	}
	

    weights();

	for(i=0; i<=33; i++)
	{
		sum[i]=0.;
//		fc[i]=10.*pow(2.,double(i)/3.);
        fband[i]=fc[i]/(pow(2.,(1./6.)));

//		printf(" fc[%ld]=%8.4g Hz\n",i,fc[i]);
	}

    for(i=1; i<m/2; i++)
    {
		freq=i*df;
        mag[i]*=df;


		for(j=0; j<=32; j++)
		{
			if( freq >= fband[j] && freq < fband[j+1] )
			{
			   sum[j]+=mag[i];
			   break;
			}
        }
	}

	printf("\n\n  fc(Hz)  SPL(dB) \n");

	for(i=0; i<= 31; i++)
	{
		if( sum[i] >= 1.0e-50)
		{
			sum[i]=20.*log10( sqrt(sum[i])/reference );
        }

		if(iwn==2)
		{
			if(iscale==1)
			{
				sum[i]+=aw[i];
			}
            if(iscale==2)
			{
				sum[i]+=bw[i];
			}
			if(iscale==3)
			{
				sum[i]+=cw[i];
			}
		}
		if( sum[i] >= 1.0e-50)
		{
			fprintf(pFile[2],"%10.4e %10.4e \n",fc[i],sum[i]);
        }
		printf("%10.4g \t %10.4g \n",fc[i],sum[i]);
	}
	printf("\n\n");

	fclose( pFile[2] );


	printf("\n m= %ld", m );
	printf("\n df = %10.3e ",df);
	printf("\n nave = %ld ",nave);

    rms=sqrt(rms*df);

	double dB = 20.*log10(rms/reference);

	printf("\n\n overall rms = %10.4g dB \n",dB);

	double psi_rms = (2.9e-09)*pow(10,(dB/20.));

    printf("\n  pressure in air = %8.4g psi rms  ",psi_rms);
    printf("\n                  = %8.4g Pa  rms  \n",psi_rms*(6891.2));
}


void meanremove()
{
	printf("\n mean removal");
	for(i=0; i < nsegment; i++)
    {
	    y[i].r-=float(ave);
	}
}

void hanning()
{
	ae=sqrt(8./3.);

	printf("\n Hanning window ");
	for(i=0; i < nsegment; i++)
    {
	    y[i].r*=float( ae*pow(sin( (i*pi/nsegment) ),2) );
	}
}

void files(void)
{
  // strcpy(filename[0], "fft.in");
  // pFile[0]=fopen(filename[0], "rb");

	if(mflag==0)
	{
		printf("\n The input file must have two columns:  time(sec) & pressure ");

		printf("\n Enter input filename\n");
		scanf("%s",&filename[1]);
	}
	
    pFile[1]=fopen(filename[1], "rb");

//   strcpy(filename[2], "a.out");
//   pFile[2]=fopen(filename[2], "w");


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


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

   long ipress;
   
   printf("\n\n Select input data pressure unit. ");
   printf("\n 1=psi   2=Pascal   3= micro Pascal  \n\n");
   scanf("%ld",&ipress);

   if(ipress==1)
   {
	   reference = 2.9e-09;
   }
   if(ipress==2)
   {
	   reference = 20.e-06;
   }
   if(ipress==3)
   {
	   reference = 20.;
   }

} 

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

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

       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;
		  }
          
          kk=0;                      
          ia=0;

          ilk=unsigned 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=ia + md2*2;
          
          }
          md2=(md2/2);   
		  
        }
 }

void advise(void)
{
	unsigned long ith=1;

	unsigned long idd;
	unsigned long nf, ns;


	double dt;
	double dur,durmin;

	double rinn;
	double sr;

	double t, t1,a;
   

    printf("\n\n reading input file.....\n\n");

	fscanf(pFile[1],"%lf %lf",&t1,&a);

    for(i=2; i < 99999999; i++)
	{
		if ( fscanf(pFile[1],"%lf %lf",&t,&a) <= 0){break;}
	}
 
       dur=t-t1;
       ns=i-1;

       dt=dur/ns;
       
	   sr=1./dt;

       if(ns>0){dt=dur/ns;}

       printf("\n time %lf to %lf sec",t1,t);
   
       printf("\n dur = %11.4f sec\tns = %ld samples", dur, ns );        

       rinn=log(ns)/log(2);
       idd=unsigned long (rinn);
       nf =unsigned long ( pow(2,idd) );

       nsegment=nf;

//	   seg=1;
//	   while( 1 > 0 )
//	   {
//	      if(nf > MAXP1)
//	      {
//              nf/=2;
//			  seg*=2;
//	      }
//	      else
//	      {
//		      break;
//	      }
//	   }
       
    //   printf("\n rinn=%lf idd=%ld nf=%ld",rinn,idd,nf);

//       durmax=dur; 
       durmin=nf*dt;    
       
       printf("\n dt = %11.6f sec\tsr = %11.3f samples/sec ",dt,sr);

       df=1./durmin;

//      idof=unsigned long(2.*df*durmin);
//	   se=1./sqrt(df*dur);
       
//       ed=durmin;
//       ii=1;
      
//	   printf("\n");
//       printf("\n Choice\t    DF\t      Samples\tSegments  Dof\t Duration\tStd");
//       printf("\n       \t    (Hz)      per seg  \t\t\t\t        Error");

//       printf("\n %ld\t%10.3f\t%ld\t   %ld\t   %ld\t%10.4f  %10.4f",ii,df,nf,seg,idof,ed,se);

//      while( ith ==1)
//      {
//	   df*=2.;
      
//       se=1./sqrt(df*dur);       
//	   nf=unsigned long(nf/2);
       
//       seg=unsigned long(ns/nf);
//	   ed=seg*nf*dt;
//       idof=unsigned long(2.*seg);

//	   if(nf > 4)
//	   {
//	   ii++;
//	   printf("\n %ld\t%10.3f\t%ld\t   %ld\t   %ld\t%10.4f  %10.4f",ii,df,nf,seg,idof,ed,se);
//	   }
//	   else
//	   {
//		   break;
//	   }
//	   }
//	   printf("\n");

}
void weights(void)
{
	fc[0]=10.;
    aw[0]=-70.4;
	bw[0]=-38.2;
	cw[0]=-14.3;

	fc[1]=12.5;
    aw[1]=-63.4;
	bw[1]=-33.2;
	cw[1]=-11.2;

	fc[2]=16.;
    aw[2]=-56.7;
	bw[2]=-28.5;
	cw[2]=-8.5;

	fc[3]=20.;
    aw[3]=-50.5;
	bw[3]=-24.2;
	cw[3]=-6.2;

	fc[4]=25.;
    aw[4]=-44.7;
	bw[4]=-20.4;
	cw[4]=-4.4;

	fc[5]=31.5;
    aw[5]=-39.4;
	bw[5]=-17.1;
	cw[5]=-3.0;

	fc[6]=40.;
    aw[6]=-34.6;
	bw[6]=-14.2;
	cw[6]=-2.0;

	fc[7]=50.;
    aw[7]=-30.2;
	bw[7]=-11.6;
	cw[7]=-1.3;

    fc[8]=63.;
    aw[8]=-26.2;
	bw[8]=-9.3;
	cw[8]=-0.8;

	fc[9]=80.;
    aw[9]=-22.5;
	bw[9]=-7.4;
	cw[9]=-0.5;

	fc[10]=100.;
    aw[10]=-19.1;
	bw[10]=-5.6;
	cw[10]=-0.3;

	fc[11]=125.;
    aw[11]=-16.1;
	bw[11]=-4.4;
	cw[11]=-0.2;

	fc[12]=160.;
    aw[12]=-13.4;
	bw[12]=-3.0;
	cw[12]=-0.1;

	fc[13]=200.;
    aw[13]=-10.9;
	bw[13]=-2.0;
	cw[13]=0.;

	fc[14]=250.;
    aw[14]=-8.6;
	bw[14]=-1.3;
	cw[14]=0.;

	fc[15]=315.;
    aw[15]=-6.6;
	bw[15]=-0.8;
	cw[15]=0.;

    fc[16]=400.;
    aw[16]=-4.8;
	bw[16]=-0.5;
	cw[16]=0.;

	fc[17]=500.;
    aw[17]=-3.2;
	bw[17]=-0.3;
	cw[17]=0.;

	fc[18]=630.;
    aw[18]=-1.9;
	bw[18]=-0.1;
	cw[18]=0.;

	fc[19]=800.;
    aw[19]=-0.8;
	bw[19]=0.;
	cw[19]=0.;

	fc[20]=1000.;
    aw[20]=0.;
	bw[20]=0.;
	cw[20]=0.;

	fc[21]=1250.;
    aw[21]=0.6;
	bw[21]=0.;
	cw[21]=0.;

	fc[22]=1600.;
    aw[22]=1.0;
	bw[22]=0.;
	cw[22]=-0.1;

	fc[23]=2000.;
    aw[23]=1.2;
	bw[23]=-0.1;
	cw[23]=-0.2;

    fc[24]=2500.;
    aw[24]=1.3;
	bw[24]=-0.2;
	cw[24]=-0.3;

	fc[25]=3150.;
    aw[25]=1.2;
	bw[25]=-0.4;
	cw[25]=-0.5;

	fc[26]=4000.;
    aw[26]=1.0;
	bw[26]=-0.7;
	cw[26]=-0.8;

	fc[27]=5000.;
    aw[27]=0.5;
	bw[27]=-1.2;
	cw[27]=-1.3;

	fc[28]=6300.;
    aw[28]=-0.1;
	bw[28]=-1.9;
	cw[28]=-2.0;

	fc[29]=8000.;
    aw[29]=-1.1;
	bw[29]=-2.9;
	cw[29]=-3.0;

	fc[30]=10000.;
    aw[30]=-2.5;
	bw[30]=-4.3;
	cw[30]=-4.4;

	fc[31]=12500.;
    aw[31]=-4.3;
	bw[31]=-6.1;
	cw[31]=-6.2;

	fc[32]=16000.;
    aw[32]=-6.6;
	bw[32]=-8.4;
	cw[32]=-4.4;

	fc[33]=20000.;
    aw[33]=-9.3;
	bw[33]=-11.1;
	cw[33]=-11.2;

	fc[34]=25000.;

	fc[35]=31500.;
}