!************************************************************** !* 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