ch12n Subroutine

subroutine ch12n ( n, z, nm, chf1, chd1, chf2, chd2 )

  !*****************************************************************************80
  !
  !! CH12N computes Hankel functions of first and second kinds, complex argument.
  !
  !  Discussion:
  !
  !    Both the Hankel functions and their derivatives are computed.
  !
  !  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:
  !
  !    26 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, integer ( kind = 4 ) N, the order of the functions.
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, integer ( kind = 4 ) NM, the highest order computed.
  !
  !    Output, complex ( kind = 8 ) CHF1(0:n), CHD1(0:n), CHF2(0:n), CHD2(0:n), 
  !    the values of Hn(1)(z), Hn(1)'(z), Hn(2)(z), Hn(2)'(z).
  !
  implicit none

  integer ( kind = 4 ) n

  complex ( kind = 8 ) cbi(0:250)
  complex ( kind = 8 ) cbj(0:250)
  complex ( kind = 8 ) cbk(0:250)
  complex ( kind = 8 ) cby(0:250)
  complex ( kind = 8 ) cdi(0:250)
  complex ( kind = 8 ) cdj(0:250)
  complex ( kind = 8 ) cdk(0:250)
  complex ( kind = 8 ) cdy(0:250)
  complex ( kind = 8 ) chd1(0:n)
  complex ( kind = 8 ) chd2(0:n)
  complex ( kind = 8 ) cf1
  complex ( kind = 8 ) cfac
  complex ( kind = 8 ) chf1(0:n)
  complex ( kind = 8 ) chf2(0:n)
  complex ( kind = 8 ) ci
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nm
  real ( kind = 8 ) pi
  complex ( kind = 8 ) z
  complex ( kind = 8 ) zi

  ci = cmplx ( 0.0D+00, 1.0D+00, kind = 8 )
  pi = 3.141592653589793D+00

  if ( imag ( z ) < 0.0D+00 ) then

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

     do k = 0, nm
        chf1(k) = cbj(k) + ci * cby(k)
        chd1(k) = cdj(k) + ci * cdy(k)
     end do

     zi = ci * z
     call ciknb ( n, zi, nm, cbi, cdi, cbk, cdk )
     cfac = -2.0D+00 / ( pi * ci )

     do k = 0, nm
        chf2(k) = cfac * cbk(k)
        chd2(k) = cfac * ci * cdk(k)
        cfac = cfac * ci
     end do

  else if ( 0.0D+00 < imag ( z ) ) then

     zi = - ci * z
     call ciknb ( n, zi, nm, cbi, cdi, cbk, cdk )
     cf1 = -ci
     cfac = 2.0D+00 / ( pi * ci )

     do k = 0, nm
        chf1(k) = cfac * cbk(k)
        chd1(k) = -cfac * ci * cdk(k)
        cfac = cfac * cf1
     end do

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

     do k = 0, nm
        chf2(k) = cbj(k) - ci * cby(k)
        chd2(k) = cdj(k) - ci * cdy(k)
     end do

  else

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

     do k = 0, nm
        chf1(k) = cbj(k) + ci * cby(k)
        chd1(k) = cdj(k) + ci * cdy(k)
        chf2(k) = cbj(k) - ci * cby(k)
        chd2(k) = cdj(k) - ci * cdy(k)
     end do

  end if

  return
end subroutine ch12n