#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      100
#define a      1.0
#define FTCS    1  /* -- forward in time, centered in space -- */
#define UPWIND  2  /* -- choose depending on the sign of a  -- */



#define METHOD UPWIND  /* -- either UPWIND or FTCS -- */

/* *************************************************************** */
int main()
/* 
 *
 *  Solve the linear advection equation with a first-order
 *  method.
 *
 *  Last Modified  14 Nov 2011 by A. Mignone (mignone@ph.unito.it) 
 *
 ***************************************************************** */
{
  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;

/* -- default values -- */

  xbeg     = -.5;
  xend     =  .5;
  tstop    = 1.0;
  cfl      = 0.9;
  out_freq = 10;   /* -- output frequency -- */

/* -- 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]);
  }

  dt = cfl*dx/fabs(a);

/* -- start computation -- */

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

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

  /* -- set periodic boundary conditions -- */

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

    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 = sin(2.0*PI*x);
  u = exp(-x*x/0.02);
*/

  if (fabs(x) < 0.25){
    u = 1.0;
  }else{
    u = 0.0;
  }

  return(u);

}

/* *************************************************************** */
void Integrate (double *u0, double *u1, double dtdx, int ibeg, int iend)
/*
 *
 *
 ***************************************************************** */
{
  int i; 

#if METHOD == UPWIND
  if (a > 0.0){
    for (i = ibeg; i <= iend; i++){
      u1[i] = u0[i] - a*dtdx*(u0[i] - u0[i - 1]);
    }
  }else{
    for (i = ibeg; i <= iend; i++){
      u1[i] = u0[i] - a*dtdx*(u0[i+1] - u0[i]);
    }
  }
#elif METHOD == FTCS
    for (i = ibeg; i <= iend; i++){
      u1[i] = u0[i] - 0.5*a*dtdx*(u0[i+1] - u0[i - 1]);
    }
#endif
}

/* ******************************************************* */
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);
}