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

# include "lambert_w.h"

int main ( );
void lambert_w_test01 ( );
void lambert_w_values ( int *n_data, double *x, double *fx, int *b );
void timestamp ( );

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

int main ( )

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

    lambert_w_test() tests lambert_w().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    21 June 2023

  Author:

    John Burkardt.
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "lambert_w_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test lambert_w().\n" );

  lambert_w_test01 ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "lambert_w_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void lambert_w_test01 ( )

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

    lambert_w_test01() compares lambert_w() to stored values, and a built in function

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    21 June 2023

  Author:

    John Burkardt.
*/
{
  int b;
  double fx1;
  double fx2;
  int n_data;
  double x;

  printf ( "\n" );
  printf ( "lambert_w_values_test():\n" );
  printf ( "  Compare a stored value to the lambert_w(X) function.\n" );
  printf ( "\n" );
  printf ( "                X            Branch           W(X)    lambert_w(X)\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    lambert_w_values ( &n_data, &x, &fx1, &b );

    if ( n_data == 0 )
    {
      break;
    }

    fx2 = lambert_w ( x, b, 0 );

    printf ( "  %24.16g  %2d  %24.16g  %24.16g\n", x, b, fx1, fx2 );
  }
  return;
}
/******************************************************************************/

void lambert_w_values ( int *n_data, double *x, double *fx, int *b )

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

    lambert_w_values() returns some values of the Lambert W function.

  Discussion:

    The function W(X) is defined implicitly by:

      W(X) * e^W(X) = X

    The function is defined for -1/e <= x.

    There are two branches, joining at -1/e = x.
    The lower branch extends from -1/e <= x < 0
    The upper branch extends from -1/e <= x

    The function is also known as the "Omega" function.

    In Mathematica, the function can be evaluated by:
      W = ProductLog [ X ]

    In MATLAB,
      W = lambertw ( b, x )
    where b = -1 for lower branch, 0 for upper branch.

    In Python,
      W = scipy.special.lambertw ( x, b )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    20 June 2023

  Author:

    John Burkardt

  Reference:

    Brian Hayes,
    "Why W?",
    The American Scientist,
    Volume 93, March-April 2005, pages 104-108.

    Eric Weisstein,
    "Lambert's W-Function",
    CRC Concise Encyclopedia of Mathematics,
    CRC Press, 1998.

    Stephen Wolfram,
    The Mathematica Book,
    Fourth Edition,
    Cambridge University Press, 1999,
    ISBN: 0-521-64314-7,
    LC: QA76.95.W65.

  Input:

    int *n_data: the user sets N_DATA to 0 before the first call.

  Output:

    int *n_data: the routine increments N_DATA by 1, and
    returns the corresponding data; when there is no more data, the
    output value of N_DATA will be 0 again.

    double *x: the argument of the function.

    double *fx: the value of the function.

    int b: -1 (lower branch) or 0 (upper branch).
*/
{
# define N_MAX 41

  static int branch_vec[N_MAX] = {
    -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
     0 };

  static double fx_vec[N_MAX] = {
   -4.889720169867429,
   -3.994308347002122,
   -3.439216483280204,
   -3.022313245324657,
   -2.678346990016661,
   -2.376421342062887,
   -2.097349210703492,
   -1.824388309032984,
   -1.531811608389612,
   -1.000000000000000,
   -0.608341284733432,
   -0.471671909743522,
   -0.374493134019498,
   -0.297083462446424,
   -0.231960952986534,
   -0.175356500529299,
   -0.125066982982524,
   -0.079678160511477,
   -0.038221241746799,
    0.0000000000000000E+00,
    0.3517337112491958E+00,
    0.5671432904097839E+00,
    0.7258613577662263E+00,
    0.8526055020137255E+00,
    0.9585863567287029E+00,
    0.1000000000000000E+01,
    0.1049908894964040E+01,
    0.1130289326974136E+01,
    0.1202167873197043E+01,
    0.1267237814307435E+01,
    0.1326724665242200E+01,
    0.1381545379445041E+01,
    0.1432404775898300E+01,
    0.1479856830173851E+01,
    0.1524345204984144E+01,
    0.1566230953782388E+01,
    0.1605811996320178E+01,
    0.1745528002740699E+01,
    0.3385630140290050E+01,
    0.5249602852401596E+01,
    0.1138335808614005E+02 };

  static double x_vec[N_MAX] = {
   -0.036787944117144,
   -0.073575888234288,
   -0.110363832351433,
   -0.147151776468577,
   -0.183939720585721,
   -0.220727664702865,
   -0.257515608820010,
   -0.294303552937154,
   -0.331091497054298,
   -0.367879441171442,
   -0.331091497054298,
   -0.294303552937154,
   -0.257515608820010,
   -0.220727664702865,
   -0.183939720585721,
   -0.147151776468577,
   -0.110363832351433,
   -0.073575888234288,
   -0.036787944117144,
    0.0000000000000000E+00,
    0.5000000000000000E+00,
    0.1000000000000000E+01,
    0.1500000000000000E+01,
    0.2000000000000000E+01,
    0.2500000000000000E+01,
    0.2718281828459045E+01,
    0.3000000000000000E+01,
    0.3500000000000000E+01,
    0.4000000000000000E+01,
    0.4500000000000000E+01,
    0.5000000000000000E+01,
    0.5500000000000000E+01,
    0.6000000000000000E+01,
    0.6500000000000000E+01,
    0.7000000000000000E+01,
    0.7500000000000000E+01,
    0.8000000000000000E+01,
    0.1000000000000000E+02,
    0.1000000000000000E+03,
    0.1000000000000000E+04,
    0.1000000000000000E+07 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
    *b = 0;
  }
  else
  {
    *x  = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
    *b = branch_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
/******************************************************************************/

void timestamp ( )

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

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}