jy01b Subroutine

subroutine jy01b ( x, bj0, dj0, bj1, dj1, by0, dy0, by1, dy1 )

  !*****************************************************************************80
  !
  !! JY01B computes Bessel functions J0(x), J1(x), Y0(x), Y1(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:
  !
  !    02 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 ) BJ0, DJ0, BJ1, DJ1, BY0, DY0, BY1, DY1,
  !    the values of J0(x), J0'(x), J1(x), J1'(x), Y0(x), Y0'(x), Y1(x), Y1'(x).
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) bj0
  real ( kind = 8 ) bj1
  real ( kind = 8 ) by0
  real ( kind = 8 ) by1
  real ( kind = 8 ) dj0
  real ( kind = 8 ) dj1
  real ( kind = 8 ) dy0
  real ( kind = 8 ) dy1
  real ( kind = 8 ) p0
  real ( kind = 8 ) p1
  real ( kind = 8 ) pi
  real ( kind = 8 ) q0
  real ( kind = 8 ) q1
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) ta0
  real ( kind = 8 ) ta1
  real ( kind = 8 ) x

  pi = 3.141592653589793D+00

  if ( x == 0.0D+00 ) then

     bj0 = 1.0D+00
     bj1 = 0.0D+00
     dj0 = 0.0D+00
     dj1 = 0.5D+00
     by0 = -1.0D+300
     by1 = -1.0D+300
     dy0 = 1.0D+300
     dy1 = 1.0D+300
     return

  else if ( x <= 4.0D+00 ) then

     t = x / 4.0D+00
     t2 = t * t

     bj0 = (((((( &
          - 0.5014415D-03 * t2 &
          + 0.76771853D-02 ) * t2 &
          - 0.0709253492D+00 ) * t2 &
          + 0.4443584263D+00 ) * t2 &
          - 1.7777560599D+00 ) * t2 &
          + 3.9999973021D+00 ) * t2 &
          - 3.9999998721D+00 ) * t2 &
          + 1.0D+00

     bj1 = t * ((((((( &
          - 0.1289769D-03 * t2 &
          + 0.22069155D-02 ) * t2 &
          - 0.0236616773D+00 ) * t2 &
          + 0.1777582922D+00 ) * t2 &
          - 0.8888839649D+00 ) * t2 &
          + 2.6666660544D+00 ) * t2 &
          - 3.9999999710D+00 ) * t2 &
          + 1.9999999998D+00 )

     by0 = ((((((( &
          - 0.567433D-04 * t2 &
          + 0.859977D-03 ) * t2 &
          - 0.94855882D-02 ) * t2 &
          + 0.0772975809D+00 ) * t2 &
          - 0.4261737419D+00 ) * t2 &
          + 1.4216421221D+00 ) * t2 &
          - 2.3498519931D+00 ) * t2 &
          + 1.0766115157D+00 ) * t2 &
          + 0.3674669052D+00

     by0 = 2.0D+00 / pi * log ( x / 2.0D+00 ) * bj0 + by0

     by1 = (((((((( &
          0.6535773D-03 * t2 &
          - 0.0108175626D+00 ) * t2 &
          + 0.107657606D+00 ) * t2 &
          - 0.7268945577D+00 ) * t2 &
          + 3.1261399273D+00 ) * t2 &
          - 7.3980241381D+00 ) * t2 &
          + 6.8529236342D+00 ) * t2 &
          + 0.3932562018D+00 ) * t2 &
          - 0.6366197726D+00 ) / x

     by1 = 2.0D+00 / pi * log ( x / 2.0D+00 ) * bj1 + by1

  else

     t = 4.0D+00 / x
     t2 = t * t
     a0 = sqrt ( 2.0D+00 / ( pi * x ) )

     p0 = (((( &
          - 0.9285D-05 * t2 &
          + 0.43506D-04 ) * t2 &
          - 0.122226D-03 ) * t2 &
          + 0.434725D-03 ) * t2 &
          - 0.4394275D-02 ) * t2 &
          + 0.999999997D+00

     q0 = t * ((((( &
          0.8099D-05 * t2 &
          - 0.35614D-04 ) * t2 &
          + 0.85844D-04 ) * t2 &
          - 0.218024D-03 ) * t2 &
          + 0.1144106D-02 ) * t2 &
          - 0.031249995D+00 )

     ta0 = x - 0.25D+00 * pi
     bj0 = a0 * ( p0 * cos ( ta0 ) - q0 * sin ( ta0 ) )
     by0 = a0 * ( p0 * sin ( ta0 ) + q0 * cos ( ta0 ) )

     p1 = (((( &
          0.10632D-04 * t2 &
          - 0.50363D-04 ) * t2 &
          + 0.145575D-03 ) * t2 &
          - 0.559487D-03 ) * t2 &
          + 0.7323931D-02 ) * t2 &
          + 1.000000004D+00

     q1 = t * ((((( &
          - 0.9173D-05      * t2 &
          + 0.40658D-04 )   * t2 &
          - 0.99941D-04 )   * t2 &
          + 0.266891D-03 )  * t2 &
          - 0.1601836D-02 ) * t2 &
          + 0.093749994D+00 )

     ta1 = x - 0.75D+00 * pi
     bj1 = a0 * ( p1 * cos ( ta1 ) - q1 * sin ( ta1 ) )
     by1 = a0 * ( p1 * sin ( ta1 ) + q1 * cos ( ta1 ) )

  end if

  dj0 = - bj1
  dj1 = bj0 - bj1 / x
  dy0 = - by1
  dy1 = by0 - by1 / x

  return
end subroutine jy01b