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