cjynb Subroutine

subroutine cjynb ( n, z, nm, cbj, cdj, cby, cdy )

  !*****************************************************************************80
  !
  !! CJYNB: Bessel functions, derivatives, Jn(z) and Yn(z) of complex argument.
  !
  !  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, integer ( kind = 4 ) N, the order of Jn(z) and Yn(z).
  !
  !    Input, complex ( kind = 8 ) Z, the argument of Jn(z) and Yn(z).
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, complex ( kind = 8 ) CBJ(0:N), CDJ(0:N), CBY(0:N), CDY(0:N), 
  !    the values of Jn(z), Jn'(z), Yn(z), Yn'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ), save, dimension ( 4 ) :: a = (/ &
       -0.7031250000000000D-01, 0.1121520996093750D+00, &
       -0.5725014209747314D+00, 0.6074042001273483D+01 /)
  real ( kind = 8 ) a0
  real ( kind = 8 ), save, dimension ( 4 ) :: a1 = (/ &
       0.1171875000000000D+00,-0.1441955566406250D+00, &
       0.6765925884246826D+00,-0.6883914268109947D+01 /)
  real ( kind = 8 ), save, dimension ( 4 ) :: b = (/  &
       0.7324218750000000D-01,-0.2271080017089844D+00, &
       0.1727727502584457D+01,-0.2438052969955606D+02 /)
  real ( kind = 8 ), save, dimension ( 4 ) :: b1 = (/ &
       -0.1025390625000000D+00,0.2775764465332031D+00, &
       -0.1993531733751297D+01,0.2724882731126854D+02 /)
  complex ( kind = 8 ) cbj(0:n)
  complex ( kind = 8 ) cbj0
  complex ( kind = 8 ) cbj1
  complex ( kind = 8 ) cbjk
  complex ( kind = 8 ) cbs
  complex ( kind = 8 ) cby(0:n)
  complex ( kind = 8 ) cby0
  complex ( kind = 8 ) cby1
  complex ( kind = 8 ) cdj(0:n)
  complex ( kind = 8 ) cdy(0:n)
  complex ( kind = 8 ) ce
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) cp0
  complex ( kind = 8 ) cp1
  complex ( kind = 8 ) cq0
  complex ( kind = 8 ) cq1
  complex ( kind = 8 ) cs0
  complex ( kind = 8 ) csu
  complex ( kind = 8 ) csv
  complex ( kind = 8 ) ct1
  complex ( kind = 8 ) ct2
  complex ( kind = 8 ) cu
  complex ( kind = 8 ) cyy
  real ( kind = 8 ) el
  integer ( kind = 4 ) k
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  real ( kind = 8 ) r2p
  real ( kind = 8 ) y0
  complex ( kind = 8 ) z

  el = 0.5772156649015329D+00
  pi = 3.141592653589793D+00
  r2p = 0.63661977236758D+00
  y0 = abs ( imag ( z ) )
  a0 = abs ( z )
  nm = n

  if ( a0 < 1.0D-100 ) then
     do k = 0, n
        cbj(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cdj(k) = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        cby(k) = - cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
        cdy(k) = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
     end do
     cbj(0) = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     cdj(1) = cmplx ( 0.5D+00, 0.0D+00, kind = 8 )
     return
  end if

  if ( a0 <= 300.0D+00 .or. 80 < n ) then

     if ( n == 0 ) then
        nm = 1
     end if
     m = msta1 ( a0, 200 )
     if ( m < nm ) then
        nm = m
     else
        m = msta2 ( a0, nm, 15 )
     end if

     cbs = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     csu = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     csv = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     cf2 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
     cf1 = cmplx ( 1.0D-30, 0.0D+00, kind = 8 )

     do k = m, 0, -1
        cf = 2.0D+00 * ( k + 1.0D+00 ) / z * cf1 - cf2
        if ( k <= nm ) then
           cbj(k) = cf
        end if
        if ( k == 2 * int ( k / 2 ) .and. k .ne. 0 ) then
           if ( y0 <= 1.0D+00 ) then
              cbs = cbs + 2.0D+00 * cf
           else
              cbs = cbs + ( -1.0D+00 ) ** ( k / 2 ) * 2.0D+00 * cf
           end if
           csu = csu + ( -1.0D+00 ) ** ( k / 2 ) * cf / k
        else if ( 1 < k ) then
           csv = csv + ( -1.0D+00 ) ** ( k / 2 ) * k / ( k * k - 1.0D+00 ) * cf
        end if
        cf2 = cf1
        cf1 = cf
     end do

     if ( y0 <= 1.0D+00 ) then
        cs0 = cbs + cf
     else
        cs0 = ( cbs + cf ) / cos ( z )
     end if

     do k = 0, nm
        cbj(k) = cbj(k) / cs0
     end do

     ce = log ( z / 2.0D+00 ) + el
     cby(0) = r2p * ( ce * cbj(0) - 4.0D+00 * csu / cs0 )
     cby(1) = r2p * ( - cbj(0) / z + ( ce - 1.0D+00 ) * cbj(1) &
          - 4.0D+00 * csv / cs0 )

  else

     ct1 = z - 0.25D+00 * pi
     cp0 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, 4
        cp0 = cp0 + a(k) * z ** ( - 2 * k )
     end do
     cq0 = -0.125D+00 / z
     do k = 1, 4
        cq0 = cq0 + b(k) * z ** ( - 2 * k - 1 )
     end do
     cu = sqrt ( r2p / z )
     cbj0 = cu * ( cp0 * cos ( ct1 ) - cq0 * sin ( ct1 ) )
     cby0 = cu * ( cp0 * sin ( ct1 ) + cq0 * cos ( ct1 ) )
     cbj(0) = cbj0
     cby(0) = cby0
     ct2 = z - 0.75D+00 * pi
     cp1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, 4
        cp1 = cp1 + a1(k) * z ** ( - 2 * k )
     end do
     cq1 = 0.375D+00 / z
     do k = 1, 4
        cq1 = cq1 + b1(k) * z ** ( - 2 * k - 1 )
     end do
     cbj1 = cu * ( cp1 * cos ( ct2 ) - cq1 * sin ( ct2 ) )
     cby1 = cu * ( cp1 * sin ( ct2 ) + cq1 * cos ( ct2 ) )
     cbj(1) = cbj1
     cby(1) = cby1
     do k = 2, nm
        cbjk = 2.0D+00 * ( k - 1.0D+00 ) / z * cbj1 - cbj0
        cbj(k) = cbjk
        cbj0 = cbj1
        cbj1 = cbjk
     end do
  end if

  cdj(0) = -cbj(1)
  do k = 1, nm
     cdj(k) = cbj(k-1) - k / z * cbj(k)
  end do

  if ( 1.0D+00 < abs ( cbj(0) ) ) then
     cby(1) = ( cbj(1) * cby(0) - 2.0D+00 / ( pi * z ) ) / cbj(0)
  end if

  do k = 2, nm
     if ( abs ( cbj(k-2) ) <= abs ( cbj(k-1) ) ) then
        cyy = ( cbj(k) * cby(k-1) - 2.0D+00 / ( pi * z ) ) / cbj(k-1)
     else
        cyy = ( cbj(k) * cby(k-2) - 4.0D+00 * ( k - 1.0D+00 ) &
             / ( pi * z * z ) ) / cbj(k-2)
     end if
     cby(k) = cyy
  end do

  cdy(0) = -cby(1)
  do k = 1, nm
     cdy(k) = cby(k-1) - k / z * cby(k)
  end do

  return
end subroutine cjynb