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

# include "flame_exact.h"

/******************************************************************************/

void flame_exact ( int n, double t[n], double *y )

/******************************************************************************/
/*
  Purpose:

    flame_exact() evaluates the exact solution for the flame ODE.

  Discussion:

    Evaluating the exact solution requires the Lambert W function.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 April 2021

  Author:

    John Burkardt

  Reference:

    Cleve Moler,
    Cleve's Corner: Stiff Differential Equations,
    MATLAB News and Notes,
    May 2003, pages 12-13.

  Input:

    int n: the number of times.

    double t[n]: the times.

  Output:

    double y[n]: the exact solution values.
*/
{
  double a;
  int i;
  int l;
  int nb;
  double t0;
  double tstop;
  double x;
  double y0;
/*
  Retrieve current parameter values.
*/
  flame_parameters ( NULL, NULL, NULL, &t0, &y0, &tstop );

  a = ( 1.0 - y0 ) / y0;

  nb = 0;
  l = 0;

  for ( i = 0; i < n; i++ )
  {
    x = a * exp ( a - ( t[i] - t0 ) );
    y[i] = 1.0 / ( lambert_w ( x, nb, l ) + 1.0 );
  }

  return;
}
/******************************************************************************/

void flame_parameters ( double *t0_in, double *y0_in, double *tstop_in, 
  double *t0_out, double *y0_out, double *tstop_out  )

/******************************************************************************/
/*
  Purpose:

    flame_parameters() returns parameters for flame_ode().

  Discussion:

    If input values are specified, this resets the default parameters.
    Otherwise, the output will be the current defaults.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 April 2024

  Author:

    John Burkardt

  Input:

    double *t0_in: the initial time.

    double *y0_in: the initial condition at time T0.

    double *tstop_in: the final time.

  Output:

    double *t0_out: the initial time.

    double *y0_out: the initial condition at time T0.

    double *tstop_out: the final time.
*/
{
  static double t0_default = 0.0;
  static double tstop_default = 200.0;
  static double y0_default = 0.01;
/*
  New values, if supplied on input, overwrite the current values.
*/
  if ( t0_in )
  {
    t0_default = *t0_in;
  }

  if ( y0_in )
  {
    y0_default = *y0_in;
  }

  if ( tstop_in )
  {
    tstop_default = *tstop_in;
  }
/*
  The current values are copied to the output.
*/
  if ( t0_out )
  {
    *t0_out = t0_default;
  }

  if ( y0_out )
  {
    *y0_out = y0_default;
  }

  if ( tstop_out )
  {
    *tstop_out = tstop_default;
  }

  return;
}
/******************************************************************************/

double lambert_w ( double x, int nb, int l )

/******************************************************************************/
/*
  Purpose:

    lambert_w() approximates the W function.

  Discussion:

    The call will fail if the input value X is out of range.
    The range requirement for the upper branch is:
      -exp(-1) <= X.
    The range requirement for the lower branch is:
      -exp(-1) < X < 0.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 June 2014

  Author:

    Original FORTRAN77 version by Andrew Barry, S. J. Barry, 
    Patricia Culligan-Hensley.
    This version by John Burkardt.

  Reference:

    Andrew Barry, S. J. Barry, Patricia Culligan-Hensley,
    Algorithm 743: WAPR - A Fortran routine for calculating real 
    values of the W-function,
    ACM Transactions on Mathematical Software,
    Volume 21, Number 2, June 1995, pages 172-181.

  Input:

    double X, the argument.

    int NB, indicates the desired branch.
    * 0, the upper branch;
    * nonzero, the lower branch.

    int L, indicates the interpretation of X.
    * 1, X is actually the offset from -(exp-1), so compute W(X-exp(-1)).
    * not 1, X is the argument; compute W(X);

  Output:

    double lambert_w, the approximate value of W(X).
*/
{
  double an2;
  static double an3;
  static double an4;
  static double an5;
  static double an6;
  static double c13;
  static double c23;
  static double d12;
  double delx;
  static double em;
  static double em2;
  static double em9;
  double eta;
  int i;
  static int init = 0;
  static int nbits;
  static int niter = 1;
  double reta;
  static double s2;
  static double s21;
  static double s22;
  static double s23;
  double t;
  static double tb;
  double temp;
  double temp2;
  double ts;
  double value;
  static double x0;
  static double x1;
  double xx;
  double zl;
  double zn;

  value = 0.0;

  if ( init == 0 )
  {
    init = 1;

    nbits = 52;

    if ( 56 <= nbits )
    {
      niter = 2;
    }
/*
  Various mathematical constants.
*/
    em = -exp ( -1.0 );
    em9 = -exp ( -9.0 );
    c13 = 1.0 / 3.0;
    c23 = 2.0 * c13;
    em2 = 2.0 / em;
    d12 = -em2;
    tb = pow ( 0.5, nbits );
    x0 = pow ( tb, 1.0 / 6.0 ) * 0.5;
    x1 = ( 1.0 - 17.0 * pow ( tb, 2.0 / 7.0 ) ) * em;
    an3 = 8.0 / 3.0;
    an4 = 135.0 / 83.0;
    an5 = 166.0 / 39.0;
    an6 = 3167.0 / 3549.0;
    s2 = sqrt ( 2.0 );
    s21 = 2.0 * s2 - 3.0;
    s22 = 4.0 - 3.0 * s2;
    s23 = s2 - 2.0;
  }

  if ( l == 1 )
  {
    delx = x;

    if ( delx < 0.0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "WAPR - Fatal error!\n" );
      fprintf ( stderr, "  The offset X is negative.\n" );
      fprintf ( stderr, "  It must be nonnegative.\n" );
      exit ( 1 );
    }

    xx = x + em;
  }
  else
  {
    if ( x < em )
    {
      return value;
    }
    else if ( x == em )
    {
      value = -1.0;
      return value;
    }
    xx = x;
    delx = xx - em;
  }
/*
  Calculations for Wp.
*/
  if ( nb == 0 )
  {
    if ( fabs ( xx ) <= x0 )
    {
      value = xx / ( 1.0 + xx / ( 1.0 + xx 
        / ( 2.0 + xx / ( 0.6 + 0.34 * xx ))));
      return value;
    }
    else if ( xx <= x1 )
    {
      reta = sqrt ( d12 * delx );
      value = reta / ( 1.0 + reta / ( 3.0 + reta / ( reta 
        / ( an4 + reta / ( reta * an6 + an5 ) ) + an3 ) ) ) 
        - 1.0;
      return value;
    }
    else if ( xx <= 20.0 )
    {
      reta = s2 * sqrt ( 1.0 - xx / em );
      an2 = 4.612634277343749 * sqrt ( sqrt ( reta + 
        1.09556884765625 ));
      value = reta / ( 1.0 + reta / ( 3.0 + ( s21 * an2 
        + s22 ) * reta / ( s23 * ( an2 + reta )))) - 1.0;
    }
    else
    {
      zl = log ( xx );
      value = log ( xx / log ( xx 
        / pow ( zl, exp ( -1.124491989777808 / 
        ( 0.4225028202459761 + zl )) ) ));
    }
  }
/*
  Calculations for Wm.
*/
  else
  {
    if ( 0.0 <= xx )
    {
      return value;
    }
    else if ( xx <= x1 )
    {
      reta = sqrt ( d12 * delx );
      value = reta / ( reta / ( 3.0 + reta / ( reta / ( an4 
        + reta / ( reta * an6 - an5 ) ) - an3 ) ) - 1.0 ) - 1.0;
      return value;
    }
    else if ( xx <= em9 )
    {
      zl = log ( -xx );
      t = -1.0 - zl;
      ts = sqrt ( t );
      value = zl - ( 2.0 * ts ) / ( s2 + ( c13 - t 
        / ( 270.0 + ts * 127.0471381349219 )) * ts );
    }
    else
    {
      zl = log ( -xx );
      eta = 2.0 - em2 * xx;
      value = log ( xx / log ( -xx / ( ( 1.0 
        - 0.5043921323068457 * ( zl + 1.0 ) ) 
        * ( sqrt ( eta ) + eta / 3.0 ) + 1.0 )));
    }

  }

  for ( i = 1; i <= niter; i++ )
  {
    zn = log ( xx / value ) - value;
    temp = 1.0 + value;
    temp2 = temp + c23 * zn;
    temp2 = 2.0 * temp * temp2;
    value = value * ( 1.0 + ( zn / temp ) * ( temp2 - zn ) 
      / ( temp2 - 2.0 * zn ) );
  }

  return value;
}