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