cjyna Subroutine

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

  !*****************************************************************************80
  !
  !! CJYNA: Bessel functions and 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:
  !
  !    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, 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 ) a0
  complex ( kind = 8 ) cbj(0:n)
  complex ( kind = 8 ) cbj0
  complex ( kind = 8 ) cbj1
  complex ( kind = 8 ) cby(0:n)
  complex ( kind = 8 ) cby0
  complex ( kind = 8 ) cby1
  complex ( kind = 8 ) cdj(0:n)
  complex ( kind = 8 ) cdj0
  complex ( kind = 8 ) cdj1
  complex ( kind = 8 ) cdy(0:n)
  complex ( kind = 8 ) cdy0
  complex ( kind = 8 ) cdy1
  complex ( kind = 8 ) cf
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cf2
  complex ( kind = 8 ) cg0
  complex ( kind = 8 ) cg1
  complex ( kind = 8 ) ch0
  complex ( kind = 8 ) ch1
  complex ( kind = 8 ) ch2
  complex ( kind = 8 ) cj0
  complex ( kind = 8 ) cj1
  complex ( kind = 8 ) cjk
  complex ( kind = 8 ) cp11
  complex ( kind = 8 ) cp12
  complex ( kind = 8 ) cp21
  complex ( kind = 8 ) cp22
  complex ( kind = 8 ) cs
  complex ( kind = 8 ) cyk
  complex ( kind = 8 ) cyl1
  complex ( kind = 8 ) cyl2
  complex ( kind = 8 ) cylk
  integer ( kind = 4 ) k
  integer ( kind = 4 ) lb
  integer ( kind = 4 ) lb0
  integer ( kind = 4 ) m
    ! integer ( kind = 4 ) msta1
  ! integer ( kind = 4 ) msta2
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  real ( kind = 8 ) wa
  real ( kind = 8 ) ya0
  real ( kind = 8 ) ya1
  real ( kind = 8 ) yak
  complex ( kind = 8 ) z

  pi = 3.141592653589793D+00
  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

  call cjy01 ( z, cbj0, cdj0, cbj1, cdj1, cby0, cdy0, cby1, cdy1 )
  cbj(0) = cbj0
  cbj(1) = cbj1
  cby(0) = cby0
  cby(1) = cby1
  cdj(0) = cdj0
  cdj(1) = cdj1
  cdy(0) = cdy0
  cdy(1) = cdy1

  if ( n <= 1 ) then
     return
  end if

  if ( n < int ( 0.25D+00 * a0 ) ) then

     cj0 = cbj0
     cj1 = cbj1
     do k = 2, n
        cjk = 2.0D+00 * ( k - 1.0D+00 ) / z * cj1 - cj0
        cbj(k) = cjk
        cj0 = cj1
        cj1 = cjk
     end do

  else

     m = msta1 ( a0, 200 )

     if ( m < n ) then
        nm = m
     else
        m = msta2 ( a0, n, 15 )
     end if

     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
        cf2 = cf1
        cf1 = cf
     end do

     if ( abs ( cbj1 ) < abs ( cbj0 ) ) then
        cs = cbj0 / cf
     else
        cs = cbj1 / cf2
     end if

     do k = 0, nm
        cbj(k) = cs * cbj(k)
     end do

  end if

  do k = 2, nm
     cdj(k) = cbj(k-1) - k / z * cbj(k)
  end do
  ya0 = abs ( cby0 )
  lb = 0
  cg0 = cby0
  cg1 = cby1
  do k = 2, nm
     cyk = 2.0D+00 * ( k - 1.0D+00 ) / z * cg1 - cg0
     if ( abs ( cyk ) <= 1.0D+290 ) then         
        yak = abs ( cyk )
        ya1 = abs ( cg0 )
        if ( yak < ya0 .and. yak < ya1 ) then
           lb = k
        end if
        cby(k) = cyk
        cg0 = cg1
        cg1 = cyk
     end if
  end do

  if ( 4 < lb  .and. imag ( z ) /= 0.0D+00 ) then

     do

        if ( lb == lb0 ) then
           exit
        end if

        ch2 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        ch1 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        lb0 = lb
        do k = lb, 1, -1
           ch0 = 2.0D+00 * k / z * ch1 - ch2
           ch2 = ch1
           ch1 = ch0
        end do
        cp12 = ch0
        cp22 = ch2
        ch2 = cmplx ( 0.0D+00, 0.0D+00, kind = 8 )
        ch1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
        do k = lb, 1, -1
           ch0 = 2.0D+00 * k / z * ch1 - ch2
           ch2 = ch1
           ch1 = ch0
        end do
        cp11 = ch0
        cp21 = ch2

        if ( lb == nm ) then
           cbj(lb+1) = 2.0D+00 * lb / z * cbj(lb) - cbj(lb-1)
        end if

        if ( abs ( cbj(1) ) < abs ( cbj(0) ) ) then
           cby(lb+1) = ( cbj(lb+1) * cby0 - 2.0D+00 * cp11 / ( pi * z ) ) / cbj(0)
           cby(lb) = ( cbj(lb) * cby0 + 2.0D+00 * cp12 / ( pi * z ) ) / cbj(0)
        else
           cby(lb+1) = ( cbj(lb+1) * cby1 - 2.0D+00 * cp21 / ( pi * z ) ) / cbj(1)
           cby(lb) = ( cbj(lb) * cby1 + 2.0D+00 * cp22 / ( pi * z ) )  / cbj(1)
        end if

        cyl2 = cby(lb+1)
        cyl1 = cby(lb)
        do k = lb - 1, 0, -1
           cylk = 2.0D+00 * ( k + 1.0D+00 ) / z * cyl1 - cyl2
           cby(k) = cylk
           cyl2 = cyl1
           cyl1 = cylk
        end do

        cyl1 = cby(lb)
        cyl2 = cby(lb+1)
        do k = lb + 1, nm - 1
           cylk = 2.0D+00 * k / z * cyl2 - cyl1
           cby(k+1) = cylk
           cyl1 = cyl2
           cyl2 = cylk
        end do

        do k = 2, nm
           wa = abs ( cby(k) )
           if ( wa < abs ( cby(k-1) ) ) then
              lb = k
           end if
        end do

     end do

  end if

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

  return
end subroutine cjyna