Lambert Function

  1. module_numerics.f90
  2. numerics
  3. Lambert

public function Lambert(z)

Function to calculate Lambert's W-function.

Algorithm found at http://en.citizendium.org/wiki/Lambert_W_function#Numerical_calculation

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in) :: z

The input value.

Return Value real(kind=wp)

The return value Lambert(z).


Calls

proc~~lambert~~CallsGraph proc~lambert Lambert interface~eps eps proc~lambert->interface~eps proc~eps_prec_s eps_prec_s interface~eps->proc~eps_prec_s proc~eps_prec_q eps_prec_q interface~eps->proc~eps_prec_q proc~eps_prec_d eps_prec_d interface~eps->proc~eps_prec_d

Called by

proc~~lambert~~CalledByGraph proc~lambert Lambert proc~invert_tortoise invert_tortoise proc~invert_tortoise->proc~lambert proc~tdc_set_coefficients tdc_set_coefficients proc~tdc_set_coefficients->proc~invert_tortoise proc~scal_schw_init scal_schw_init proc~scal_schw_init->proc~invert_tortoise proc~output_coords output_coords proc~output_coords->proc~invert_tortoise interface~output_coords output_coords interface~output_coords->proc~output_coords interface~scal_schw_init scal_schw_init interface~scal_schw_init->proc~scal_schw_init interface~tdc_set_coefficients tdc_set_coefficients interface~tdc_set_coefficients->proc~tdc_set_coefficients

Contents

Source Code

Lambert

Source Code

  function Lambert ( z )
  !! Function to calculate Lambert's W-function.
  !!
  !! Algorithm found at
  !! http://en.citizendium.org/wiki/Lambert_W_function#Numerical_calculation

    implicit none

    real(wp), intent(in) :: z
    !! The input value.
    real(wp) :: Lambert
    !! The return value Lambert(z).
    real(wp) :: wcurrent, wnew, expw, diff
    integer(ip) :: iter
!    real(wp), parameter :: eps = 1.0e-12_wp

    iter = 0
    wcurrent = 1.0_wp

    loop: do
      iter = iter+1
      expw = exp(wcurrent)
      diff = wcurrent*expw - z
      wnew = wcurrent - diff / ( expw*(wcurrent+1) - &
             ((wcurrent+2)*diff/(2*wcurrent+2)))
      if ( abs(wnew - wcurrent) < eps(1.0_wp) ) exit loop
      if (iter>maxiter) stop 'Function Lambert failed to converge. Aborting'
      wcurrent = wnew
    end do loop
    Lambert = wnew
  end function Lambert