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

#define MAX   262144
#define MAXP1 131072



void files(void);
void read_data(void);
void input(void);


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

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


floatcomplex c;

floatcomplex y[MAX],H[MAX];

float freq[MAX];

doublecomplex num, den;

double f, ff, fn, damp, rho, z;

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

long i, idamp,last;


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


void main()
{

	 printf("\n transfer_fft.cpp   version 1.0,  July 20, 2003 \n");
	 printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

	 printf("\n This program multiplies the Fourier transform of a ");
	 printf("\n base input function by the transfer function of a ");
	 printf("\n single-degree-of-freedom system. \n\n");

	 printf("\n The Fourier transform must be complex and full-sided. \n");

	 files();

	 input();

     read_data();


	 long half = long(double(last)/double(2) );


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

        ff = freq[i];

	    rho=ff/fn;

        num.r=1.;
	    num.i=2.*damp*rho;

        den.r=1-pow(rho,2.);
	    den.i=num.i;


	    COMPLEXDIV(num.r, num.i, den.r, den.i, H[i].r, H[i].i) 

	 }

	 for(i=0; i<half; i++)
	 {
		 H[last-i-1].r =  H[i].r; 
		 H[last-i-1].i = -H[i].i;
	 }


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

	 	COMPLEXMUL(y[i].r, y[i].i, H[i].r, H[i].i, c.r, c.i);
		
        fprintf(pFile[1]," %12.5g  %12.4g  %12.4g\n",freq[i],c.r,c.i);

     }
	 fclose( pFile[1] );


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


}
void files(void)
{

    printf("\n The Fourier transform must have 2^N points, where N is an integer. \n");

    printf( "\n The transform must contain three columns: \n");
    printf( " frequency(Hz)  real accel(G)  imaginary accel(G)    \n");


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


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

 
    i=0;


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


    printf("\n Enter the output filename: \n");

	scanf("%s",filename[1]);
	pFile[1]=fopen(filename[1], "w");

}
void read_data(void)
{
	printf("\n\n reading Fourier transform ");
     
	i=0;
    while( ( fscanf(pFile[0]," %f %f %f",&freq[i], &y[i].r, &y[i].i) ) > 0 )
    {
		i++;
    }
	fclose( pFile[0] );

	last=i;

	printf("\n\n points read = %ld \n",last);

	long iflag = 0;

    long number;

	number=2;

	while(1)
	{
    
		if(number >= MAX ){break;}


		if(last==number){iflag=1;}

	    number*=2;

	}

	if( iflag == 0 )
	{

        printf("\n Error:  Number of points is not equal to 2^N \n");

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

		exit(1);
	}

}
void input(void)
{

	printf("\n\n Enter natural frequency (Hz) \n");
	scanf("%lf",&fn);

   	printf("\n Enter damping format. \n 1= damping ratio   2= Q \n");
	scanf("%ld",&idamp);

	if(idamp == 1 )
	{
		printf("\n Enter damping ratio (typically 0.05) ");
		scanf("%lf",&damp);
	}
	else
    {
		printf("\n Enter amplification factor Q (typically 10) ");
		scanf("%lf",&damp);
		damp = 1./(2.*damp);
	}
}