!****************************************************************** !* Purpose: This program computes the confluent * !* hypergeometric function M(a,b,x) using * !* subroutine CHGM * !* Input : a --- Parameter * !* b --- Parameter ( b <> 0,-1,-2,... ) * !* x --- Argument * !* Output: HG --- M(a,b,x) * !* Example: * !* a b x M(a,b,x) * !* ----------------------------------------- * !* 1.5 2.0 20.0 .1208527185D+09 * !* 4.5 2.0 20.0 .1103561117D+12 * !* -1.5 2.0 20.0 .1004836854D+05 * !* -4.5 2.0 20.0 -.3936045244D+03 * !* 1.5 2.0 50.0 .8231906643D+21 * !* 4.5 2.0 50.0 .9310512715D+25 * !* -1.5 2.0 50.0 .2998660728D+16 * !* -4.5 2.0 50.0 -.1807604439D+13 * !* -------------------------------------------------------------- * !* 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 MCHGM IMPLICIT DOUBLE PRECISION (A-H,O-Z) WRITE(*,*)' ' WRITE(*,*)' a b x M(a,b,x)' WRITE(*,*)' -----------------------------------------' A=1.5D0; B=2.D0; X=20.D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=4.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=-1.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=-4.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=1.5D0; X=50.D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=4.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=-1.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG A=-4.5D0 CALL CHGM(A,B,X,HG) WRITE(*,10)A,B,X,HG 10 FORMAT(1X,F5.1,3X,F5.1,3X,F5.1,D20.10) END SUBROUTINE CHGM(A,B,X,HG) ! =================================================== ! Purpose: Compute confluent hypergeometric function ! M(a,b,x) ! Input : a --- Parameter ! b --- Parameter ( b <> 0,-1,-2,... ) ! x --- Argument ! Output: HG --- M(a,b,x) ! Routine called: GAMMA for computing â(x) ! =================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) PI=3.141592653589793D0 A0=A A1=A X0=X HG=0.0D0 IF (B.EQ.0.0D0.OR.B.EQ.-ABS(INT(B))) THEN HG=1.0D+300 ELSE IF (A.EQ.0.0D0.OR.X.EQ.0.0D0) THEN HG=1.0D0 ELSE IF (A.EQ.-1.0D0) THEN HG=1.0D0-X/B ELSE IF (A.EQ.B) THEN HG=DEXP(X) ELSE IF (A-B.EQ.1.0D0) THEN HG=(1.0D0+X/B)*DEXP(X) ELSE IF (A.EQ.1.0D0.AND.B.EQ.2.0D0) THEN HG=(DEXP(X)-1.0D0)/X ELSE IF (A.EQ.INT(A).AND.A.LT.0.0D0) THEN M=INT(-A) R=1.0D0 HG=1.0D0 DO 10 K=1,M R=R*(A+K-1.0D0)/K/(B+K-1.0D0)*X 10 HG=HG+R ENDIF IF (HG.NE.0.0D0) RETURN IF (X.LT.0.0D0) THEN A=B-A A0=A X=DABS(X) ENDIF IF (A.LT.2.0D0) NL=0 IF (A.GE.2.0D0) THEN NL=1 LA=INT(A) A=A-LA-1.0D0 ENDIF DO 30 N=0,NL IF (A0.GE.2.0D0) A=A+1.0D0 IF (X.LE.30.0D0+DABS(B).OR.A.LT.0.0D0) THEN HG=1.0D0 RG=1.0D0 DO 15 J=1,500 RG=RG*(A+J-1.0D0)/(J*(B+J-1.0D0))*X HG=HG+RG IF (DABS(RG/HG).LT.1.0D-15) GO TO 25 15 CONTINUE ELSE CALL GAMMA(A,TA) CALL GAMMA(B,TB) XG=B-A CALL GAMMA(XG,TBA) SUM1=1.0D0 SUM2=1.0D0 R1=1.0D0 R2=1.0D0 DO 20 I=1,8 R1=-R1*(A+I-1.0D0)*(A-B+I)/(X*I) R2=-R2*(B-A+I-1.0D0)*(A-I)/(X*I) SUM1=SUM1+R1 20 SUM2=SUM2+R2 HG1=TB/TBA*X**(-A)*DCOS(PI*A)*SUM1 HG2=TB/TA*DEXP(X)*X**(A-B)*SUM2 HG=HG1+HG2 ENDIF 25 IF (N.EQ.0) Y0=HG IF (N.EQ.1) Y1=HG 30 CONTINUE IF (A0.GE.2.0D0) THEN DO 35 I=1,LA-1 HG=((2.0D0*A-B+X)*Y1+(B-A)*Y0)/A Y0=Y1 Y1=HG 35 A=A+1.0D0 ENDIF IF (X0.LT.0.0D0) HG=HG*DEXP(X0) A=A1 X=X0 RETURN END SUBROUTINE GAMMA(X,GA) ! ================================================== ! Purpose: Compute gamma function â(x) ! Input : x --- Argument of â(x) ! ( x is not equal to 0,-1,-2,úúú) ! Output: GA --- â(x) ! ================================================== IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION G(26) PI=3.141592653589793D0 IF (X.EQ.INT(X)) THEN IF (X.GT.0.0D0) THEN GA=1.0D0 M1=X-1 DO 10 K=2,M1 10 GA=GA*K ELSE GA=1.0D+300 ENDIF ELSE IF (DABS(X).GT.1.0D0) THEN Z=DABS(X) M=INT(Z) R=1.0D0 DO 15 K=1,M 15 R=R*(Z-K) Z=Z-M ELSE Z=X ENDIF DATA G/1.0D0,0.5772156649015329D0, & -0.6558780715202538D0, -0.420026350340952D-1, & 0.1665386113822915D0,-.421977345555443D-1, & -.96219715278770D-2, .72189432466630D-2, & -.11651675918591D-2, -.2152416741149D-3, & .1280502823882D-3, -.201348547807D-4, & -.12504934821D-5, .11330272320D-5, & -.2056338417D-6, .61160950D-8, & .50020075D-8, -.11812746D-8, & .1043427D-9, .77823D-11, & -.36968D-11, .51D-12, & -.206D-13, -.54D-14, .14D-14, .1D-15/ GR=G(26) DO 20 K=25,1,-1 20 GR=GR*Z+G(K) GA=1.0D0/(GR*Z) IF (DABS(X).GT.1.0D0) THEN GA=GA*R IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X)) ENDIF ENDIF RETURN END !end of file mchgm.f90