!**************************************************************      

!*       Calculate Function Gamma with a complex argument     * 

!* ---------------------------------------------------------- *
!* EXPLANATION:                                               *

!* Purpose: This program computes the gamma function G(z)     * 

!*          or Ln[G(z)] for a complex argument using          *

!*          subroutine CGAMA                                  *

!* Input :  x  --- Real part of z                             *

!*          y  --- Imaginary part of z                        *

!*          KF --- Function code                              *

!*          KF=0 for Ln[G(z)]                                 *

!*          KF=1 for G(z)                                     *

!* Output:  GR --- Real part of Ln[G(z)] or G(z)              *

!*          GI --- Imaginary part of Ln[G(z)] or G(z)         *

!* Examples:                                                  *

!*    x         y           Re[G(z)]           Im[G(z)]       *

!*  --------------------------------------------------------  *

!*   2.50      5.00     .2267360319D-01    -.1172284404D-01   *
!*   5.00     10.00     .1327696517D-01     .3639011746D-02   *

!*   2.50     -5.00     .2267360319D-01     .1172284404D-01   *

!*   5.00    -10.00     .1327696517D-01    -.3639011746D-02   *

!*                                                            *             

!*    x         y          Re[LnG(z)]         Im[LnG(z)]      *

!*  --------------------------------------------------------  *

!*   2.50      5.00    -.3668103262D+01     .5806009801D+01   *

!*   5.00     10.00    -.4285507444D+01     .1911707090D+02   *

!*   2.50     -5.00    -.3668103262D+01    -.5806009801D+01   *

!*   5.00    -10.00    -.4285507444D+01    -.1911707090D+02   *

!* ---------------------------------------------------------- *
!* SAMPLE RUNS:                                               *
!*                                                            *
!* Please enter KF, x and y: 1 2.5 5.0                        *
!*                                                            *
!*    x         y           Re[G(z)]           Im[G(z)]       *
!*  --------------------------------------------------------  *
!*   2.50      5.00    0.2267360319D-01   -0.1172284404D-01   *
!*                                                            *
!* Please enter KF, zx and zy: 0 2.5 5.0                      *
!*                                                            *
!*    x         y          Re[LnG(z)]         Im[LnG(z)]      *
!*  --------------------------------------------------------  *
!*   2.50      5.00   -0.3668103262D+01    0.5806009801D+01   *
!*                                                            *
!* ---------------------------------------------------------- *
!* 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 MCGAMA
    
DOUBLE PRECISION X,Y,GR,GI
    PRINT *,' '
    
    WRITE(*,10,advance='no')
    
    READ *, KF, X, Y
    
    WRITE(*,*)
    
    IF (KF.EQ.1) THEN
        
      WRITE(*,*)'       x         y           Re[G(z)]',  &

      '           Im[G(z)]'
    
    ELSE
        
      WRITE(*,*)'       x         y          Re[LnG(z)]', &
                  
      '         Im[LnG(z)]'
    
    ENDIF
    
    WRITE(*,*)'    ------------------------------------',   &
              
    '---------------------'
    
    CALL CGAMA(X,Y,KF,GR,GI)
    
    WRITE(*,20)X,Y,GR,GI
    PRINT *,' '
    STOP

10  FORMAT('  Please enter KF, x and y: ')

20  FORMAT(1X,2F10.2,2D20.10)
    
    END


    

    SUBROUTINE CGAMA(X,Y,KF,GR,GI)

    ! ===========================================================
    
!       Purpose: Compute the gamma function G(z) or Ln[G(z)]
    
!                for a complex argument
    
!       Input :  x  --- Real part of z
    
!                y  --- Imaginary part of z
    
!                KF --- Function code
    
!                       KF=0 for Ln[G(z)]
    
!                       KF=1 for G(z)
    
!       Output:  GR --- Real part of Ln[G(z)] or G(z)

    !                GI --- Imaginary part of Ln[G(z)] or G(z)
    
! ===========================================================
    
    IMPLICIT DOUBLE PRECISION (A-H,O-Z)
    
    DIMENSION A(10)
    
    PI=3.141592653589793D0
    
    DATA A/8.333333333333333D-02,-2.777777777777778D-03,  &
           
           7.936507936507937D-04,-5.952380952380952D-04,  &
           
           8.417508417508418D-04,-1.917526917526918D-03,  &
           
           6.410256410256410D-03,-2.955065359477124D-02,  &
           
           1.796443723688307D-01,-1.39243221690590D+00/
    
    IF (Y.EQ.0.0D0.AND.X.EQ.INT(X).AND.X.LE.0.0D0) THEN
       
      GR=1.0D+300
       
      GI=0.0D0
       
      RETURN
    
    ELSE IF (X.LT.0.0D0) THEN
       
      X1=X
       
      Y1=Y
       
      X=-X
       
      Y=-Y
    
    ENDIF
    
    X0=X
    
    IF (X.LE.7.0) THEN
       
      NA=INT(7-X)
       
      X0=X+NA
    
    ENDIF
    
    Z1=DSQRT(X0*X0+Y*Y)
    
    TH=DATAN(Y/X0)
    
    GR=(X0-.5D0)*DLOG(Z1)-TH*Y-X0+0.5D0*DLOG(2.0D0*PI)
    
    GI=TH*(X0-0.5D0)+Y*DLOG(Z1)-Y
    
    DO K=1,10
       
      T=Z1**(1-2*K)
       
      GR=GR+A(K)*T*DCOS((2.0D0*K-1.0D0)*TH)
       
      GI=GI-A(K)*T*DSIN((2.0D0*K-1.0D0)*TH)
    END DO
    
    IF (X.LE.7.0) THEN
       
      GR1=0.0D0
       
      GI1=0.0D0
       
      DO J=0,NA-1
          
        GR1=GR1+.5D0*DLOG((X+J)**2+Y*Y)
          
        GI1=GI1+DATAN(Y/(X+J))
      END DO
       
      GR=GR-GR1
       
      GI=GI-GI1
    
    ENDIF
    
    IF (X1.LT.0.0D0) THEN
       
      Z1=DSQRT(X*X+Y*Y)
       
      TH1=DATAN(Y/X)
       
      SR=-DSIN(PI*X)*DCOSH(PI*Y)
       
      SI=-DCOS(PI*X)*DSINH(PI*Y)
       
      Z2=DSQRT(SR*SR+SI*SI)
       
      TH2=DATAN(SI/SR)
       
      IF (SR.LT.0.0D0) TH2=PI+TH2
       
      GR=DLOG(PI/(Z1*Z2))-GR
       
      GI=-TH1-TH2-GI
       
      X=X1
       
      Y=Y1
    
    ENDIF
    
    IF (KF.EQ.1) THEN
       
      G0=DEXP(GR)
       
      GR=G0*DCOS(GI)
       
      GI=G0*DSIN(GI)
    
    ENDIF
    
    RETURN
    
    END


!end of file mcgama.f90