#include <stdio.h>
#include <math.h>
#define pi 3.141592653589793

/* User-defined complex structure */
typedef struct FCOMPLEX {double r,i;} fcomplex;

/* Complex functions */
fcomplex csum(fcomplex a, fcomplex b);
fcomplex cmul(fcomplex a, fcomplex b);
fcomplex complexx(double re, double im);
fcomplex csimpson(double a, double b, fcomplex (*func)(double, double), int N, double w);
fcomplex cftintegrand(double w, double t);
fcomplex c2,c3,c4;

/* Additional math routines NOT in math.h */
int iseven(int i);	/* Checks whether integer argument is even(=1) or odd(=0)  */ 
double dabs(double x);   /* Returns the absolute value of a double argument */
int main()
{
   FILE *out;
   int samples_per_second=2; /* defines T as 1/(samples_per_sec) */
   int n=1,N;                /* n=time index, N=number of samples total */
   double dt,w;              
   int i;
   
// Initial and final time values could be input from screen 
   double t_ll=-100.0, t_ul=100.0;  
// Set values used in integration method and elsewhere - these are unchanged 
   c2=complexx(2.0,0.0);c3=complexx(3.0,0.0);c4=complexx(4.0,0.0);
// H is Fourier transform of given function h(cintegrand.c)
   fcomplex H;
   H=complexx(0.0,0.0);

// Prepare to write for plotting
   printf("\nBrute force complex Fourier transform algorithm");
   if ((out=fopen("output.dat","w"))==NULL) {
      printf("\nCannot open file for output\n");
      getchar();
      return(0);  		
   }

// Establish time and frequency controls
   N=samples_per_second*(int)(t_ul-t_ll)+1; /* # of data points */
   dt=(t_ul-t_ll)/(double)N;  // Time sampling period, T

// Loop through the frequencies and evaluate time integrals   
   n=-1.0*(N/2.0);
   for (i=1;i<=N;i++)
   {
      w=2.0*pi*n/(N*dt);
      H=csimpson(t_ll,t_ul,cftintegrand,N,w);
      fprintf(out,"%5.3f %5.3f\n",w,H.r);
      if ((i%(N/10))==0) printf("\nThrough frequency %6.3f",w);
      n++;
   }
   fclose(out);
   printf("\nProgram complete without known error.\n");
   return(1);
}

fcomplex csimpson(double a, double b, fcomplex    
   (*func)(double, double), int N, double w)
/* Computes simpsons rule integration for a complex integrand */
{
  double dt,t;
  int i;
  fcomplex integral, cdt3;
  integral=complexx(0.0,0.0);

  dt = (b-a)/(double)(N);
  cdt3=complexx(dt/3.0,0.0);

  t = a;
  
  for(i=1;i<=N;i++) {
    if(i == 1 || i == N) 
  integral = csum(integral,(*func)(w,t));
    else if(iseven(i)) 
  integral = csum(integral,cmul(c4,(*func)(w,t)));
    else integral = csum(integral,cmul(c2,(*func)(w,t)));
    t += dt;
  }
  integral = cmul(integral,cdt3);
  return(integral);
} 

fcomplex cftintegrand(double w, double t)
{
   fcomplex fterm, h;
   
/* Unit step function */
   if ((t<=5.0)&&(t>=-5.0)) h=complexx(1.0,0.0);
   else h=complexx(0.0,0.0);
   
   //printf("%f   %f   %f\n",t,h.r,h.i);
   fterm=complexx(cos(w*t),sin(w*t));
   return(cmul(h,fterm));
}

fcomplex complexx(double re, double im)
{
    fcomplex c;
    c.r = re;
    c.i = im;
    return c;
}

fcomplex cmul(fcomplex a, fcomplex b)
{
    fcomplex c;
    c.r=a.r*b.r-a.i*b.i;
    c.i=a.i*b.r+a.r*b.i;
    return c;
}

fcomplex csum(fcomplex a, fcomplex b)
{
    fcomplex c;
    c.r=a.r+b.r;
    c.i=a.i+b.i;
    return c;
}

int iseven(int i)
{
   if (i%2==0) return(1);
   else return(0);
}

double dabs(double x)
{
    if (x>=0.0) return(x);
    else return(-1.0*x);
}