HYP_2F1 Function

  RECURSIVE FUNCTION HYP_2F1(A,B,C,Z) RESULT(RES)
    !! CPC Michel Gauss hypergeometric function \({}_2F_1(a, b; c; z)\).
    !
    ! Calculation of F(z) for arbitrary z using previous routines
    !
    ! Firstly, it is checked if a,b and c are negative integers.
    ! If neither a nor b is negative integer but c is, 
    ! F(z) is undefined so that the program stops with an error message.
    ! If a and c are negative integers with c < a, 
    ! or b and c are negative integers with b < a, 
    ! or c is not negative integer integer but a or b is, 
    ! one is in the polynomial case.
    ! In this case, if |z| < |z/(z-1)| or z = 1, 
    ! hyp_PS_zero is used directly, as then |z| <= 2 
    ! and no instability arises with hyp_PS_zero 
    ! as long the degree of the polynomial is small (<= 10 typically).
    ! If not, one uses the transformation 
    ! F(z) = (1-z)^{-a} 2F1(a,c-b,c,z/(z-1)) if a is negative integer 
    ! or F(z) = (1-z)^{-b} 2F1(b,c-a,c,z/(z-1)) if b is negative integer 
    ! along with hyp_PS_zero.
    ! Indeed, 2F1(a,c-b,c,X) is a polynomial if a is negative integer, 
    ! and so is 2F1(b,c-a,c,X) if b is negative integer, 
    ! so that one has here |z/(z-1)| <= 2 
    ! and the stability of the method is the same 
    ! as for the |z| < |z/(z-1)| case.
    ! If one is in the non-polynomial case, one checks if z >= 1. 
    ! If it is, one is the cut of F(z) 
    ! so that z is replaced by z - 10^{-307}i.
    ! Then, using F(z) = 2F1(b,a,c,z) 
    ! and F(z) = (1-z)^{c-a-b} 2F1(c-a,c-b,c,z), 
    ! one replaces a,b,c parameters by combinations of them 
    ! so that Re(b-a) >= 0 and Re(c-a-b) >= 0.
    ! Exchanging a and b does not change convergence properties, 
    ! while having Re(c-a-b) >= 0 accelerates it 
    ! (In hyp_PS_zero, t[n] z^n ~ z^n/(n^{c-a-b}) for n -> +oo).
    ! If |1-z| < 1E-5, one uses hyp_PS_one 
    ! as the vicinity of the singular point z = 1 is treated properly.
    ! After that, one compares |z| and |z/(z-1)| 
    ! to R in {0.5,0.6,0.7,0.8,0.9}. 
    ! If one of them is smaller than R, 
    ! one uses hyp_PS_zero without transformation
    ! or with the transformation F(z) = (1-z)^{-a} 2F1(a,c-b,c,z/(z-1)).
    ! Then, if both of them are larger than 0.9, 
    ! one compares |1/z|, |(z-1)/z|, |1-z| and |1/(1-z)| 
    ! to R in {0.5,0.6,0.7,0.8,0.9}. 
    ! If one of them is found smaller than R, 
    ! with the condition that |c-b|oo < 5 for (z-1)/z transformation, 
    ! |a,b,c|oo < 5 for |1-z| transformation 
    ! and |a,c-b,c|oo < 5 for |1/(1-z)| transformation,
    ! the corresponding transformation is used. 
    ! If none of them was smaller than 0.9, 
    ! one is in the missing zones of transformation theory 
    ! so that the Taylor series of hyp_PS_complex_plane_rest is used.
    !
    ! Variables
    ! ---------
    ! a,b,c,z: a,b,c parameters and z argument of 2F1(a,b,c,z)
    ! Re_a,Re_b,Re_c,na,nb,nc,is_a_neg_int,is_b_neg_int,is_c_neg_int: 
    ! real parts of a,b,c, closest integers to a,b,c, 
    ! true if a,b,c is negative integers and false if not.
    ! zm1,z_over_zm1,z_shift: z-1, z/(z-1), z - 10^{-307}i in case z >= 1.
    ! ab_condition, cab_condition: true if Re(b-a) >= 0 and false if not, 
    ! true if Re(c-a-b) >= 0 and false if not.
    ! abs_zm1,abz_z,abs_z_inv,abs_z_over_zm1,abs_zm1_inv,abs_zm1_over_z: 
    ! |z-1|, |z|, |1/z|, |z/(z-1)|, |1/(z-1)|, |(z-1)/z|
    ! are_ac_small: true if |a|oo < 5 and |c|oo < 5, false if not.
    ! is_cmb_small: true if |c-b|oo < 5, false if not.
    ! are_abc_small: true if |a|oo < 5, |b|oo < 5 and |c|oo < 5, 
    ! false if not.
    ! are_a_cmb_c_small: true if |a|oo < 5, |c-b|oo < 5 and |c|oo < 5, 
    ! false if not.
    ! R_tab,R: table of radii {0.5,0.6,0.7,0.8,0.9}, one of these radii.
    ! res: returned result
    !----------------------------------------------------------------------

    COMPLEX(wp),INTENT(IN) :: A,B,C,Z
    INTEGER(i4) :: NA,NB,NC,I
    REAL(wp)     :: RE_A,RE_B,RE_C,ABS_Z,ABS_ZM1,ABS_Z_OVER_ZM1
    REAL(wp)     :: ABS_ZM1_OVER_Z,ABS_ZM1_INV,R_TABLE(1:5),R,ABS_Z_INV
    COMPLEX(wp)  :: RES,Z_SHIFT
    COMPLEX(wp)  :: Z_OVER_ZM1,ZM1
    LOGICAL      :: IS_A_NEG_INT,IS_B_NEG_INT,IS_C_NEG_INT
    LOGICAL      :: AB_CONDITION,CAB_CONDITION,ARE_A_CMB_C_SMALL
    LOGICAL      :: IS_CMB_SMALL,ARE_AC_SMALL,ARE_ABC_SMALL
    !
    RE_A=REAL(A,wp); RE_B=REAL(B,wp); RE_C=REAL(C,wp);
    NA=NINT(RE_A); NB=NINT(RE_B); NC=NINT(RE_C);
    IS_A_NEG_INT=A.EQ.NA.AND.NA.LE.0
    IS_B_NEG_INT=B.EQ.NB.AND.NB.LE.0
    IS_C_NEG_INT=C.EQ.NC.AND.NC.LE.0
    ZM1=Z-ONE
    IF(IS_C_NEG_INT) THEN
       ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1
       ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1)
       IF(IS_A_NEG_INT.AND.(NC.LT.NA)) THEN
          IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN
             RES=HYP_PS_ZERO(A,B,C,Z)
             RETURN
          ELSE
             RES=((-ZM1)**(-A))*HYP_PS_ZERO(A,C-B,C,Z_OVER_ZM1)
             RETURN
          ENDIF
       ELSE IF(IS_B_NEG_INT.AND.(NC.LT.NB)) THEN
          IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN
             RES=HYP_PS_ZERO(A,B,C,Z)
             RETURN
          ELSE
             RES=((-ZM1)**(-B))*HYP_PS_ZERO(B,C-A,C,Z_OVER_ZM1)
             RETURN
          ENDIF
       ELSE
          STOP '2F1 UNDEFINED'
       ENDIF
    ENDIF
    IF(IS_A_NEG_INT) THEN
       ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1
       ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1)
       IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN
          RES=HYP_PS_ZERO(A,B,C,Z)
          RETURN
       ELSE
          RES=((-ZM1)**(-A))*HYP_PS_ZERO(A,C-B,C,Z_OVER_ZM1)
          RETURN
       ENDIF
    ELSE IF(IS_B_NEG_INT) THEN
       ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1
       ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1)
       IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN
          RES=HYP_PS_ZERO(A,B,C,Z)
          RETURN
       ELSE
          RES=((-ZM1)**(-B))*HYP_PS_ZERO(B,C-A,C,Z_OVER_ZM1)
          RETURN
       ENDIF
    ENDIF
    IF((REAL(Z,wp).GE.ONE).AND.(AIMAG(Z).EQ.ZERO)) THEN
       Z_SHIFT=CMPLX(REAL(Z,wp),-1.0e-307_wp,wp)
       RES=HYP_2F1(A,B,C,Z_SHIFT)
       RETURN
    ENDIF
    AB_CONDITION = (RE_B.GE.RE_A)
    CAB_CONDITION = (RE_C.GE.RE_A + RE_B)
    IF ((.NOT.AB_CONDITION).OR.(.NOT.CAB_CONDITION)) THEN
       IF ((.NOT.AB_CONDITION).AND.(CAB_CONDITION)) THEN
          RES=HYP_2F1(B,A,C,Z)
          RETURN
       ELSE IF((.NOT.CAB_CONDITION).AND.(AB_CONDITION)) THEN
          RES=((-ZM1)**(C-A-B))*HYP_2F1(C-B,C-A,C,Z)
          RETURN
       ELSE
          RES=((-ZM1)**(C-A-B))*HYP_2F1(C-A,C-B,C,Z)
          RETURN 
       ENDIF
    ENDIF
    ABS_ZM1=ABS(ZM1)
    IF(ABS_ZM1.LT.1.0e-5_wp) THEN 
       RES=HYP_PS_ONE (A,B,C,-ZM1)
       RETURN
    ENDIF
    ABS_Z=ABS(Z); ABS_Z_OVER_ZM1=ABS_Z/ABS_ZM1; ABS_Z_INV=ONE/ABS_Z
    ABS_ZM1_OVER_Z=ONE/ABS_Z_OVER_ZM1; ABS_ZM1_INV=ONE/ABS_ZM1
    IS_CMB_SMALL = INF_NORM(C-B).LT.5.0_wp; 
    ARE_AC_SMALL = (INF_NORM(A).LT.5.0_wp).AND.(INF_NORM(C).LT.5.0_wp)
    ARE_ABC_SMALL = ARE_AC_SMALL.AND.(INF_NORM(B).LT.5.0_wp)
    ARE_A_CMB_C_SMALL = ARE_AC_SMALL.AND.IS_CMB_SMALL
    R_TABLE=[0.5_wp,0.6_wp,0.7_wp,0.8_wp,0.9_wp]
    DO I=1,5
       R=R_TABLE(I)
       IF(ABS_Z.LE.R) THEN 
          RES=HYP_PS_ZERO (A,B,C,Z)
          RETURN
       ENDIF
       IF(IS_CMB_SMALL.AND.(ABS_Z_OVER_ZM1.LE.R)) THEN
          RES=((-ZM1)**(-A))*HYP_PS_ZERO (A,C-B,C,Z/ZM1)
          RETURN
       ENDIF
    ENDDO
    DO I=1,5
       R=R_TABLE(I)
       IF(ABS_Z_INV.LE.R) THEN 
          RES=HYP_PS_INFINITY (A,B,C,Z)
          RETURN 
       ENDIF
       IF(IS_CMB_SMALL.AND.(ABS_ZM1_OVER_Z.LE.R)) THEN 
          RES=((-ZM1)**(-A))*HYP_PS_INFINITY (A,C-B,C,Z/ZM1)
          RETURN
       ENDIF
       IF(ARE_ABC_SMALL.AND.(ABS_ZM1.LE.R)) THEN 
          RES=HYP_PS_ONE (A,B,C,-ZM1)
          RETURN
       ENDIF
       IF(ARE_A_CMB_C_SMALL.AND.(ABS_ZM1_INV.LE.R)) THEN 
          RES=((-ZM1)**(-A))*HYP_PS_ONE (A,C-B,C,-ONE/ZM1)
          RETURN
       ENDIF
    ENDDO
    RES=HYP_PS_COMPLEX_PLANE_REST (A,B,C,Z)
    RETURN
  END FUNCTION HYP_2F1