/*  prog3e.c -- Burgers Equation using Lax Wendroff Method
    Written May 7, 2003 by Eric Olson for Mathematics 484
 
    For information about the Lax Wendroff scheme used here see 6.7 in 
    Hyperbolic Systems of Conservation Laws and the Mathematical Theory 
    of Shock Waves by Peter Lax and chapter 8 in Numerical Analysis for 
    Applied Mathematics Science and Engineering by Greenspan and Casulli.  
    Information about Burgers equations may be found in chapter 14 in 
    A First Course in the Numerical Analysis of Differential Equations 
    by Iserles. */

#define N 256
#include <math.h>
#include <stdio.h>

double u[N],f[N],dx,dt;

init(double u[N]){
    int i;
    for(i=0;i<N;i++){
        double x=dx*i;
        u[i]=2+3*sin(x);
    }
}

double flux(double x){
    return x*x;
}

force(double f[N],double u[N]){
    double uavg,fluxjm1,fluxj;
    int i;
    uavg=((u[N-1]+u[0])+dt*(flux(u[N-1])-flux(u[0]))/dx)/2;
    fluxjm1=flux(uavg);
    for(i=0;i<N;i++) {
        int ip1=(i+1)%N; 
        uavg=((u[i]+u[ip1])+dt*(flux(u[i])-flux(u[ip1]))/dx)/2;
        fluxj=flux(uavg);
        f[i]=(fluxjm1-fluxj)/dx;
        fluxjm1=fluxj;
    }
    return;
}

euler(double u[N],double tn,int n){
    int i,j;
    dt=tn/n;
    for(j=0;j<n;j++) {
        force(f,u);
        for(i=0;i<N;i++) {
            u[i]=u[i]+dt*f[i];
            if(fabs(u[i])>1e7){
                fprintf(stderr,"Blowup!\n");
                exit(1);
            }
        }
    }
    return;
}

vectwr(double u[N],double t){
    int i;
    printf("\n\n#At time %g...\n",t);
    for(i=0;i<N;i++){
        printf("%g %g\n",dx*i,u[i]);
    }
}

main(){
    double t;
    int i;
    dx=2*M_PI/N;
    printf("#prog3e.c -- Burgers Equation Using Lax Wendroff Method.\n");
    
    init(u);
    vectwr(u,0.0);

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