klvna Subroutine

subroutine klvna ( x, ber, bei, ger, gei, der, dei, her, hei )

  !*****************************************************************************80
  !
  !! KLVNA: Kelvin functions ber(x), bei(x), ker(x), and kei(x), 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:
  !
  !    03 August 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 ) X, the argument.
  !
  !    Output, real ( kind = 8 ) BER, BEI, GER, GEI, DER, DEI, HER, HEI, 
  !    the values of ber x, bei x, ker x, kei x, ber'x, bei'x, ker'x, kei'x.
  !
  implicit none

  real ( kind = 8 ) bei
  real ( kind = 8 ) ber
  real ( kind = 8 ) cn0
  real ( kind = 8 ) cp0
  real ( kind = 8 ) cs
  real ( kind = 8 ) dei
  real ( kind = 8 ) der
  real ( kind = 8 ) el
  real ( kind = 8 ) eps
  real ( kind = 8 ) fac
  real ( kind = 8 ) gei
  real ( kind = 8 ) ger
  real ( kind = 8 ) gs
  real ( kind = 8 ) hei
  real ( kind = 8 ) her
  integer ( kind = 4 ) k
  integer ( kind = 4 ) km
  integer ( kind = 4 ) m
  real ( kind = 8 ) pi
  real ( kind = 8 ) pn0
  real ( kind = 8 ) pn1
  real ( kind = 8 ) pp0
  real ( kind = 8 ) pp1
  real ( kind = 8 ) qn0
  real ( kind = 8 ) qn1
  real ( kind = 8 ) qp0
  real ( kind = 8 ) qp1
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r1
  real ( kind = 8 ) rc
  real ( kind = 8 ) rs
  real ( kind = 8 ) sn0
  real ( kind = 8 ) sp0
  real ( kind = 8 ) ss
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) x4
  real ( kind = 8 ) xc1
  real ( kind = 8 ) xc2
  real ( kind = 8 ) xd
  real ( kind = 8 ) xe1
  real ( kind = 8 ) xe2
  real ( kind = 8 ) xt

  pi = 3.141592653589793D+00
  el = 0.5772156649015329D+00
  eps = 1.0D-15

  if ( x == 0.0D+00 ) then
     ber = 1.0D+00
     bei = 0.0D+00
     ger = 1.0D+300
     gei = -0.25D+00 * pi
     der = 0.0D+00
     dei = 0.0D+00
     her = -1.0D+300
     hei = 0.0D+00
     return
  end if

  x2 = 0.25D+00 * x * x
  x4 = x2 * x2

  if ( abs ( x ) < 10.0D+00 ) then

     ber = 1.0D+00
     r = 1.0D+00
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2.0D+00 * m - 1.0D+00 ) ** 2 * x4
        ber = ber + r
        if ( abs ( r ) < abs ( ber ) * eps ) then
           exit
        end if
     end do

     bei = x2
     r = x2
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2.0D+00 * m + 1.0D+00 ) ** 2 * x4
        bei = bei + r
        if ( abs ( r ) < abs ( bei ) * eps ) then
           exit
        end if
     end do

     ger = - ( log ( x / 2.0D+00 ) + el ) * ber + 0.25D+00 * pi * bei
     r = 1.0D+00
     gs = 0.0D+00
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2.0D+00 * m - 1.0D+00 ) ** 2 * x4
        gs = gs + 1.0D+00 / ( 2.0D+00 * m - 1.0D+00 ) + 1.0D+00 / ( 2.0D+00 * m )
        ger = ger + r * gs
        if ( abs ( r * gs ) < abs ( ger ) * eps ) then
           exit
        end if
     end do

     gei = x2 - ( log ( x / 2.0D+00 ) + el ) * bei - 0.25D+00 * pi * ber
     r = x2
     gs = 1.0D+00
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2.0D+00 * m + 1.0D+00 ) ** 2 * x4
        gs = gs + 1.0D+00 / ( 2.0D+00 * m ) + 1.0D+00 / ( 2.0D+00 * m + 1.0D+00 )
        gei = gei + r * gs
        if ( abs ( r * gs ) < abs ( gei ) * eps ) then
           exit
        end if
     end do

     der = -0.25D+00 * x * x2
     r = der
     do m = 1, 60
        r = -0.25D+00 * r / m / ( m + 1.0D+00 ) &
             / ( 2.0D+00 * m + 1.0D+00 ) ** 2 * x4
        der = der + r
        if ( abs ( r ) < abs ( der ) * eps ) then
           exit
        end if
     end do

     dei = 0.5D+00 * x
     r = dei
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2.0D+00 * m - 1.0D+00 ) &
             / ( 2.0D+00 * m + 1.0D+00 ) * x4
        dei = dei + r
        if ( abs ( r ) < abs ( dei ) * eps ) then
           exit
        end if
     end do

     r = -0.25D+00 * x * x2
     gs = 1.5D+00
     her = 1.5D+00 * r - ber / x &
          - ( log ( x / 2.0D+00 ) + el ) * der + 0.25D+00 * pi * dei
     do m = 1, 60
        r = -0.25D+00 * r / m / ( m + 1.0D+00 ) &
             / ( 2.0D+00 * m + 1.0D+00 ) ** 2 * x4
        gs = gs + 1.0D+00 / ( 2 * m + 1.0D+00 ) + 1.0D+00 &
             / ( 2 * m + 2.0D+00 )
        her = her + r * gs
        if ( abs ( r * gs ) < abs ( her ) * eps ) then
           exit
        end if
     end do

     r = 0.5D+00 * x
     gs = 1.0D+00
     hei = 0.5D+00 * x - bei / x &
          - ( log ( x / 2.0D+00 ) + el ) * dei - 0.25D+00 * pi * der
     do m = 1, 60
        r = -0.25D+00 * r / ( m * m ) / ( 2 * m - 1.0D+00 ) &
             / ( 2 * m + 1.0D+00 ) * x4
        gs = gs + 1.0D+00 / ( 2.0D+00 * m ) + 1.0D+00 &
             / ( 2 * m + 1.0D+00 )
        hei = hei + r * gs
        if ( abs ( r * gs ) < abs ( hei ) * eps ) then 
           return
        end if
     end do

  else

     pp0 = 1.0D+00
     pn0 = 1.0D+00
     qp0 = 0.0D+00
     qn0 = 0.0D+00
     r0 = 1.0D+00

     if ( abs ( x ) < 40.0D+00 ) then
        km = 18
     else
        km = 10
     end if

     fac = 1.0D+00
     do k = 1, km
        fac = -fac
        xt = 0.25D+00 * k * pi - int ( 0.125D+00 * k ) * 2.0D+00 * pi
        cs = cos ( xt )
        ss = sin ( xt )
        r0 = 0.125D+00 * r0 * ( 2.0D+00 * k - 1.0D+00 ) ** 2 / k / x
        rc = r0 * cs
        rs = r0 * ss
        pp0 = pp0 + rc
        pn0 = pn0 + fac * rc
        qp0 = qp0 + rs
        qn0 = qn0 + fac * rs
     end do

     xd = x / sqrt (2.0D+00 )
     xe1 = exp ( xd )
     xe2 = exp ( - xd )
     xc1 = 1.0D+00 / sqrt ( 2.0D+00 * pi * x )
     xc2 = sqrt ( 0.5D+00 * pi / x )
     cp0 = cos ( xd + 0.125D+00 * pi )
     cn0 = cos ( xd - 0.125D+00 * pi )
     sp0 = sin ( xd + 0.125D+00 * pi )
     sn0 = sin ( xd - 0.125D+00 * pi )
     ger = xc2 * xe2 * (  pn0 * cp0 - qn0 * sp0 )
     gei = xc2 * xe2 * ( -pn0 * sp0 - qn0 * cp0 )
     ber = xc1 * xe1 * (  pp0 * cn0 + qp0 * sn0 ) - gei / pi
     bei = xc1 * xe1 * (  pp0 * sn0 - qp0 * cn0 ) + ger / pi
     pp1 = 1.0D+00
     pn1 = 1.0D+00
     qp1 = 0.0D+00
     qn1 = 0.0D+00
     r1 = 1.0D+00
     fac = 1.0D+00

     do k = 1, km
        fac = -fac
        xt = 0.25D+00 * k * pi - int ( 0.125D+00 * k ) * 2.0D+00 * pi
        cs = cos ( xt )
        ss = sin ( xt )
        r1 = 0.125D+00 * r1 &
             * ( 4.0D+00 - ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) / k / x
        rc = r1 * cs
        rs = r1 * ss
        pp1 = pp1 + fac * rc
        pn1 = pn1 + rc
        qp1 = qp1 + fac * rs
        qn1 = qn1 + rs
     end do

     her = xc2 * xe2 * ( - pn1 * cn0 + qn1 * sn0 )
     hei = xc2 * xe2 * (   pn1 * sn0 + qn1 * cn0 )
     der = xc1 * xe1 * (   pp1 * cp0 + qp1 * sp0 ) - hei / pi
     dei = xc1 * xe1 * (   pp1 * sp0 - qp1 * cp0 ) + her / pi

  end if

  return
end subroutine klvna