lamv Subroutine

subroutine lamv ( v, x, vm, vl, dl )

  !*****************************************************************************80
  !
  !! LAMV computes lambda functions and derivatives of arbitrary order.
  !
  !  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:
  !
  !    31 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, real ( kind = 8 ) V, the order.
  !
  !    Input, real ( kind = 8 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) VM, the highest order computed.
  !
  !    Output, real ( kind = 8 ) VL(0:*), DL(0:*), the Lambda function and 
  !    derivative, of orders N+V0.
  !
  implicit none

  real ( kind = 8 ) v

  real ( kind = 8 ) a0
  real ( kind = 8 ) bjv0
  real ( kind = 8 ) bjv1
  real ( kind = 8 ) bk
  real ( kind = 8 ) ck
  real ( kind = 8 ) cs
  real ( kind = 8 ) dl(0:int(v))
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  real ( kind = 8 ) fac
  real ( kind = 8 ) ga
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) n
  real ( kind = 8 ) pi
  real ( kind = 8 ) px
  real ( kind = 8 ) qx
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) rc
  real ( kind = 8 ) rp
  real ( kind = 8 ) rp2
  real ( kind = 8 ) rq
  real ( kind = 8 ) sk
  real ( kind = 8 ) uk
  real ( kind = 8 ) v0
  real ( kind = 8 ) vk
  real ( kind = 8 ) vl(0:int(v))
  real ( kind = 8 ) vm
  real ( kind = 8 ) vv
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) xk

  pi = 3.141592653589793D+00
  rp2 = 0.63661977236758D+00
  x = abs ( x )
  x2 = x * x
  n = int ( v )
  v0 = v - n
  vm = v

  if ( x <= 12.0D+00 ) then

     do k = 0, n

        vk = v0 + k
        bk = 1.0D+00
        r = 1.0D+00

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

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

        dl(k) = - 0.5D+00 * x / ( vk + 1.0D+00 ) * uk

     end do

     return

  end if

  if ( x < 35.0D+00 ) then
     k0 = 11
  else if ( x < 50.0D+00 ) then
     k0 = 10
  else
     k0 = 8
  end if

  do j = 0, 1
     vv = 4.0D+00 * ( j + v0 ) * ( j + v0 )
     px = 1.0D+00
     rp = 1.0D+00
     do k = 1, k0
        rp = - 0.78125D-02 * rp * ( vv - ( 4.0D+00 * k - 3.0D+00 ) ** 2 ) &
             * ( vv - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
             / ( k * ( 2.0 * k - 1.0D+00 ) * x2 )
        px = px + rp
     end do
     qx = 1.0D+00
     rq = 1.0D+00
     do k = 1, k0
        rq = - 0.78125D-02 * rq * ( vv - ( 4.0D+00 * k - 1.0D+00 ) ** 2 ) &
             * ( vv - ( 4.0D+00 * k + 1.0D+00 ) ** 2 ) &
             / ( k * ( 2.0D+00 * k + 1.0D+00 ) * x2 )
        qx = qx + rq
     end do
     qx = 0.125D+00 * ( vv - 1.0D+00 ) * qx / x
     xk = x - ( 0.5D+00 * ( j + v0 ) + 0.25D+00 ) * pi
     a0 = sqrt ( rp2 / x )
     ck = cos ( xk )
     sk = sin ( xk )
     if ( j == 0 ) then
        bjv0 = a0 * ( px * ck - qx * sk )
     else
        bjv1 = a0 * ( px * ck - qx * sk )
     end if
  end do

  if ( v0 == 0.0D+00 ) then
     ga = 1.0D+00
  else
     call gam0 ( v0, ga )
     ga = v0 * ga
  end if

  fac = ( 2.0D+00 / x ) ** v0 * ga
  vl(0) = bjv0
  dl(0) = - bjv1 + v0 / x * bjv0
  vl(1) = bjv1
  dl(1) = bjv0 - ( 1.0D+00 + v0 ) / x * bjv1
  r0 = 2.0D+00 * ( 1.0D+00 + v0 ) / x

  if ( n <= 1 ) then
     vl(0) = fac * vl(0)
     dl(0) = fac * dl(0) - v0 / x * vl(0)
     vl(1) = fac * r0 * vl(1)
     dl(1) = fac * r0 * dl(1) - ( 1.0D+00 + v0 ) / x * vl(1)
     return
  end if

  if ( 2 <= n .and. n <= int ( 0.9D+00 * x ) ) then

     f0 = bjv0
     f1 = bjv1
     do k = 2, n
        f = 2.0D+00 * ( k + v0 - 1.0D+00 ) / x * f1 - f0
        f0 = f1
        f1 = f
        vl(k) = f
     end do

  else if ( 2 <= n ) then

     m = msta1 ( x, 200 )
     if ( m < n ) then
        n = m
     else
        m = msta2 ( x, n, 15 )
     end if
     f2 = 0.0D+00
     f1 = 1.0D-100
     do k = m, 0, -1
        f = 2.0D+00 * ( v0 + k + 1.0D+00 ) / x * f1 - f2
        if ( k <= n ) then
           vl(k) = f
        end if
        f2 = f1
        f1 = f
     end do

     if ( abs ( bjv0 ) <= abs ( bjv1 ) ) then
        cs = bjv1 / f2
     else
        cs = bjv0 / f
     end if

     do k = 0, n
        vl(k) = cs * vl(k)
     end do

  end if

  vl(0) = fac * vl(0)
  do j = 1, n
     rc = fac * r0
     vl(j) = rc * vl(j)
     dl(j-1) = - 0.5D+00 * x / ( j + v0 ) * vl(j)
     r0 = 2.0D+00 * ( j + v0 + 1 ) / x * r0
  end do
  dl(n) = 2.0D+00 * ( v0 + n ) * ( vl(n-1) - vl(n) ) / x
  vm = n + v0

  return
end subroutine lamv