#include <fstream.h>
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include <stdlib.h>

// prototype of the exact solution
float burgex(float x, float t, float u0);

int main(void)

 /**********************************************************************\
 *                                                                      *
 * Solve the Burger's equation problem:                                 *
 *                                                                      *
 *   u_t + u u_x = 0    u(x,0) = u0 + sin(x)   u(x,t) = u(x + 2 pi,t)   *
 *                                                                      *
 * This program uses a nonconservative upwind scheme                    *
 *                                                                      *
 * Compile using: cpp -o burgnc burgnc.cpp burgex.cpp                   *
 *                                                                      *
 * Run using:          burgnc                                           *
 *                                                                      *
 * Plot results using:                                                  *
 * ~dommelen/bin/xyplot burgnc.plt 2 ' ' ' ' ' ' 2 '0,100,0,100,.4'     *
 * which creates a printable postscript file named dump.ps              *
 *                                                                      *
 * To use your own plotting software, change the format of file         *
 * burgnc.plt below as needed.                                          *
 *                                                                      *
 \**********************************************************************/

{

  // Declarations

  // declared storage size in the j-direction
  const int J_DIM=1000;

  // average initial velocity
  float u0;

  // time step index n
  int n;
  // index at which the velocities u^n are actually stored
  int n_s;
  // index at which the velocities u^{n+1} are actually stored
  int np1_s;
  // next value of n at which to print out diagnostics
  int nprint;

  // time
  float t;

  // actual number of mesh cells
  int j_max;
  // mesh point index;
  int j;

  // x-positions of the mesh points
  float x[J_DIM];

  // velocity values.  only two time levels are stored at any point
  float u[J_DIM][2];

  // mesh cell size Delta x and time step Delta t
  float dx,dt;

  // maximum Courant number C = u Delta t / Delta x
  float C;

  // output time interval;
  float t_out;
  // output time index and number of output times
  int out,outmx;

  // pointwise and rms errors in the solution
  float err,errrms;

  // minimum and maximum velocity
  float u_min,u_max;

  // input file
  ifstream infile;
  // output and plot files
  ofstream outfile, plotfile;

  // character temporary to flush comments from the input file
  char tmp;

  // constants
  float pi;

  pi=4.*atan(1.);

  // executable statements

  // open files

  // input file
  infile.open("burg.in");

  // output file
  outfile.open("burgnc.out");

  // plot file
  plotfile.open("burgnc.plt");

  // read input file and write to output

  infile >> u0;
  outfile << "average initial velocity: " << u0 << endl;
  while(infile.get(tmp) && tmp!='\n');

  infile >> j_max;
  outfile << "number of mesh cells j_max: " << j_max << endl;
  if(j_max>J_DIM-2)
    {
      cout << "*** declared dimension exceeded" << endl;
      exit(1);
    }
  if(j_max<2)
    {
      cout << "*** number of mesh cells should be at least 2" << endl;
      exit(1);
    }
  while(infile.get(tmp) && tmp!='\n');

  infile >> C;
  outfile << "maximum Courant number u Delta t/Delta x: " << C << endl;
  if(C<=0.)
    {
      cout << "*** Courant number must be positive" << endl;
      exit(1);
    }
  while(infile.get(tmp) && tmp!='\n');

  infile >> t_out;
  outfile << "output time interval: " << t_out << endl;
  if(t_out<=0.)
    {
      cout << "*** output time interval must be positive" << endl;
      exit(1);
    }
  while(infile.get(tmp) && tmp!='\n');

  infile >> outmx;
  outfile << "number of output time intervals to compute: " << outmx << endl;
  if(outmx<=0)
    {
      cout << "*** number of output times must be positive" << endl;
      exit(1);
    }
  while(infile.get(tmp) && tmp!='\n');

  infile.close();

  // initialize the parameters

  // mesh size
  dx=2.*pi/j_max;

  // output count
  out=1;

  // print value of n
  nprint=1;

  // initial velocities

  // set all initial velocities and x-positions
  for (j=0; j<=j_max; j++)
    {
      x[j]=j*dx;
      u[j][0]=u0+sin(x[j]);
    }

  // initialize time
  t=0.;

  // main loop over all time steps

  for (n=0; n<=123456789; n++)
    {

      // where to store the velocities at n and n+1
      n_s=n%2;
      np1_s=(n+1)%2;

      // perform output;

      while (t>=out*t_out)
	{

	  // write computed velocities to the plot file;
	  plotfile << -j_max-1 << endl;
	  for (j=0; j<=j_max; j++)
	    plotfile << x[j] << ',' << u[j][n_s] << endl;

	  // write the exact solution to the plot file as a 201 point curve;
	  plotfile << "201    " << j_max << t << endl;
	  for (j=0; j<=200; j++) plotfile
	    << 0.005*2.*pi*j << ',' << burgex(0.005*2.*pi*j,t,u0) << endl;

	  out=out+1;
	}

      if (out>outmx) break;

      // time step
      dt=1.e30;
      for (j=0; j<j_max; j++)
	{
	  if (u[j][n_s]!=0.)
	    {
	      if (C*dx/fabs(u[j][n_s])<dt) dt=C*dx/fabs(u[j][n_s]);
	    }
	}

      // increment time
      t=t+dt;

      // set the j=jmax+1 value of u using periodicity
      u[j_max+1][n_s]=u[1][n_s];

      // compute the velocities
      for (j=1; j<=j_max; j++)
	{
	  if(u[j][n_s]>=0.)
	    u[j][np1_s]=u[j][n_s]-dt/dx*u[j][n_s]*(u[j][n_s]-u[j-1][n_s]);
	  else
	    u[j][np1_s]=u[j][n_s]-dt/dx*u[j][n_s]*(u[j+1][n_s]-u[j][n_s]);
	}

      // set the boundary velocity at 0 using periodicity
      u[0][np1_s]=u[j_max][np1_s];

      // print diagnostics only occasionally
      if (n+1!=nprint) continue;
      nprint=nprint*2;

      // compute the rms error and minimum and maximum velocity
      errrms=0.;
      u_min=u[0][np1_s];
      u_max=u[0][np1_s];
      for (j=0; j<=j_max-1; j++)
	{
	  err=fabs(u[j][np1_s]-burgex(x[j],t,u0));
	  errrms=errrms+err*err;
	  if (u[j][np1_s]<u_min) u_min=u[j][np1_s];
	  if (u[j][np1_s]>u_max) u_max=u[j][np1_s];
	}
      errrms=sqrt(errrms/j_max);
      cout
	<< "n =" << setw(3) << n+1
	<< setiosflags(ios::fixed) << setprecision(3)
	<< "   t =" << setw(6) << t
	<< setprecision(5)
	<< "   u-range: " << setw(8) << u_min << " to" << setw(8) << u_max
	<< "   error: " << setw(7) << errrms << endl
	<< resetiosflags(ios::fixed) << setprecision(0);
      outfile
	<< "n = " << n+1
	<< "   t = " << t
	<< "   u-range: " << u_min << '-' << u_max
	<< "   error: " << errrms << endl;

    }

  return 0;
}