#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

double Initial_Condition (double x);
void   Integrate (double *u0, double *u1, double dtdx, int ibeg, int iend);

#define PI      3.14159265358979
#define NGHOST   2
#define NX      4000
/* *************************************************************** */
int main()
/* 
 *
 * 
 *
 ***************************************************************** */
{
  int    i, nstep, out_freq;
  int    ibeg, iend;
  double xbeg, xend;
  double x [NX + 2*NGHOST], dx;
  double u0[NX + 2*NGHOST], u1[NX + 2*NGHOST];
  double t, tstop, dt, cfl, dtdx;
  double umax;

/* -- default values -- */

  xbeg     = -5.0;
  xend     = 5.0;
  tstop    = 2.0;
  cfl      = 0.9;
  out_freq = 400;

/* -- make grid & initial conditions -- */

  ibeg = NGHOST;
  iend = ibeg + NX - 1;

  dx = (xend - xbeg)/(double)NX;
  for (i = 0; i <= iend + NGHOST; i++){
    x[i]  = xbeg + (0.5 + i - ibeg)*dx;
    u0[i] = Initial_Condition (x[i]);
  }

/* -- start computation -- */

  t     = 0.0; nstep = 0;
  while (t < tstop){

  /* -- set time step -- */

    umax = 0.0;
    for (i = ibeg; i <= iend; i++){
      if (fabs(u0[i]) > umax) umax = fabs(u0[i]);
    }

    dt = dx/umax;

    printf ("step %d, dt = %f, t = %f\n",nstep, dt, t);
    if (nstep%out_freq == 0) Output (x, u0, ibeg, iend);

  /* -- set outflow boundary conditions -- */

    for (i = 1; i <= NGHOST; i++){
      u0[ibeg - i] = u0[ibeg];
      u0[iend + i] = u0[iend];
    } 

    Integrate (u0, u1, dt/dx, ibeg, iend);

    t += dt; nstep++;
    for (i = ibeg; i <= iend; i++) u0[i] = u1[i];
  }
  Output (x, u0, ibeg, iend);

  return (0);
}

/* *************************************************************** */
double Initial_Condition (double x)
/* 
 *
 *
 ***************************************************************** */
{
  int i;
  double u;

  u = exp(-x*x);
/*
  if (x < 0.) u = 0.0;
  else u = 1.0;
*/
  return(u);
}

/* *************************************************************** */
void Integrate (double *u0, double *u1, double dtdx, int ibeg, int iend)
/*
 *
 *
 ***************************************************************** */
{
  int i; 
  double flux[NX + 2*NGHOST]; /* -- conservative flux -- */
  double S;  /* -- shock speed -- */
  double us; /* -- solution to the Riemann problem -- */
  
/* -- compute fluxes by solving Riemann problem -- */

  for (i = ibeg - 1; i <= iend; i++){
    if (u0[i] > u0[i+1]) {  /* -- this is a shock wave -- */

      S = 0.5*(u0[i] + u0[i+1]);
      if (S >= 0.0) {
        us = u0[i];
      }else{
        us = u0[i+1];
      } 

    }else{  /* -- this is a rarefaction wave -- */
 
      if (u0[i] >= 0.0)  
        us = u0[i];
      else if (u0[i] <= 0.0)        
        us = u0[i+1];
      else 
        us = 0.0;
 
    }

    flux[i] = 0.5*us*us;
  }

/* -- update solution -- */

  for (i = ibeg; i <= iend; i++){
    u1[i] = u0[i] - dtdx*(flux[i] - flux[i-1]);
  }
}

/* ******************************************************* */
int Output (double *x, double *u, int ibeg, int iend)
/* 
 *
 *
 *
 ********************************************************* */
{
  int i;
  char name[128];
  FILE *fout;
  static int nfile = 0;

  sprintf (name,"data.%04d.out",nfile);
  printf ("Opening %s for writing...\n",name);

  fout = fopen(name,"w");
  for (i = ibeg; i <= iend; i++){
    fprintf (fout, "%12.6e  %12.6e\n",x[i], u[i]);
  }
  fclose(fout);
 
  nfile++;
  return (nfile);
}