cpsi Subroutine

subroutine cpsi ( x, y, psr, psi )

  !*****************************************************************************80
  !
  !! CPSI computes the psi function for a 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:
  !
  !    16 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, Y, the real and imaginary parts 
  !    of the argument.
  !
  !    Output, real ( kind = 8 ) PSR, PSI, the real and imaginary parts
  !    of the function value.
  !
  implicit none

  real ( kind = 8 ), save, dimension ( 8 ) :: a = (/ &
       -0.8333333333333D-01, 0.83333333333333333D-02, &
       -0.39682539682539683D-02, 0.41666666666666667D-02, &
       -0.75757575757575758D-02, 0.21092796092796093D-01, &
       -0.83333333333333333D-01, 0.4432598039215686D+00 /)
  real ( kind = 8 ) ct2
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  real ( kind = 8 ) pi
  real ( kind = 8 ) psi
  real ( kind = 8 ) psr
  real ( kind = 8 ) ri
  real ( kind = 8 ) rr
  real ( kind = 8 ) th
  real ( kind = 8 ) tm
  real ( kind = 8 ) tn
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1
  real ( kind = 8 ) y
  real ( kind = 8 ) y1
  real ( kind = 8 ) z0
  real ( kind = 8 ) z2

  pi = 3.141592653589793D+00

  if ( y == 0.0D+00 .and. x == int ( x ) .and. x <= 0.0D+00 ) then

     psr = 1.0D+300
     psi = 0.0D+00

  else

     if ( x < 0.0D+00 ) then
        x1 = x
        y1 = y
        x = -x
        y = -y
     end if

     x0 = x

     if ( x < 8.0D+00 ) then
        n = 8 - int ( x )
        x0 = x + n
     end if

     if ( x0 == 0.0D+00 ) then
        if ( y /= 0.0D+00 ) then
           th = 0.5D+00 * pi
        else
           th = 0.0D+00
        end if
     else
        th = atan ( y / x0 )
     end if

     z2 = x0 * x0 + y * y
     z0 = sqrt ( z2 )
     psr = log ( z0 ) - 0.5D+00 * x0 / z2
     psi = th + 0.5D+00 * y / z2
     do k = 1, 8
        psr = psr + a(k) * z2 ** ( - k ) * cos ( 2.0D+00 * k * th )
        psi = psi - a(k) * z2 ** ( - k ) * sin ( 2.0D+00 * k * th )
     end do

     if ( x < 8.0D+00 ) then
        rr = 0.0D+00
        ri = 0.0D+00
        do k = 1, n
           rr = rr + ( x0 - k ) / ( ( x0 - k ) ** 2.0D+00 + y * y )
           ri = ri + y / ( ( x0 - k ) ** 2.0D+00 + y * y )
        end do
        psr = psr - rr
        psi = psi + ri
     end if

     if ( x1 < 0.0D+00 ) then
        tn = tan ( pi * x )
        tm = tanh ( pi * y )
        ct2 = tn * tn + tm * tm
        psr = psr + x / ( x * x + y * y ) + pi * ( tn - tn * tm * tm ) / ct2
        psi = psi - y / ( x * x + y * y ) - pi * tm * ( 1.0D+00 + tn * tn ) / ct2
        x = x1
        y = y1
     end if

  end if

  return
end subroutine cpsi