      SUBROUTINE GRATIO(A,X,ANS,QANS,IOP3)
      REAL J, L, ACC0(3), E00(3), X00(3)
      REAL D0(22),D2(20),D4(12),D6(8),D8(5),D10(1)
C     ------------------------------------------------------------------
C        EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
C                      P(A,X) AND Q(A,X)
C
C     IOP3 SPECIFIES THE RELATIVE ERROR TOLERANCE ACC TO BE USED.
C     IF IOP3 = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
C     POSSIBLE (UP TO 12 DIGITS). IF IOP3 = 1 THEN 6 DIGIT ACCURACY
C     IS REQUESTED. OTHERWISE, 3 DIGIT ACCURACY IS REQUESTED.
C
C     E0 AND X0 DEPEND ON THE TOLERANCE ACC.
C     ------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS
C        NAVAL SURFACE WEAPONS CENTER
C        DAHLGREN, VIRGINIA
C        MARCH, 1984
C     -------------------
      DATA ACC0(1)/5.E-13/, ACC0(2)/5.E-7/, ACC0(3)/5.E-4/
      DATA E00(1)/1.25E-3/, E00(2)/3.85E-2/, E00(3)/.22/
      DATA X00(1)/31.0/, X00(2)/17.0/, X00(3)/9.7/
C     -------------------
C     ALOG10 = LN(10)
C     RT2PIN = 1/SQRT(2*PI)
C     RTPI   = SQRT(PI)
C     -------------------
      DATA ALOG10/2.30258509299405/
      DATA RT2PIN/.398942280401433/
      DATA RTPI  /1.77245385090552/
      DATA THIRD /.333333333333333/
C     -------------------
         DATA D0(2)  / -.833333333333333E-01/
         DATA D0(3)  / -.425925925925926E-01/
         DATA D0(4)  / -.272376543209877E-01/
         DATA D0(5)  / -.194775132275132E-01/
         DATA D0(6)  / -.148962007642563E-01/
         DATA D0(7)  / -.119154786400157E-01/
         DATA D0(8)  / -.984223052003723E-02/
         DATA D0(9)  / -.832809351218023E-02/
         DATA D0(10) / -.718034838506986E-02/
         DATA D0(11) / -.628441927210121E-02/
         DATA D0(12) / -.556825261788511E-02/
         DATA D0(13) / -.498444568441567E-02/
         DATA D0(14) / -.450063635725766E-02/
         DATA D0(15) / -.409403566946016E-02/
         DATA D0(16) / -.374816818889724E-02/
         DATA D0(17) / -.345085037697095E-02/
         DATA D0(18) / -.319289361180614E-02/
         DATA D0(19) / -.296724737735710E-02/
         DATA D0(20) / -.276841853775599E-02/
         DATA D0(21) /  .154857396777124E-01/
         DATA D0(22) /  .269935709445088E-02/
         DATA D2(1)  /  .185185185185185E-02/
         DATA D2(2)  / -.347222222222222E-02/
         DATA D2(3)  / -.380291005291005E-02/
         DATA D2(4)  / -.342904908877131E-02/
         DATA D2(5)  / -.298813198118754E-02/
         DATA D2(6)  / -.259725819983343E-02/
         DATA D2(7)  / -.227086554281877E-02/
         DATA D2(8)  / -.200161263181312E-02/
         DATA D2(9)  / -.177894274461660E-02/
         DATA D2(10) / -.159341392427589E-02/
         DATA D2(11) / -.143743058374835E-02/
         DATA D2(12) / -.130507556863837E-02/
         DATA D2(13) / -.119177327190732E-02/
         DATA D2(14) / -.109397653898372E-02/
         DATA D2(15) / -.100891680192332E-02/
         DATA D2(16) / -.934414630830751E-03/
         DATA D2(17) / -.868738608276550E-03/
         DATA D2(18) / -.810500575064140E-03/
         DATA D2(19) /  .534748475828742E-01/
         DATA D2(20) / -.214654105710934E-01/
         DATA D4(1)  /  .413359788359788E-02/
         DATA D4(2)  /  .268132716049383E-02/
         DATA D4(3)  /  .166538065843621E-02/
         DATA D4(4)  /  .103376307441701E-02/
         DATA D4(5)  /  .638902257745493E-03/
         DATA D4(6)  /  .386095165779183E-03/
         DATA D4(7)  /  .220253534240928E-03/
         DATA D4(8)  /  .109183361550253E-03/
         DATA D4(9)  /  .335805190004565E-04/
         DATA D4(10) / -.184745824613574E-04/
         DATA D4(11) / -.545519921343409E-04/
         DATA D4(12) / -.795842560994557E-04/
         DATA D6(1)  / -.649434156378601E-03/
         DATA D6(2)  /  .229472093621399E-03/
         DATA D6(3)  /  .545680192269055E-03/
         DATA D6(4)  /  .625133202086059E-03/
         DATA D6(5)  /  .608954623450347E-03/
         DATA D6(6)  /  .557237674593796E-03/
         DATA D6(7)  /  .495639378848119E-03/
         DATA D6(8)  /  .434915394155868E-03/
         DATA D8(1)  / -.861888290916712E-03/
         DATA D8(2)  / -.784039221720067E-03/
         DATA D8(3)  / -.560418887543212E-03/
         DATA D8(4)  / -.350369879167589E-03/
         DATA D8(5)  / -.187647900119735E-03/
         DATA D10(1) /  .336798553366358E-03/
C     -------------------
C     ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
C            FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
C
         E = SPMPAR(1)
C     -------------------
      IF (A .LT. 0.0 .OR. X .LT. 0.0) GO TO 400
      IF (A .EQ. 0.0 .AND. X .EQ. 0.0) GO TO 400
      IF (A*X .EQ. 0.0) GO TO 331
C
   10 IOP = IOP3 + 1
      IF (IOP .NE. 1 .AND. IOP .NE. 2) IOP = 3
      ACC = AMAX1(ACC0(IOP),E)
      E0 = E00(IOP)
      X0 = X00(IOP)
C
C       SELECT THE APPROPRIATE ALGORITHM
C
   20 IF (A .GE. 1.0) GO TO 30
      IF (A .EQ. 0.5) GO TO 320
      IF (X .LT. 1.5) GO TO 110
      T1 = A*ALOG(X) - X
      U = A*EXP(T1)
      IF (U .EQ. 0.0) GO TO 310
      R = U/GAMMA(A + 1.0)
      GO TO 170
C
   30 IF (A .GE. 30.0) GO TO 40
      IF (A .GT. X .OR. X .GE. X0) GO TO 31
      TWOA = A + A
      M = INT(TWOA)
      IF (TWOA .NE. FLOAT(M)) GO TO 31
      I = M/2
      IF (A .EQ. FLOAT(I)) GO TO 140
      GO TO 150
   31 T1 = A*ALOG(X) - X
      R = EXP(T1)/GAMMA(A)
      GO TO 50
C
   40 L = X/A
      IF (L .EQ. 0.0) GO TO 300
      S = 0.5 + (0.5 - L)
      Y = RLOG(L)
      IF (Y .GE. 700.0/A) GO TO 330
      Y = A*Y
      RTA = SQRT(A)
      IF (ABS(S) .LE. E0/RTA) GO TO 180
      IF (IOP .EQ. 1 .AND. ABS(S) .LE. 1.E-3) GO TO 185
      IF (A .GE. 75.0 .AND. ABS(S) .LE. THIRD) GO TO 200
      T = 1.0/(A*A)
      T1 = (((0.75*T - 1.0)*T + 3.5)*T - 105.0)/(A*1260.0)
      T1 = T1 - Y
      R = RT2PIN*RTA*EXP(T1)
C
   50 IF (R .EQ. 0.0) GO TO 331
      IF (X .LE. AMAX1(A,ALOG10)) GO TO 100
      IF (X .LT. X0) GO TO 170
      GO TO 130
C
C                 FORMULA (84)
C
  100 HACC = 0.5*ACC
      SUM = 1.0
      T = 1.0
      C = A
  101 C = C + 1.0
      T = T*(X/C)
      SUM = SUM + T
      IF (T .GT. HACC) GO TO 101
      ANS = (R/A)*SUM
      QANS = 1.0 - ANS
      RETURN
C
C                 FORMULA (88)
C
  110 T = X
      C = 1.0
      SUM = T/(A + C)
      TOL = (0.1*ACC)*SUM
  111 C = C + 1.0
      T = -T*(X/C)
      AK = T/(A + C)
      SUM = SUM + AK
      IF (ABS(AK) .GT. TOL) GO TO 111
      J = A*(SUM - 0.5*AK)
C
      T1 = A*ALOG(X)
      U = EXP(T1)
      H = RGAM(A,ACC)
      G = 1.0 + H
      ANS = U*G*(0.5 + (0.5 - J))
      QANS = 1.0 - ANS
      IF (ANS .LE. 0.9) RETURN
C
C                 FORMULA (90)
C
      TOL = 0.1*ACC
      T = 0.5*T1
      C = 2.0
      SUM = T
  120 C = C + 1.0
      T = T*(T1/C)
      SUM = SUM + T
      IF (ABS(T) .GT. TOL) GO TO 120
      L = T1*(1.0 + SUM)
C
      QANS = (U*J - L)*G - H
      IF (QANS .LE. 0.0) GO TO 310
      RETURN
C
C                 FORMULA (85)
C
  130 SUM = 1.0
      T = 1.0
      C = A
  131 C = C - 1.0
      T = T*(C/X)
      SUM = SUM + T
      IF (ABS(T) .GE. ACC) GO TO 131
      QANS = (R/X)*SUM
      ANS = 1.0 - QANS
      RETURN
C
C                 FORMULA (86)
C
  140 SUM = EXP(-X)
      T = SUM
      N = 1
      C = 0.0
      GO TO 160
C
  150 RTX = SQRT(X)
      SUM = ERFC1(0,RTX)
      T = EXP(-X)/(RTPI*RTX)
      N = 0
      C = -0.5
C
  160 IF (N .EQ. I) GO TO 161
      N = N + 1
      C = C + 1.0
      T = (X*T)/C
      SUM = SUM + T
      GO TO 160
  161 QANS = SUM
      ANS = 1.0 - QANS
      RETURN
C
C                 FORMULA (87)
C
  170 A2NM1 = 1.0
      A2N = 1.0
      B2NM1 = X
      B2N = X + (1.0 - A)
      C = 1.0
  171 A2NM1 = X*A2N + C*A2NM1
      B2NM1 = X*B2N + C*B2NM1
      AM0 = A2NM1/B2NM1
      C = C + 1.0
      CMA = C - A
      A2N = A2NM1 + CMA*A2N
      B2N = B2NM1 + CMA*B2N
      AN0 = A2N/B2N
      IF (ABS(AN0 - AM0) .GE. ACC*AN0) GO TO 171
C
      QANS = R*AN0
      ANS = 1.0 - QANS
      RETURN
C
C             FORMULAS (95) AND (96)
C
  180 IF (A*E*E .GT. 3.28E-3) GO TO 400
      C = 0.5 + (0.5 - Y)
      W = (0.5 - SQRT(Y)*(1.0 - Y/3.0)/RTPI)/C
      T = 1.0/A
      IF (IOP - 2) 190,191,192
C
  185 IF (ABS(S) .LE. 2.0*E .AND. A*E*E .GT. 3.28E-3) GO TO 400
      C = EXP(-Y)
      W = 0.5*ERFC1(1,SQRT(Y))
      T = 1.0/A
C
  190 Q3 = ((((-D10(1)*T + (D8(2)*S + D8(1)))*T - (D6(2)*S
     1        + D6(1)))*T + ((D4(3)*S + D4(2))*S + D4(1)))*T
     2        - ((D2(3)*S + D2(2))*S + D2(1)))*T + ((D0(4)*S
     3        + D0(3))*S + D0(2))*S - THIRD
      GO TO 240
  191 Q3 = (D4(1)*T - D2(1))*T + (D0(3)*S + D0(2))*S - THIRD
      GO TO 240
  192 Q3 = D0(2)*S - THIRD
      GO TO 240
C
C                 FORMULA (97)
C
  200 IF (ABS(S) .LE. 2.0*E .AND. A*E*E .GT. 3.28E-3) GO TO 400
      C = EXP(-Y)
      W = 0.5*ERFC1(1,SQRT(Y))
      T = 1.0/A
      S1 = S
      IF (IOP - 2) 210,220,230
C
  210 C0 = (((((((((((D0(22) * S1 + D0(21)) * S1 + D0(20))
     1     * S1 + D0(19)) * S1 + D0(18)) * S1 + D0(17))
     2     * S1 + D0(16)) * S1 + D0(15)) * S1 + D0(14))
     3     * S1 + D0(13)) * S1 + D0(12)) * S1 + D0(11))
     4     * S1
      C0 = (((((((((C0 + D0(10)) * S1 + D0(9)) * S1
     1     + D0(8)) * S1 + D0(7)) * S1 + D0(6)) * S1
     2     + D0(5)) * S1 + D0(4)) * S1 + D0(3)) * S1
     3     + D0(2)) * S1 - THIRD
      C2 = (((((((((((D2(20) * S1 + D2(19)) * S1 + D2(18))
     1     * S1 + D2(17)) * S1 + D2(16)) * S1 + D2(15))
     2     * S1 + D2(14)) * S1 + D2(13)) * S1 + D2(12))
     3     * S1 + D2(11)) * S1 + D2(10)) * S1 + D2(9))
     4     * S1
      C2 = (((((((C2 + D2(8)) * S1 + D2(7)) * S1 + D2(6))
     1     * S1 + D2(5)) * S1 + D2(4)) * S1 + D2(3))
     2     * S1 + D2(2)) * S1 + D2(1)
      C4 = ((((((((((D4(12) * S1 + D4(11)) * S1 + D4(10))
     1     * S1 + D4(9)) * S1 + D4(8)) * S1 + D4(7))
     2     * S1 + D4(6)) * S1 + D4(5)) * S1 + D4(4))
     3     * S1 + D4(3)) * S1 + D4(2)) * S1 + D4(1)
      C6 = ((((((D6(8) * S1 + D6(7)) * S1 + D6(6))
     1     * S1 + D6(5)) * S1 + D6(4)) * S1 + D6(3))
     2     * S1 + D6(2)) * S1 + D6(1)
      C8 = (((D8(5) * S1 + D8(4)) * S1 + D8(3)) * S1
     1     + D8(2)) * S1 + D8(1)
      C10 = D10(1)
      Q3 = ((((-C10 * T + C8) * T - C6) * T + C4)
     1     * T - C2) * T + C0
      GO TO 240
C
  220 C0 = (((((((D0(9) * S1 + D0(8)) * S1 + D0(7))
     1     * S1 + D0(6)) * S1 + D0(5)) * S1 + D0(4))
     2     * S1 + D0(3)) * S1 + D0(2)) * S1 - THIRD
      C2 = ((D2(4) * S1 + D2(3)) * S1 + D2(2)) * S1 + D2(1)
      C4 = D4(1)
      Q3 = (C4 * T - C2) * T + C0
      GO TO 240
C
  230 Q3 = ((D0(4) * S1 + D0(3)) * S1 + D0(2)) * S1 - THIRD
C
  240 IF (L .LT. 1.0) GO TO 241
      QANS = C*(W + RT2PIN*Q3/RTA)
      ANS = 1.0 - QANS
      RETURN
  241 ANS = C*(W - RT2PIN*Q3/RTA)
      QANS = 1.0 - ANS
      RETURN
C
C                SPECIAL CASES
C
  300 ANS = 0.0
      QANS = 1.0
      RETURN
C
  310 ANS = 1.0
      QANS = 0.0
      RETURN
C
  320 IF (X .GE. 0.25) GO TO 321
      ANS = ERF(SQRT(X))
      QANS = 1.0 - ANS
      RETURN
  321 QANS = ERFC1(0,SQRT(X))
      ANS = 1.0 - QANS
      RETURN
C
  330 IF (ABS(S) .LE. 2.0*E) GO TO 400
  331 IF (X .LE. A) GO TO 300
      GO TO 310
C
C                ERROR RETURN
C
  400 ANS = 2.0
      RETURN
      END
