cchg Subroutine

subroutine cchg ( a, b, z, chg )

  !*****************************************************************************80
  !
  !! CCHG computes the confluent hypergeometric function.
  !
  !  Discussion:
  !
  !    This function computes the confluent hypergeometric function
  !    M(a,b,z) with real parameters a, b and complex argument z.
  !
  !  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, real ( kind = 8 ) A, B, parameter values.
  !
  !    Input, complex ( kind = 8 ) Z, the argument.
  !
  !    Output, complex ( kind = 8 ) CHG, the value of M(a,b,z).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) a1
  real ( kind = 8 ) b
  real ( kind = 8 ) ba
  complex ( kind = 8 ) cfac
  complex ( kind = 8 ) chg
  complex ( kind = 8 ) chg1
  complex ( kind = 8 ) chg2
  complex ( kind = 8 ) chw
  complex ( kind = 8 ) ci
  complex ( kind = 8 ) cr
  complex ( kind = 8 ) cr1
  complex ( kind = 8 ) cr2
  complex ( kind = 8 ) crg
  complex ( kind = 8 ) cs1
  complex ( kind = 8 ) cs2
  complex ( kind = 8 ) cy0
  complex ( kind = 8 ) cy1
  real ( kind = 8 ) g1
  real ( kind = 8 ) g2
  real ( kind = 8 ) g3
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) la
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nl
  integer ( kind = 4 ) ns
  real ( kind = 8 ) phi
  real ( kind = 8 ) pi
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) y
  complex ( kind = 8 ) z
  complex ( kind = 8 ) z0

  pi = 3.141592653589793D+00
  ci = cmplx ( 0.0D+00, 1.0D+00 )
  a0 = a
  a1 = a
  z0 = z

  if ( b == 0.0D+00 .or. b == - int ( abs ( b ) ) ) then
     chg = cmplx ( 1.0D+30, 0.0D+00, kind = 8 )
  else if ( a == 0.0D+00 .or. z == 0.0D+00 ) then
     chg = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
  else if ( a == -1.0D+00 ) then
     chg = 1.0D+00 - z / b
  else if ( a == b ) then
     chg = exp ( z )
  else if ( a - b == 1.0D+00 ) then
     chg = ( 1.0D+00 + z / b ) * exp ( z )
  else if ( a == 1.0D+00 .and. b == 2.0D+00 ) then
     chg = ( exp ( z ) - 1.0D+00 ) / z
  else if ( a == int ( a ) .and. a < 0.0D+00 ) then
     m = int ( - a )
     cr = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     chg = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
     do k = 1, m
        cr = cr * ( a + k - 1.0D+00 ) / k / ( b + k - 1.0D+00 ) * z
        chg = chg + cr
     end do
  else

     x0 = real ( z, kind = 8 )
     if ( x0 < 0.0D+00 ) then
        a = b - a
        a0 = a
        z = - z
     end if

     if ( a < 2.0D+00 ) then
        nl = 0
     else
        nl = 1
        la = int ( a )
        a = a - la - 1.0D+00
     end if

     do n = 0, nl

        if ( 2.0D+00 <= a0 ) then
           a = a + 1.0D+00
        end if

        if ( cdabs ( z ) < 20.0D+00 + abs ( b ) .or. a < 0.0D+00 ) then

           chg = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           crg = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           do j = 1, 500
              crg = crg * ( a + j - 1.0D+00 ) / ( j * ( b + j - 1.0D+00 ) ) * z
              chg = chg + crg
              if ( abs ( ( chg - chw ) / chg ) < 1.0D-15 ) then
                 exit
              end if
              chw = chg
           end do

        else

           call gammaf ( a, g1 )
           call gammaf ( b, g2 )
           ba = b - a
           call gammaf ( ba, g3 )
           cs1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           cs2 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           cr1 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )
           cr2 = cmplx ( 1.0D+00, 0.0D+00, kind = 8 )

           do i = 1, 8
              cr1 = - cr1 * (     a + i - 1.0D+00 ) * ( a - b + i ) / ( z * i )
              cr2 =   cr2 * ( b - a + i - 1.0D+00 ) * ( i - a ) / ( z * i )
              cs1 = cs1 + cr1
              cs2 = cs2 + cr2
           end do

           x = real ( z, kind = 8 )
           y = imag ( z )

           if ( x == 0.0D+00 .and. 0.0D+00 <= y ) then
              phi = 0.5D+00 * pi
           else if ( x == 0.0D+00 .and. y <= 0.0D+00 ) then
              phi = -0.5D+00 * pi
           else
              phi = atan ( y / x )
           end if

           if ( -1.5D+00 * pi < phi .and. phi <= -0.5 * pi ) then
              ns = -1
           else if ( -0.5D+00 * pi < phi .and. phi < 1.5D+00 * pi ) then
              ns = 1
           end if

           if ( y == 0.0D+00 ) then
              cfac = cos ( pi * a )
           else
              cfac = exp ( ns * ci * pi * a )
           end if

           chg1 = g2 / g3 * z ** ( - a ) * cfac * cs1
           chg2 = g2 / g1 * exp ( z ) * z ** ( a - b ) * cs2
           chg = chg1 + chg2

        end if

        if ( n == 0 ) then
           cy0 = chg
        else if ( n == 1 ) then
           cy1 = chg
        end if

     end do

     if ( 2.0D+00 <= a0 ) then
        do i = 1, la - 1
           chg = ( ( 2.0D+00 * a - b + z ) * cy1 + ( b - a ) * cy0 ) / a
           cy0 = cy1
           cy1 = chg
           a = a + 1.0D+00
        end do
     end if

     if ( x0 < 0.0D+00 ) then
        chg = chg * exp ( - z )
     end if

  end if

  a = a1
  z = z0

  return
end subroutine cchg