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

float exact

 /********************************************************************\
 * The exact solution to the unidirectional flow Navier-Stokes        *
 * equation:                                                          *
 *                                                                    *
 *                          v_t = nu v_xx                             *
 *                                                                    *
 * with boundary conditions:                                          *
 *                                                                    *
 *              v(0,t) = V_1            v(l,t) = V_2                  *
 *                                                                    *
 * and initial condition:                                             *
 *                                                                    *
 *             v(x,0) = V_0  4 (x/l) (1 - (x/l)) + V_2                *
 *                                                                    *
 * the exact solution is normally found by summing the Separation     *
 * of Variables series to the desired accuracy tol.                   *
 *                                                                    *
 \********************************************************************/

(
 float x,     // value of x at which the velocity is desired
 float t,     // value of t at which the velocity is desired
 float nu,    // coefficient of kinematic viscosity
 float l,     // length of the flow region
 float V_0,   // constant in the given initial condition
 float V_1,   // left boundary velocity
 float V_2,   // right boundary velocity
 float tol    // desired accuracy (specify 0. for maximum precision)
)

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

  // Fourier mode number
  int k;
  // the value of k as a float
  float rk;

  // the amplitude of the k-th Fourier mode
  float amplit;

  // the total k-th term of the Fourier series
  float term;

  // the sum of all the terms of the Fourier series so far
  float sum;

  // whether we are done summing
  int done;

  // pi  (the correct value 3.14... of pi will be set later)
  static float pi=0.;

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

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

  // on the boundaries, the solution is given
  if(x<=0.)return V_1;
  if(x>=l)return V_2;

  // else, at t=0, the solution is also given
  if(t==0.)return 4.*V_0*(x/l)*(1.-(x/l))+V_2;

  // initialize the exact solution to zero
  sum=0.;

  // we will be computing the exact solution in 3 parts; the first
  // part is the SOV series due to the V_0 part of the initial
  // condition alone:

  // initial value of k as a float number (why do we need it?)
  rk=1.;

  // not done yet
  done=0;

  // sum up to 1,000,000 Fourier modes together
  for(k=1; k<2000000 && !done; k=k+2)
    {

      // the k-th term of the sum, except for the sin(k pi x/l) factor
      amplit=(32.*V_0/(pi*pi*pi))/(rk*rk*rk)*exp(-t*nu*(pi*pi/(l*l))*rk*rk);

      // the total k-th term to add (why is pi*x/l between parentheses??)
      term=amplit*sin(rk*(pi*x/l));

      // done if we have have the requested accuracy
      if(fabs(amplit*rk)<=tol)done=1;
      // Why look at the amplitude alone, not the complete 'term' below??
      // Why is the *rk there???

      // done in any case if maximum possible accuracy has been achieved
      if(amplit+sum==sum)done=1;
      // Does this make sense?????????????????????????  Explain!

      // add the k-th term to the Fourier sum
      sum=sum+term;

      // increment the float value of k (is this correct?)
      rk=rk+2.;

    }

  if(!done)
    {
      cout<<"*** exact: First sum not converged after 1,000,000 terms";
      exit(1);
    }

  // the second part of the solution is that due to the V_2 initial
  // condition, V_2 boundary condition at x=l, and also assuming a
  // V_2 boundary condition at x=0.  The solution is simply V_2.
  
  sum=sum+V_2;

  // the last part is the solution due to the missing V_1-V_2 boundary
  // condition at x=0, with zero initial conditions and zero boundary
  // condition at x=l.  This solution itself consists of two parts;
  // the first is the large time solution satisfying the inhomogeneous
  // boundary condition:

  sum=sum+(V_1-V_2)*(1.-(x/l));

  // the other, final part is the SOV solution to the problem with
  // homogeneous boundary conditions and a -(V_1-V_2)*(1.-(x/l))
  // initial condition

  // initial value of k as a float number (why do we need it?)
  rk=1.;

  // not done yet
  done=0;

  // sum up to 1,000,000 Fourier modes together
  for(k=1; k<1000000 && !done; k=k+1)
    {

      // the k-th term of the sum, except for the sin(k pi x/l) factor
      amplit=-(2.*(V_1-V_2)/pi)/rk*exp(-t*nu*(pi*pi/(l*l))*rk*rk);

      // the total k-th term to add (why is pi*x/l between parentheses??)
      term=amplit*sin(rk*(pi*x/l));

      // done if we have have the requested accuracy
      if(fabs(amplit*rk)<=tol)done=1;

      // done in any case if maximum possible accuracy has been achieved
      if(amplit+sum==sum)done=1;

      // add the k-th term to the Fourier sum
      sum=sum+term;

      // increment the float value of k (is this correct?)
      rk=rk+1.;

    }

  if(!done)
    {
      cout<<"*** exact: 2nd sum not converged after 1,000,000 terms";
      exit(1);
    }

  return sum;
}