C
C Subroutine file 'hg.for' to calculate some hypergeometric functions C C Version: 6.20 C Date: 2008, June 4 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm C C....................................................................... C C This file consists of the following external procedures: C HG2F1...Subroutine designed to calculate hypergeometric function C 2F1(A,B;C;X). C HG2F1 C HGF1... Subroutine designed to calculate hypergeometric function C F1(A,B1,1,C;X1,X2). C HGF1 C HGF2... Subroutine designed to calculate hypergeometric function C F2(A,B1,1,C1,C2;X1,X2). C HGF2 C HGFM2...Subroutine designed to calculate modified hypergeometric C function F~2(A,B1,1,C1,1;X1,X2). C HGFM2 C C======================================================================= C C C SUBROUTINE HG2F1(A,B,C,X,ERR,F) REAL A,B,C,X,ERR,F C C Subroutine designed to calculate hypergeometric function 2F1(A,B;C;X). C C Input: C A,B,C.. Parameters of hypergeometric function 2F1(A,B;C;X). C The parameters are assumed to be positive. C X... Variable of hypergeometric function 2F1(A,B;C;X). C The subroutine is designed for non-negative X sufficiently C smaller than 1. C ERR... Absolute error of the output value. C Output: C F... Value of hypergeometric function 2F1(A,B;C;X). C C Date: 2008, June 3 C Coded by: Ludek Klimes C C----------------------------------------------------------------------- C C Local storage locations: INTEGER I REAL AA,BB,CC,DD,XX,ERR1,RATIO,TERM,SUM C C....................................................................... C C Check of input values IF(A.LE.0..OR.B.LE.0..OR.C.LE.0..OR.X.LT.0..OR.X.GE.1. * .OR.ERR.LE.0.) THEN C HG-01 CALL ERROR('HG-01: Wrong input values in HG2F1') END IF C C Summing the Gauss series: AA=A BB=B CC=C DD=1. XX=X ERR1=ERR C Zero and first terms TERM=XX*AA*BB/CC SUM=1.+TERM C Second and subsequent terms DO 1 I=2,1000 AA=AA+1. BB=BB+1. CC=CC+1. DD=DD+1. RATIO=XX*AA*BB/CC/DD IF(RATIO.EQ.0.) THEN F=SUM RETURN END IF TERM=TERM*RATIO SUM=SUM+TERM IF(TERM.LT.ERR1/RATIO-ERR1) THEN F=SUM RETURN END IF 1 CONTINUE C C HG-02 CALL ERROR('HG-02: Convergence failure in HG2F1') RETURN END C C======================================================================= C C C SUBROUTINE HGF1(A,B1,C,X1,X2,ERR,F) REAL A,B1,C,X1,X2,ERR,F C C Subroutine designed to calculate hypergeometric function C F1(A,B1,1,C;X1,X2). C C Input: C A,B1,C..Parameters of hypergeometric function F1(A,B1,1,C;X1,X2). C The parameters are assumed to be positive. C X1,X2...Variables of hypergeometric function F1(A,B1,1,C;X1,X2). C The subroutine is designed for non-negative X1 and X2 C sufficiently smaller than 1. C ERR... Absolute error of the output value. C Output: C F... Value of hypergeometric function F1(A,B1,1,C;X1,X2). C C Date: 2008, June 3 C Coded by: Ludek Klimes C C----------------------------------------------------------------------- C C Local storage locations: INTEGER I1,I2 REAL AA1,BB1,CC1,DD1,AA2,CC2,XX1,XX2,ERR1 REAL RATIO1,RATIO2,TERM1,TERM2,SUM1,SUM2 C C....................................................................... C C Check of input values IF(A.LE.0..OR.B1.LE.0..OR.C.LE.0..OR.X1.LT.0..OR.X1.GE.1. * .OR.ERR.LE.0..OR.X2.LT.0..OR.X2.GE.1.) THEN C HG-03 CALL ERROR('HG-03: Wrong input values in HGF1') END IF IF(X1.EQ.0.) THEN CALL HG2F1(A,1.,C,X2,ERR,F) RETURN END IF IF(X2.EQ.0.) THEN CALL HG2F1(A,B1,C,X1,ERR,F) RETURN END IF C C Estimating the number of iterations with respect to X1 ERR1=ALOG((1.-X1)*ERR)/ALOG(X1) C Error for iterations with respect to X2 ERR1=ERR/(AMAX1(0.,ERR1)+1.) C C Summing the Gauss series: AA1=A BB1=B1 CC1=C DD1=1. XX1=X1 XX2=X2 TERM1=1. SUM1=0. C Loop over powers in X1 DO 11 I1=0,1000 AA2=AA1 CC2=CC1 C Zero and first terms with respect to X2 TERM2=TERM1*XX2*AA2/CC2 SUM2=TERM1+TERM2 C Second and subsequent terms with respect to X2 DO 2 I2=2,1000 AA2=AA2+1. CC2=CC2+1. RATIO2=XX2*AA2/CC2 IF(RATIO2.EQ.0.) THEN GO TO 10 END IF TERM2=TERM2*RATIO2 SUM2=SUM2+TERM2 IF(TERM2.LT.ERR1/RATIO2-ERR1) THEN GO TO 10 END IF 2 CONTINUE C HG-04 CALL ERROR('HG-04: Convergence failure in inner loop in HGF1') 10 CONTINUE SUM1=SUM1+SUM2 RATIO1=XX1*AA1*BB1/CC1/DD1 IF(RATIO1.EQ.0.) THEN F=SUM1 RETURN END IF IF(SUM2.LT.ERR1/RATIO1-ERR1) THEN F=SUM1 RETURN END IF C Values for the next iteration with respect to X1 TERM1=TERM1*RATIO1 AA1=AA1+1. BB1=BB1+1. CC1=CC1+1. DD1=DD1+1. 11 CONTINUE C C HG-05 CALL ERROR('HG-05: Convergence failure in outer loop in HGF1') RETURN END C C======================================================================= C C C SUBROUTINE HGF2(A,B1,C1,C2,X1,X2,ERR,F) REAL A,B1,C1,C2,X1,X2,ERR,F C C Subroutine designed to calculate hypergeometric function C F2(A,B1,1,C1,C2;X1,X2). C C Input: C A,B1,C1,C2... Parameters of hypergeometric function C F2(A,B1,1,C1,C2;X1,X2). C The parameters are assumed to be positive. C X1,X2...Variables of hypergeometric function C F2(A,B1,1,C1,C2;X1,X2). C The subroutine is designed for non-negative X1 and X2 C with X1+X2 sufficiently smaller than 1. C ERR... Absolute error of the output value. C Output: C F... Value of hypergeometric function F2(A,B1,1,C1,C2;X1,X2). C C Date: 2008, June 3 C Coded by: Ludek Klimes C C----------------------------------------------------------------------- C C Local storage locations: INTEGER I1,I2 REAL AA1,BB1,CC1,DD1,AA2,CC2,XX1,XX2,ERR1 REAL RATIO1,RATIO2,TERM1,TERM2,SUM1,SUM2 C C....................................................................... C C Check of input values IF(A.LE.0..OR.B1.LE.0..OR.C1.LE.0..OR.X1.LT.0..OR.X1+X2.GE.1. * .OR.ERR.LE.0..OR.C2.LE.0..OR.X2.LT.0.) THEN C HG-06 CALL ERROR('HG-06: Wrong input values in HGF2') END IF IF(X1.EQ.0.) THEN CALL HG2F1(A,1.,C2,X2,ERR,F) RETURN END IF IF(X2.EQ.0.) THEN CALL HG2F1(A,B1,C1,X1,ERR,F) RETURN END IF C C Estimating the number of iterations with respect to X1 ERR1=ALOG((1.-X1-X2)*ERR)/ALOG(X1+X2) C Error for iterations with respect to X2 ERR1=ERR/(AMAX1(0.,ERR1)+1.) C C Summing the Gauss series: AA1=A BB1=B1 CC1=C1 DD1=1. XX1=X1 XX2=X2 TERM1=1. SUM1=0. C Loop over powers in X1 DO 11 I1=0,1000 AA2=AA1 CC2=C2 C Zero and first terms with respect to X2 TERM2=TERM1*XX2*AA2/CC2 SUM2=TERM1+TERM2 C Second and subsequent terms with respect to X2 DO 2 I2=2,1000 AA2=AA2+1. CC2=CC2+1. RATIO2=XX2*AA2/CC2 IF(RATIO2.EQ.0.) THEN GO TO 10 END IF TERM2=TERM2*RATIO2 SUM2=SUM2+TERM2 IF(TERM2.LT.ERR1/RATIO2-ERR1) THEN GO TO 10 END IF 2 CONTINUE C HG-07 CALL ERROR('HG-07: Convergence failure in inner loop in HGF2') 10 CONTINUE SUM1=SUM1+SUM2 RATIO1=XX1*AA1*BB1/CC1/DD1 IF(RATIO1.EQ.0.) THEN F=SUM1 RETURN END IF IF(SUM2.LT.ERR1/RATIO1-ERR1) THEN F=SUM1 RETURN END IF C Values for the next iteration with respect to X1 TERM1=TERM1*RATIO1 AA1=AA1+1. BB1=BB1+1. CC1=CC1+1. DD1=DD1+1. 11 CONTINUE C C HG-08 CALL ERROR('HG-08: Convergence failure in outer loop in HGF2') RETURN END C C======================================================================= C C C SUBROUTINE HGFM2(A,B1,C1,X1,X2,ERR,F) REAL A,B1,C1,X1,X2,ERR,F C C Subroutine designed to calculate modified hypergeometric function C F~2(A,B1,1,C1,1;X1,X2). C C Input: C A,B1,C1... Parameters of modified hypergeometric function C F~2(A,B1,1,C1,1;X1,X2). C Parameters A and C1 are assumed to be positive, C parameter B1 is assumed to be greater than -1. C X1,X2...Variables of modified hypergeometric function C F~2(A,B1,1,C1,1;X1,X2). C The subroutine is designed for non-negative X1 and X2 C sufficiently smaller than 1. C ERR... Absolute error of the output value. C Output: C F... Value of hypergeometric function F~2(A,B1,1,C1,1;X1,X2). C C Date: 2008, June 3 C Coded by: Ludek Klimes C C----------------------------------------------------------------------- C C Local storage locations: INTEGER I1,I2 REAL AA1,BB1,CC1,DD1,AA2,DD2,XX1,XX2,ERR1 REAL RATIO1,RATIO2,TERM1,TERM2,SUM1,SUM2 C C....................................................................... C C Check of input values IF(A.LE.0..OR.B1.LE.-1..OR.C1.LE.0..OR.X1.LT.0..OR.X1.GE.1. * .OR.ERR.LE.0..OR.X2.LT.0..OR.X2.GE.1.) THEN C HG-09 CALL ERROR('HG-09: Wrong input values in HGFM2') END IF IF(X1.EQ.0.) THEN CALL HG2F1(A,1.,1.,X2,ERR,F) RETURN END IF IF(X2.EQ.0.) THEN CALL HG2F1(A,B1,C1,X1,ERR,F) RETURN END IF C C Estimating the number of iterations with respect to X1 ERR1=ALOG((1.-X1)*ERR)/ALOG(X1) C Error for iterations with respect to X2 ERR1=ERR/(AMAX1(0.,ERR1)+1.) C C Summing the Gauss series: AA1=A BB1=B1 CC1=C1 DD1=1. XX1=X1 XX2=X2 TERM1=1. SUM1=0. C Loop over powers in X1 DO 11 I1=0,1000 AA2=AA1 DD2=DD1 C Zero and first terms with respect to X2 TERM2=TERM1*XX2*AA2/DD2 SUM2=TERM1+TERM2 C Second and subsequent terms with respect to X2 DO 2 I2=2,1000 AA2=AA2+1. DD2=DD2+1. RATIO2=XX2*AA2/DD2 IF(RATIO2.EQ.0.) THEN GO TO 10 END IF TERM2=TERM2*RATIO2 SUM2=SUM2+TERM2 IF(ABS(TERM2).LT.ERR1/RATIO2-ERR1) THEN GO TO 10 END IF 2 CONTINUE C HG-10 CALL ERROR('HG-10: Convergence failure in inner loop in HGFM2') 10 CONTINUE SUM1=SUM1+SUM2 RATIO1=XX1*AA1*BB1/CC1/DD1 IF(RATIO1.EQ.0.) THEN F=SUM1 RETURN END IF IF(ABS(SUM2).LT.ERR1/RATIO1-ERR1) THEN F=SUM1 RETURN END IF C Values for the next iteration with respect to X1 TERM1=TERM1*RATIO1 AA1=AA1+1. BB1=BB1+1. CC1=CC1+1. DD1=DD1+1. 11 CONTINUE C C HG-11 CALL ERROR('HG-11: Convergence failure in outer loop in HGFM2') RETURN END C C======================================================================= C