//      intlog, ver 1.1

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

#define MAX1 20000
#define MAX2 100000


void files(void);


double fr[MAX1],r[MAX1],s[MAX1],psd[MAX2],f[MAX2];

double a,az,b,bw,fa,ff,fz;

long i,j,n,nt;

int mflag=0;


char filename[5][FILENAME_MAX];

FILE *pFile[5];

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 intlog.cpp, ver 1.3   October 21, 2008");
	   printf("\n\n by Tom Irvine ");
	   printf("\n Email:  tomirvine@aol.com ");

	   printf("\n\n This program performs a logarthmic interpolation on a frequency ");
	   printf("\n domain function, such as a transmissibilty or power spectral density");
	   printf("\n function.\n\n");

	   files();
       
       printf("\n Please enter frequency bandwidth (delta f)\n");
       scanf("%lf",&bw);

	   printf("\n Please enter starting frequency for output data.\n");
       scanf("%lf",&fa);

	   n=-1;

	   for(i=0; i< MAX2; i++)
	   {
		  if( fscanf(pFile[0],"%lf %lf",&a,&b) <= 0 ){break;}

		  if(a>= 0.000001)
		  {
			 n++;
		     fr[n]=a;
		     r[n]=b;
		  }

		  if(a<0. || b<=0.)
		  {
			  printf("\n Input error. \n Each input point must be > 0 \n");
		  }
       }

//  Calculate the number of points

       
       fz=fr[n];
       nt=long((fz-fa)/bw)+1;

//	   printf("nt=%ld  \n",nt);

//  Calculate slopes between input points

       for(i=0; i<= (n-1); i++)
	   {
          a=log(r[i+1])-log(r[i]);
          b=log(fr[i+1])-log(fr[i]);
          s[i]=a/b;

//		  printf("%lf\n",s[i]);
       }

       psd[0]=r[0];
       f[0]=fr[0];

//       fprintf(pFile[1],"%15.6e %15.6e\n",f[0],psd[0] );


//  Interpolate

       if(nt > MAX2)nt=MAX2;
   
	   for(i=0; i<= nt; i++)
	   {
		  ff=fa+bw*i;

		  if( ff >= fr[0])
          {
			for( j=1; j<=n; j++)
			{
				if(ff < fr[j])
				{
					f[i]=ff;

					az=log10(r[j-1]);
					az+=s[j-1]*(log10(ff)-log10(fr[j-1]));
					psd[i]=pow(10.,az);

					fprintf(pFile[1],"%15.6e %15.6e\n",ff,psd[i] );
					break;
				}
				if(ff == fr[j])
				{
					f[i]=ff;
					psd[i]=r[j];

					fprintf(pFile[1],"%15.6e %15.6e\n",ff,psd[i] );
					break;
				}
			}
		  }
       }


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


}

void files(void)
{

	if(mflag==0)
	{
		printf("\n The input file must have two columns:  freq(Hz) & amplitude(units) ");

		printf("\n\n Enter input filename\n");
		scanf("%s",&filename[0]);
	
	}

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


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

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