/*  prog2a.c -- Solve viscous Burgers Equation
    Written Apr 6, 2003 by Eric Olson for Mathematics 484
 
    For information about the Discrete Fourier Transform please see
    Chapter 12 in Numerical Recipies by Press, Flannery, Teukolsky and
    Vetterling.  Information about the heat and Burgers equations may
    be found in chapters 13 and 14 in A First Course in the Numerical
    Analysis of Differential Equations by Iserles. */

/* Note N must be a power of 2 */
#define N 8
#include <math.h>
#include <stdio.h>

typedef struct { double a, b; } complex;
complex u[N],uu[N],f[N];

extern void fft_(complex *,int *,int *);
void fft(complex u[N],int n,int k){
    fft_(u,&n,&k);
}

init(complex u[N]){
    int i;
    for(i=0;i<N;i++){
        double x=2*M_PI*i/N;
        double xx=M_PI-x; xx*=xx;
        u[i].a=8*exp(-xx);
        u[i].b=0.0;
    }
    fft(u,N,1);
}

force(complex f[N],complex u[N]){
    int i;
    for(i=0;i<N;i++) uu[i]=u[i];
    fft(uu,N,-1);
    for(i=0;i<N;i++) {
        uu[i].a*=uu[i].a;
        uu[i].b=0.0;
    } 
    fft(uu,N,1);
    for(i=0;i<N;i++){
        int k;
        if(2*i>N) k=i-N; else k=i;
        f[i].a=-k*(k*u[i].a-uu[i].b);
        if(2*i==N) f[i].b=0.0;
        else f[i].b=-k*(k*u[i].b+uu[i].a);
    }
    return;
}

euler(complex u[N],double tn,int n){
    double dt=tn/n;
    int i,j;
    for(j=0;j<n;j++) {
        force((complex *)f,u);
        for(i=0;i<N;i++) {
            u[i].a=u[i].a+dt*f[i].a;
            u[i].b=u[i].b+dt*f[i].b;
        }
    }
    return;
}

vectwr(complex u[N],double t){
    int i;
    for(i=0;i<N;i++) uu[i]=u[i];
    fft(uu,N,-1);
    printf("\n\n#At time %g...\n",t);
    for(i=0;i<N;i++){
        int k;
        if(2*i>N) k=i-N; else k=i;
        printf("%g %g %g %d %g %g\n",
            2*M_PI*i/N,uu[i].a,uu[i].b,k,u[i].a,u[i].b);
    }
}

main(){
    double t;
    complex f[N];
    int i;
    printf("#prog2a.c -- Compute viscous Burger's Equation\n");
    
    init(u);
    vectwr(u,0.0);

    t=0.0;
    for(i=1;i<=8;i++){
        euler(u,0.125,1024);
        vectwr(u,i*0.125);
    }
    return 0;
}