/*        ATK IV - Numerical Programming (2006)       */
/*     Excercise 12.2 -  FFT of Gaussian Function     */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>
#define N 512

float g(float t, float t0, float s){return exp(-pow((t-t0)/s,2)/2);}

int main(void){
  float t0 = 0.0, t, t1, t2, dt, df, f;
  float s = 1.0;
  float data[2*N+1];
  int i;
  float sum=0.0;
  float powert,powerf;
  FILE *file;

  file=fopen("print.txt","w");

  t1=t0-s*5;   /* I just tested that this kind of offset from the mean is sufficient*/
  t2=t0+s*5;
  dt = (t2-t1)/N;
  df = 1.0/(N*dt);

   for(i=1;i<=N;i++){
    t = t1 + (i-1)*dt;
    data[2*i-1]=g(t,t0,s);
    data[2*i]=0.0;
    sum += SQR(g(t,t0,s));
   }

   powert=dt*sum;
   four1(data,N,1);

   sum=0.0;
   for(i=N/2+1;i<=N-1;i++){
     f = (i-N)*df;
     fprintf(file,"%f %f %f\n",f,fabs(dt*data[2*i+1]),dt*data[2*i+2]);
     sum+=SQR(data[2*i+1])+SQR(data[2*i+2]);
   }

   for(i=0;i<=N/2;i++){
     f = i*df;
     fprintf(file,"%f %f %f\n",f,fabs(dt*data[2*i+1]),dt*data[2*i+2]);
     sum+=SQR(data[2*i+1])+SQR(data[2*i+2]);
   }

   powerf = sum*dt/N;

   printf("\n%f %f\n",powert,powerf);
   fclose(file);

   return 0;
}