!*********************************************************************
!*          Calculate the Psi function for a complex argument        *
!* ----------------------------------------------------------------- *
!* EXPLANATION:                                                      *
!*                                                                   *
!*      Purpose: This program computes the psi function psi(z)       *
!*               for a complex argument using subroutine CPSI        *
!*      Input :  x   --- Real part of z                              *
!*               y   --- Imaginary part of z                         *
!*      Output:  PSR --- Real part of psi(z)                         *
!*               PSI --- Imaginary part of psi(z)                    *
!*      Examples:                                                    *
!*                  x       y      Re[psi(z)]     Im[psi(z)]         *
!*                -------------------------------------------        *
!*                 3.0     2.0     1.16459152      .67080728         *
!*                 3.0    -2.0     1.16459152     -.67080728         * 
!*                -3.0     2.0     1.39536075     2.62465344         *
!*                -3.0    -2.0     1.39536075    -2.62465344         *
!*                                                                   *
!* NOTE:          Psi(z) = Gamma'(z) / Gamma(z)                      *
!* ----------------------------------------------------------------- *
!* SAMPLE RUN:                                                       *
!*                                                                   *
!*  Please enter x and y (z=x+iy): 3 2                               *
!*                                                                   *
!*     x       y      Re[Psi(z)]      Im[Psi(z)]                     *
!*   ----------------------------------------------                  *
!*    3.0     2.0    0.11645915E+01  0.67080728E+00                  *
!*                                                                   *
!* ----------------------------------------------------------------- *
!* REFERENCE:"Fortran Routines for Computation of Special Functions, *
!*            jin.ece.uiuc.edu/routines/routines.html".              *
!*                                                                   *
!*                                 F90 Release By J-P Moreau, Paris. *
!*                                        (www.jpmoreau.fr)          *
!*********************************************************************  
      
PROGRAM MCPSI
        
DOUBLE PRECISION X,Y,PSR,PSI
   
   WRITE(*,*)
        
   WRITE(*,10,advance='no'); READ *, X, Y
   WRITE(*,*)
        
   WRITE(*,*)'   x       y      Re[Psi(z)]      Im[Psi(z)]'
        
   WRITE(*,*)' ----------------------------------------------'
        
   CALL CPSI(X,Y,PSR,PSI)
        
   WRITE(*,20)X,Y,PSR,PSI
   WRITE(*,*)
   STOP
10 FORMAT(' Please enter x and y (z=x+iy): ')

20 FORMAT(1X,F5.1,3X,F5.1,2X,2E16.8)

   END

   SUBROUTINE CPSI(X,Y,PSR,PSI)
!  =============================================
!    Purpose: Compute the psi function for a
!             complex argument
!    Input :  x   --- Real part of z
!             y   --- Imaginary part of z
!    Output:  PSR --- Real part of psi(z)
!             PSI --- Imaginary part of psi(z)
!  =============================================

   IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   
   DIMENSION A(8)
   
   DATA A/-.8333333333333D-01,.83333333333333333D-02,   &
          
          -.39682539682539683D-02,.41666666666666667D-02, &
          
          -.75757575757575758D-02,.21092796092796093D-01, &
          
          -.83333333333333333D-01,.4432598039215686D0/
   
   
   PI=3.141592653589793D0

   IF (Y.EQ.0.0D0.AND.X.EQ.INT(X).AND.X.LE.0.0D0) THEN
     
     PSR=1.0D+300
     
     PSI=0.0D0
   
   ELSE 
     IF (X.LT.0.0D0) THEN
       
       X1=X
       
       Y1=Y
       
       X=-X
       
       Y=-Y
     
     ENDIF


     X0=X


     IF (X.LT.8.0D0) THEN
       
       N=8-INT(X)
       
       X0=X+N
     
     ENDIF
     
     IF (X0.EQ.0.0D0.AND.Y.NE.0.0D0) TH=0.5D0*PI
     
     IF (X0.NE.0.0D0) TH=DATAN(Y/X0)
   
     Z2=X0*X0+Y*Y
     
     Z0=DSQRT(Z2)
     
     PSR=DLOG(Z0)-0.5D0*X0/Z2
     
     PSI=TH+0.5D0*Y/Z2

     DO 10 K=1,8
       
       PSR=PSR+A(K)*Z2**(-K)*DCOS(2.0D0*K*TH)

10     PSI=PSI-A(K)*Z2**(-K)*DSIN(2.0D0*K*TH)
           
     IF (X.LT.8.0D0) THEN
       
       RR=0.0D0
       
       RI=0.0D0
       
       DO 20 K=1,N
         
         RR=RR+(X0-K)/((X0-K)**2.0D0+Y*Y)

20       RI=RI+Y/((X0-K)**2.0D0+Y*Y)
       
       PSR=PSR-RR
       
       PSI=PSI+RI
     
     ENDIF

     IF (X1.LT.0.0D0) THEN
       
       TN=DTAN(PI*X)
       
       TM=DTANH(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.0D0+TN*TN)/CT2
       
       X=X1
       
       Y=Y1
     
     ENDIF
   
   ENDIF
   
   RETURN
   
   END


!end of file mcpsi.f90