/*      ATK IV - Numerical Programming (2006)      */
/*     Excercise 12.3 -  Correlation Using FFT     */

#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);}
float h(float t, float t0, float s,  float alfa){
  long idum;
  return g(t,t0,s)+alfa*ran1(&idum);
}

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

  file=fopen("print2.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);
  alfa = 1.0;
  tlag = 0.0;

   for(i=1;i<=N;i++){
    t = t1 + (i-1)*dt;
    datag[2*i-1]=g(t,t0,s);
    datag[2*i]=0.0;
    datah[2*i-1]=h(t+tlag,t0,s,alfa);
    datah[2*i]=0.0;
   }

   four1(datag,N,1);
   four1(datah,N,1);

   for(i=1;i<=N;i++){
     datah[2*i]=-datah[2*i];
   }

   for(i=1;i<=N;i++){
     data3[2*i-1]=datag[2*i-1]*datah[2*i-1]-datag[2*i]*datah[2*i];
     data3[2*i]=datag[2*i-1]*datah[2*i]+datag[2*i]*datah[2*i-1];
   }

   four1(data3,N,-1);

   for(i=1;i<=N;i++){
     t = t1 + (i-1)*dt;
     fprintf(file,"%f %f %f\n",t,data3[2*i-1]*dt*2.0/N,data3[2*i]*dt*2.0/N);
   }

   fclose(file);

   return 0;
}