//**************************************************************************
//****  Following is a sample program  named 'sample.C'.
//****  Compile the program  (Command is : g++ sample.C -lm )
//****  and run it  (Command is : a.out  )
//****  Then the result will be stored in a file named 'outputfilename'.
//**************************************************************************

#include <iostream.h>
#include <fstream.h>
#include <stream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

//**************************************************************************
//**** Constants to be used in FDM
//**************************************************************************
const double L1=-2;
const double L2=2;
const double T=0.33;                              // final time
const int MAX=641;                            // MAX=(L2-L1)/dx+1
const double v0=2.0;
const double nv=4.0;

const double dx=(L2-L1)/(MAX-1);                // dx
const double dt=dx/nv;                          // dt

const int OSCILATION=1;
//**************************************************************************
//**** Initial data setup
//**************************************************************************
double initial_function(double x){
  if (x<-1) return(0);
  if (x<1) return(1+sin(PI*x)*0.5);
  return(0);
}
void initialize(double U[][2]){
  int i;
  double x;

  for(i=0;i<MAX;++i){
    x=L1+i*dx;
    U[i][0]=initial_function(x);
    U[i][1]=initial_function(x);
  }
}
//********************************************************
//******* 
double w(double a,double b){
  return(b-a);
}
double w(double a,double b,double c){
  return(w(b,c)-w(a,b));
}
double w(double a,double b,double c,double d){
  return(w(b,c,d)-w(a,b,c));
}
double w(double a,double b,double c,double d,double e){
  return(w(b,c,d,e)-w(a,b,c,d));
}

//**************************************************************************
//**** Average over a mesh [x_j,x_j+1]
//**************************************************************************
double AVR(double U[][2],int j,int n){
}

//**************************************************************************
//**** numerical flux function
//**************************************************************************
double nuflux(double U[][2],int j,int n){
}
//************************************************************
//**** In the main function we have I/O's and recursive  *****
//**** calls of FDM 
//************************************************************
void main(){
  filebuf fi;
  int j,m,l,next,now;
  double U[MAX][2];

  initialize(U);

  next=0;
  for(m=0;m*dt<T;m=m+2){

    now=next;
    next=(now+1)%2;
    for(j=1;j<MAX-1;++j){
      U[j][next]=AVR(U,j,now)-nv*(nuflux(U,j+1,now)-nuflux(U,j,now));
    }

    now=next;
    next=(now+1)%2;
    for(j=1;j<MAX-1;++j){
      U[j+1][next]=AVR(U,j,now)-nv*(nuflux(U,j+1,now)-nuflux(U,j,now));
    }
  }

  ofstream output("output", ios::out);
  for(j=0;j<MAX;++j){
     output<<form("%f",U[j][m%2])<<endl;
  }

return(1);
}