klvnb Subroutine

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

  !*****************************************************************************80
  !
  !! KLVNB: 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 ) csn
  real ( kind = 8 ) csp
  real ( kind = 8 ) dei
  real ( kind = 8 ) der
  real ( kind = 8 ) fxi
  real ( kind = 8 ) fxr
  real ( kind = 8 ) gei
  real ( kind = 8 ) ger
  real ( kind = 8 ) hei
  real ( kind = 8 ) her
  integer ( kind = 4 ) l
  real ( kind = 8 ) pi
  real ( kind = 8 ) pni
  real ( kind = 8 ) pnr
  real ( kind = 8 ) ppi
  real ( kind = 8 ) ppr
  real ( kind = 8 ) ssn
  real ( kind = 8 ) ssp
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) tni
  real ( kind = 8 ) tnr
  real ( kind = 8 ) tpi
  real ( kind = 8 ) tpr
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) x
  real ( kind = 8 ) yc1
  real ( kind = 8 ) yc2
  real ( kind = 8 ) yci
  real ( kind = 8 ) ye1
  real ( kind = 8 ) ye2
  real ( kind = 8 ) yei
  real ( kind = 8 ) yd

  pi = 3.141592653589793D+00

  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

  else if ( x < 8.0D+00 ) then

     t = x / 8.0D+00
     t2 = t * t
     u = t2 * t2

     ber = (((((( &
          - 0.901D-05 * u &
          + 0.122552D-02 ) * u &
          - 0.08349609D+00 ) * u &
          + 2.64191397D+00 ) * u &
          - 32.36345652D+00 ) * u &
          + 113.77777774D+00 ) * u &
          - 64.0D+00 ) * u &
          + 1.0D+00

     bei = t * t * (((((( &
          0.11346D-03 * u &
          - 0.01103667D+00 ) * u &
          + 0.52185615D+00 ) * u &
          - 10.56765779D+00 ) * u &
          + 72.81777742D+00 ) * u &
          - 113.77777774D+00 ) * u &
          + 16.0D+00 )

     ger = (((((( &
          - 0.2458D-04 * u &
          + 0.309699D-02 ) * u &
          - 0.19636347D+00 ) * u &
          + 5.65539121D+00 ) * u &
          - 60.60977451D+00 ) * u &
          + 171.36272133D+00 ) * u &
          - 59.05819744D+00 ) * u &
          - 0.57721566D+00

     ger = ger - log ( 0.5D+00 * x ) * ber + 0.25D+00 * pi * bei

     gei = t2 * (((((( &
          0.29532D-03 * u &
          - 0.02695875D+00 ) * u &
          + 1.17509064D+00 ) * u &
          - 21.30060904D+00 ) * u &
          + 124.2356965D+00 ) * u &
          - 142.91827687D+00 ) * u &
          + 6.76454936D+00 )

     gei = gei - log ( 0.5D+00 * x ) * bei - 0.25D+00 * pi * ber

     der = x * t2 * (((((( &
          - 0.394D-05 * u &
          + 0.45957D-03 ) * u &
          - 0.02609253D+00 ) * u &
          + 0.66047849D+00 ) * u &
          - 6.0681481D+00 ) * u &
          + 14.22222222D+00 ) * u &
          - 4.0D+00 )

     dei = x * (((((( &
          0.4609D-04 * u &
          - 0.379386D-02 ) * u &
          + 0.14677204D+00 ) * u &
          - 2.31167514D+00 ) * u &
          + 11.37777772D+00 ) * u &
          - 10.66666666D+00 ) * u &
          + 0.5D+00 ) 

     her = x * t2 * (((((( &
          - 0.1075D-04 * u &
          + 0.116137D-02 ) * u &
          - 0.06136358D+00 ) * u &
          + 1.4138478D+00 ) * u &
          - 11.36433272D+00 ) * u &
          + 21.42034017D+00 ) * u &
          - 3.69113734D+00 )

     her = her - log ( 0.5D+00 * x ) * der - ber / x  &
          + 0.25D+00 * pi * dei

     hei = x * (((((( &
          0.11997D-03 * u &
          - 0.926707D-02 ) * u &
          + 0.33049424D+00 ) * u &
          - 4.65950823D+00 ) * u &
          + 19.41182758D+00 ) * u &
          - 13.39858846D+00 ) * u &
          + 0.21139217D+00 )

     hei = hei - log ( 0.5D+00 * x ) * dei - bei / x  &
          - 0.25D+00 * pi * der

  else

     t = 8.0D+00 / x

     do l = 1, 2

        v = ( -1.0D+00 ) ** l * t

        tpr = (((( &
             0.6D-06 * v &
             - 0.34D-05 ) * v &
             - 0.252D-04 ) * v &
             - 0.906D-04 ) * v * v &
             + 0.0110486D+00 ) * v

        tpi = (((( &
             0.19D-05 * v &
             + 0.51D-05 ) * v * v &
             - 0.901D-04 ) * v &
             - 0.9765D-03 ) * v &
             - 0.0110485D+00 ) * v &
             - 0.3926991D+00

        if ( l == 1 ) then
           tnr = tpr
           tni = tpi
        end if

     end do

     yd = x / sqrt ( 2.0D+00 )
     ye1 = exp ( yd + tpr )
     ye2 = exp ( - yd + tnr )
     yc1 = 1.0D+00 / sqrt ( 2.0D+00 * pi * x )
     yc2 = sqrt ( pi / ( 2.0D+00 * x ) )
     csp = cos ( yd + tpi )
     ssp = sin ( yd + tpi )
     csn = cos ( - yd + tni )
     ssn = sin ( - yd + tni )
     ger = yc2 * ye2 * csn
     gei = yc2 * ye2 * ssn
     fxr = yc1 * ye1 * csp
     fxi = yc1 * ye1 * ssp
     ber = fxr - gei / pi
     bei = fxi + ger / pi

     do l = 1, 2

        v = ( -1.0D+00 ) ** l * t

        ppr = ((((( &
             0.16D-05 * v &
             + 0.117D-04 ) * v &
             + 0.346D-04 ) * v &
             + 0.5D-06 ) * v &
             - 0.13813D-02 ) * v &
             - 0.0625001D+00 ) * v &
             + 0.7071068D+00

        ppi = ((((( &
             - 0.32D-05 * v &
             - 0.24D-05 ) * v &
             + 0.338D-04 ) * v &
             + 0.2452D-03 ) * v &
             + 0.13811D-02 ) * v &
             - 0.1D-06 ) * v &
             + 0.7071068D+00

        if ( l == 1 ) then
           pnr = ppr
           pni = ppi
        end if

     end do

     her =     gei * pni - ger * pnr
     hei = - ( gei * pnr + ger * pni )
     der = fxr * ppr - fxi * ppi - hei / pi
     dei = fxi * ppr + fxr * ppi + her / pi

  end if

  return
end subroutine klvnb