function W = lambert_w ( X, B, l )

%*****************************************************************************80
%
%% lambert_w() approximates the Lambert W function.
%
%  Discussion:
%
%    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:
%
%    19 June 2023
%
%  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:
%
%    real X(*): the argument.
%
%    integer B(*): indicates the desired branch.
%    * 0, the upper branch, -exp(-1) <= X.
%    * nonzero, the lower branch, -exp(-1) < X < 0.
%
%    integer 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:
%
%    real W(*): the approximate value of W(X).
%
  W = NaN ( size ( X ) );

  nbits = 52;
  niter = 1;

  em = - exp ( -1.0 );
  em9 = - exp ( -9.0 );
  c13 = 1.0 / 3.0;
  c23 = 2.0 * c13;
  em2 = 2.0 / em;
  d12 = - em2;
  tb = 0.5 ^ nbits;
  tb2 = sqrt ( tb );
  x0 = tb ^ ( 1.0 / 6.0 ) * 0.5;
  x1 = ( 1.0 - 17.0 * 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;
%
%  Code is not vectorized.
%
  [ rows, cols ] = size ( X );

  for i = 1 : rows
    for j = 1 : cols
%
%  Check and adjust the argument.
%
      if ( l == 1 )

        delx = X(i,j);

        if ( delx < 0.0 )
          continue
        end

        xx = X(i,j) + em;

      else

        if ( X(i,j) < em )
          continue
        elseif ( X(i,j) == em )
          W(i,j) = -1.0;
          continue
        end

        xx = X(i,j);
        delx = xx - em;

      end
%
%  Calculations for the upper branch.
%
      if ( size ( B ) == 1 )
        nb = B;
      else
        nb = B(i,j);
      end

      if ( nb == 0 )

        if ( abs ( xx ) <= x0 )
          W(i,j) = xx / ( 1.0 + xx / ( 1.0 + xx ...
            / ( 2.0 + xx / ( 0.6 + 0.34 * xx ))));
          continue
        elseif ( xx <= x1 )
          reta = sqrt ( d12 * delx );
          W(i,j) = reta / ( 1.0 + reta / ( 3.0 + reta / ( reta ...
            / ( an4 + reta / ( reta * an6 + an5 ) ) + an3 ) ) ) - 1.0;
          continue
        elseif ( xx <= 20.0 )
          reta = s2 * sqrt ( 1.0 - xx / em );
          an2 = 4.612634277343749 * sqrt ( sqrt ( reta + 1.09556884765625 ));
          W(i,j) = reta / ( 1.0 + reta / ( 3.0 + ( s21 * an2 ...
            + s22 ) * reta / ( s23 * ( an2 + reta )))) - 1.0;
        else
          zl = log ( xx );
          W(i,j) = log ( xx / log ( xx ...
            / zl ^ exp ( -1.124491989777808 / ...
           ( 0.4225028202459761 + zl ))));
        end
%
%  Calculations for the lower branch.
%
      else

        if ( 0.0 <= xx )
          continue
        elseif ( xx <= x1 )
          reta = sqrt ( d12 * delx );
          W(i,j) = reta / ( reta / ( 3.0 + reta / ( reta / ( an4 ...
            + reta / ( reta * an6 - an5 ) ) - an3 ) ) - 1.0 ) - 1.0;
          continue
        elseif ( xx <= em9 )
          zl = log ( -xx );
          t = -1.0 - zl;
          ts = sqrt ( t );
          W(i,j) = zl - ( 2.0 * ts ) / ( s2 + ( c13 - t ...
            / ( 270.0 + ts * 127.0471381349219 )) * ts );
        else
          zl = log ( - xx );
          eta = 2.0 - em2 * xx;
          W(i,j) = log ( xx / log ( - xx / ( ( 1.0 ...
            - 0.5043921323068457 * ( zl + 1.0 ) ) ...
            * ( sqrt ( eta ) + eta / 3.0 ) + 1.0 )));
        end

      end

      for iter = 1 : niter
        zn = log ( xx / W(i,j) ) - W(i,j);
        temp = 1.0 + W(i,j);
        temp2 = temp + c23 * zn;
        temp2 = 2.0 * temp * temp2;
        W(i,j) = W(i,j) * ( 1.0 + ( zn / temp ) * ( temp2 - zn ) ...
          / ( temp2 - 2.0 * zn ) );
      end

    end
  end

  return
end
