ajyik Subroutine

subroutine ajyik ( x, vj1, vj2, vy1, vy2, vi1, vi2, vk1, vk2 )

  !*****************************************************************************80
  !
  !! AJYIK computes Bessel functions Jv(x), Yv(x), Iv(x), Kv(x).
  !
  !  Discussion: 
  !
  !    Compute Bessel functions Jv(x) and Yv(x), and modified Bessel functions 
  !    Iv(x) and Kv(x), and their derivatives with v = 1/3, 2/3.
  !
  !  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 ) X, the argument.  X should not be zero.
  !
  !    Output, real ( kind = 8 ) VJ1, VJ2, VY1, VY2, VI1, VI2, VK1, VK2,
  !    the values of J1/3(x), J2/3(x), Y1/3(x), Y2/3(x), I1/3(x), I2/3(x),
  !    K1/3(x), K2/3(x).
  !
  implicit none

  real ( kind = 8 ) a0
  real ( kind = 8 ) b0
  real ( kind = 8 ) c0
  real ( kind = 8 ) ck
  real ( kind = 8 ) gn
  real ( kind = 8 ) gn1
  real ( kind = 8 ) gn2
  real ( kind = 8 ) gp1
  real ( kind = 8 ) gp2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k0
  integer ( kind = 4 ) l
  real ( kind = 8 ) pi
  real ( kind = 8 ) pv1
  real ( kind = 8 ) pv2
  real ( kind = 8 ) px
  real ( kind = 8 ) qx
  real ( kind = 8 ) r
  real ( kind = 8 ) rp
  real ( kind = 8 ) rp2
  real ( kind = 8 ) rq
  real ( kind = 8 ) sk
  real ( kind = 8 ) sum
  real ( kind = 8 ) uj1
  real ( kind = 8 ) uj2
  real ( kind = 8 ) uu0
  real ( kind = 8 ) vi1
  real ( kind = 8 ) vi2
  real ( kind = 8 ) vil
  real ( kind = 8 ) vj1
  real ( kind = 8 ) vj2
  real ( kind = 8 ) vjl
  real ( kind = 8 ) vk1
  real ( kind = 8 ) vk2
  real ( kind = 8 ) vl
  real ( kind = 8 ) vsl
  real ( kind = 8 ) vv
  real ( kind = 8 ) vv0
  real ( kind = 8 ) vy1
  real ( kind = 8 ) vy2
  real ( kind = 8 ) x
  real ( kind = 8 ) x2
  real ( kind = 8 ) xk

  if ( x == 0.0D+00 ) then
     vj1 = 0.0D+00
     vj2 = 0.0D+00
     vy1 = -1.0D+300
     vy2 = 1.0D+300
     vi1 = 0.0D+00
     vi2 = 0.0D+00
     vk1 = -1.0D+300
     vk2 = -1.0D+300
     return
  end if

  pi = 3.141592653589793D+00
  rp2 = 0.63661977236758D+00
  gp1 = 0.892979511569249D+00
  gp2 = 0.902745292950934D+00
  gn1 = 1.3541179394264D+00
  gn2 = 2.678938534707747D+00
  vv0 = 0.444444444444444D+00
  uu0 = 1.1547005383793D+00
  x2 = x * x

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

  if ( x <= 12.0D+00 ) then

     do l = 1, 2
        vl = l / 3.0D+00
        vjl = 1.0D+00
        r = 1.0D+00
        do k = 1, 40
           r = -0.25D+00 * r * x2 / ( k * ( k + vl ) )
           vjl = vjl + r
           if ( abs ( r ) < 1.0D-15 ) then
              exit
           end if
        end do

        a0 = ( 0.5D+00 * x ) ** vl
        if ( l == 1 ) then
           vj1 = a0 / gp1 * vjl
        else
           vj2 = a0 / gp2 * vjl
        end if

     end do

  else

     do l = 1, 2

        vv = vv0 * l * l
        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.0D+00 * 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 * l / 3.0D+00 + 0.25D+00 ) * pi
        a0 = sqrt ( rp2 / x )
        ck = cos ( xk )
        sk = sin ( xk )
        if ( l == 1) then
           vj1 = a0 * ( px * ck - qx * sk )
           vy1 = a0 * ( px * sk + qx * ck )
        else
           vj2 = a0 * ( px * ck - qx * sk )
           vy2 = a0 * ( px * sk + qx * ck )
        end if

     end do

  end if

  if ( x <= 12.0D+00 ) then

     do l = 1, 2

        vl = l / 3.0D+00
        vjl = 1.0D+00
        r = 1.0D+00
        do k = 1, 40
           r = -0.25D+00 * r * x2 / ( k * ( k - vl ) )
           vjl = vjl + r
           if ( abs ( r ) < 1.0D-15 ) then
              exit
           end if
        end do

        b0 = ( 2.0D+00 / x ) ** vl
        if ( l == 1 ) then
           uj1 = b0 * vjl / gn1
        else
           uj2 = b0 * vjl / gn2
        end if

     end do

     pv1 = pi / 3.0D+00
     pv2 = pi / 1.5D+00
     vy1 = uu0 * ( vj1 * cos ( pv1 ) - uj1 )
     vy2 = uu0 * ( vj2 * cos ( pv2 ) - uj2 )

  end if

  if ( x <= 18.0D+00 ) then

     do l = 1, 2
        vl = l / 3.0D+00
        vil = 1.0D+00
        r = 1.0D+00
        do k = 1, 40
           r = 0.25D+00 * r * x2 / ( k * ( k + vl ) )
           vil = vil + r
           if ( abs ( r ) < 1.0D-15 ) then
              exit
           end if
        end do

        a0 = ( 0.5D+00 * x ) ** vl

        if ( l == 1 ) then
           vi1 = a0 / gp1 * vil
        else
           vi2 = a0 / gp2 * vil
        end if

     end do

  else

     c0 = exp ( x ) / sqrt ( 2.0D+00 * pi * x )

     do l = 1, 2
        vv = vv0 * l * l
        vsl = 1.0D+00
        r = 1.0D+00
        do k = 1, k0
           r = - 0.125D+00 * r &
                * ( vv - ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) / ( k * x )
           vsl = vsl + r
        end do
        if ( l == 1 ) then
           vi1 = c0 * vsl
        else
           vi2 = c0 * vsl
        end if
     end do

  end if

  if ( x <= 9.0D+00 ) then

     do l = 1, 2
        vl = l / 3.0D+00
        if ( l == 1 ) then
           gn = gn1
        else
           gn = gn2
        end if
        a0 = ( 2.0D+00 / x ) ** vl / gn
        sum = 1.0D+00
        r = 1.0D+00
        do k = 1, 60
           r = 0.25D+00 * r * x2 / ( k * ( k - vl ) )
           sum = sum + r
           if ( abs ( r ) < 1.0D-15 ) then
              exit
           end if
        end do

        if ( l == 1 ) then
           vk1 = 0.5D+00 * uu0 * pi * ( sum * a0 - vi1 )
        else
           vk2 = 0.5D+00 * uu0 * pi * ( sum * a0 - vi2 )
        end if

     end do

  else

     c0 = exp ( - x ) * sqrt ( 0.5D+00 * pi / x )

     do l = 1, 2
        vv = vv0 * l * l
        sum = 1.0D+00
        r = 1.0D+00
        do k = 1, k0
           r = 0.125D+00 * r * ( vv - ( 2.0D+00 * k - 1.0D+00 ) ** 2 ) / ( k * x )
           sum = sum + r
        end do
        if ( l == 1 ) then
           vk1 = c0 * sum
        else
           vk2 = c0 * sum
        end if
     end do

  end if

  return
end subroutine ajyik