lamn Subroutine

subroutine lamn ( n, x, nm, bl, dl )

  !*****************************************************************************80
  !
  !! LAMN computes lambda functions and derivatives.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    14 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) BL(0:N), DL(0:N), the
  !    value of the lambda function and its derivative of orders 0 through N.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) bg
  real ( kind = 8 ) bk
  real ( kind = 8 ) bl(0:n)
  real ( kind = 8 ) bs
  real ( kind = 8 ) dl(0:n)
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) i
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) uk
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

  nm = n

  if ( abs ( x ) < 1.0D-100 ) then
     do k = 0, n
        bl(k) = 0.0D+00
        dl(k) = 0.0D+00
     end do
     bl(0) = 1.0D+00
     dl(1) = 0.5D+00
     return
  end if

  if ( x <= 12.0D+00 ) then

     x2 = x * x

     do k = 0, n
        bk = 1.0D+00
        r = 1.0D+00
        do i = 1, 50
           r = -0.25D+00 * r * x2 / ( i * ( i + k ) )
           bk = bk + r
           if ( abs ( r ) < abs ( bk ) * 1.0D-15 ) then
              exit
           end if
        end do

        bl(k) = bk
        if ( 1 <= k ) then
           dl(k-1) = - 0.5D+00 * x / k * bk
        end if

     end do

     uk = 1.0D+00
     r = 1.0D+00
     do i = 1, 50
        r = -0.25D+00 * r * x2 / ( i * ( i + n + 1.0D+00 ) )
        uk = uk + r
        if ( abs ( r ) < abs ( uk ) * 1.0D-15 ) then
           exit
        end if
     end do

     dl(n) = -0.5D+00 * x / ( n + 1.0D+00 ) * uk
     return

  end if

  if ( n == 0 ) then
     nm = 1
  end if

  m = msta1 ( x, 200 )

  if ( m < nm ) then
     nm = m
  else
     m = msta2 ( x, nm, 15 )
  end if

  bs = 0.0D+00
  f0 = 0.0D+00
  f1 = 1.0D-100
  do k = m, 0, -1
     f = 2.0D+00 * ( k + 1.0D+00 ) * f1 / x - f0
     if ( k <= nm ) then
        bl(k) = f
     end if
     if ( k == 2 * int ( k / 2 ) ) then
        bs = bs + 2.0D+00 * f
     end if
     f0 = f1
     f1 = f
  end do

  bg = bs - f
  do k = 0, nm
     bl(k) = bl(k) / bg
  end do

  r0 = 1.0D+00
  do k = 1, nm
     r0 = 2.0D+00 * r0 * k / x
     bl(k) = r0 * bl(k)
  end do

  dl(0) = -0.5D+00 * x * bl(1)
  do k = 1, nm
     dl(k) = 2.0D+00 * k / x * ( bl(k-1) - bl(k) )
  end do

  return
end subroutine lamn