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

float burgex

 /**********************************************************************\
 *                                                                      *
 * 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)   *
 *                                                                      *
 \**********************************************************************/

(
 float x,   // position at which the exact solution is required
 float t,   // time at which the exact solution is required
 float u0   // average initial velocity
)

{
  /************************ local variables *************************/

  // value of u
  float u;

  // value of xbar = x - u0 t and a reduced value
  float xbar,xbar2;

  // xi-value, the change in xi_value, and the previous change 
  float xi,dxi,dxiprv;

  // iteration counter
  int iter;

  // pi
  static float pi=0.;

  /********************* executable statements **********************/

  // compute pi if it has not been done yet
  if (pi==0.) pi=4.*atan(1.);

  // if t=0, trivial
  if (t<=0.) return u0+sin(x);

  // find xbar, and reduce to the interval from 0 to 2 pi
  xbar=x-u0*t;
  while (xbar<0.) xbar=xbar+2.*pi;
  while (xbar>=2.*pi) xbar=xbar-2.*pi;

  // restrict xbar to from 0 to pi
  xbar2=xbar;
  if(xbar>pi)xbar2=2.*pi-xbar;

  // use Newton iteration to solve (x-xi)/t - sin(xi) = 0 for xi

  // initialize xi
  xi=0.;
  dxi=1.e30;

  for(iter=1; iter<=12345; iter++)
    {

      // compute change in xi
      dxiprv=dxi;
      dxi=-((xbar2-xi)/t-sin(xi))/(-1./t-cos(xi));

      // check for convergence
      if (fabs(dxi)>=fabs(dxiprv) || dxi==0.)
	{
	  u=(xbar2-xi)/t;
	  if (xbar!=xbar2) u=-u;
	  u=u+u0;
	  return u;
	}

      // increment xi
      xi=xi+dxi;

    }

  // no convergence.  print an error message and stop.
  cout << "*** No convergence in burgex\n" << endl;
  exit(1);

}