         SUBROUTINE PKILL(R,S1,S2,H,K,P)
C        ***************************************************************
C        ELLIPTICAL COVERAGE FUNCTION
C        ***************************************************************
         REAL K, K8, K9, N, X(43), Y(43)
C        ---------------------
C        RT2 = SQRT(2)
C        RTPINV = 1/SQRT(PI)
C        ---------------------
         DATA RT2/1.4142135623731/
         DATA RTPINV/.56418958354776/
C        ---------------------
C        B = A1**2
C        C = 2*(1 - EPS1)
C        ---------------------
         DATA A/4.892/, A1/5.387/, A2/3.8775/, A3/5.16/
         DATA ERR/1.E-6/, EPS1/1.04E-8/
         DATA B/29.019769/, C/1.9999999792/
C        ---------------------
         DATA X(1)  /.238619186083197/
         DATA X(2)  /.661209386466265/
         DATA X(3)  /.932469514203152/
         DATA X(4)  /.183434642495650/
         DATA X(5)  /.525532409916329/
         DATA X(6)  /.796666477413627/
         DATA X(7)  /.960289856497536/
         DATA X(8)  /.125233408511469/
         DATA X(9)  /.367831498998180/
         DATA X(10) /.587317954286617/
         DATA X(11) /.769902674194305/
         DATA X(12) /.904117256370475/
         DATA X(13) /.981560634246719/
         DATA X(14) /9.50125098376374E-2/
         DATA X(15) /.281603550779259/
         DATA X(16) /.458016777657227/
         DATA X(17) /.617876244402644/
         DATA X(18) /.755404408355003/
         DATA X(19) /.865631202387831/
         DATA X(20) /.944575023073232/
         DATA X(21) /.989400934991649/
         DATA X(22) /7.65265211334973E-2/
         DATA X(23) /.227785851141645/
         DATA X(24) /.373706088715420/
         DATA X(25) /.510867001950827/
         DATA X(26) /.636053680726515/
         DATA X(27) /.746331906460151/
         DATA X(28) /.839116971822219/
         DATA X(29) /.912234428251326/
         DATA X(30) /.963971927277914/
         DATA X(31) /.993128599185095/
         DATA X(32) /6.40568928626056E-2/
         DATA X(33) /.191118867473616/
         DATA X(34) /.315042679696163/
         DATA X(35) /.433793507626045/
         DATA X(36) /.545421471388840/
         DATA X(37) /.648093651936976/
         DATA X(38) /.740124191578554/
         DATA X(39) /.820001985973903/
         DATA X(40) /.886415527004401/
         DATA X(41) /.938274552002733/
         DATA X(42) /.974728555971309/
         DATA X(43) /.995187219997021/
C        ---------------------
         DATA Y(1)  /.467913934572691/
         DATA Y(2)  /.360761573048139/
         DATA Y(3)  /.171324492379170/
         DATA Y(4)  /.362683783378362/
         DATA Y(5)  /.313706645877887/
         DATA Y(6)  /.222381034453374/
         DATA Y(7)  /.101228536290376/
         DATA Y(8)  /.249147045813403/
         DATA Y(9)  /.233492536538355/
         DATA Y(10) /.203167426723066/
         DATA Y(11) /.160078328543346/
         DATA Y(12) /.106939325995318/
         DATA Y(13) /4.7175336386512E-2/
         DATA Y(14) /.189450610455068/
         DATA Y(15) /.182603415045924/
         DATA Y(16) /.169156519395003/
         DATA Y(17) /.149595988815577/
         DATA Y(18) /.124628971255534/
         DATA Y(19) /9.51585116824928E-2/
         DATA Y(20) /6.22535239386479E-2/
         DATA Y(21) /2.71524594117541E-2/
         DATA Y(22) /.152753387131726/
         DATA Y(23) /.149172986472604/
         DATA Y(24) /.142096109318382/
         DATA Y(25) /.131688638449177/
         DATA Y(26) /.118194531961518/
         DATA Y(27) /.101930119817240/
         DATA Y(28) /8.32767415767047E-2/
         DATA Y(29) /6.26720483341091E-2/
         DATA Y(30) /4.06014298003869E-2/
         DATA Y(31) /1.76140071391521E-2/
         DATA Y(32) /.127938195346752/
         DATA Y(33) /.125837456346828/
         DATA Y(34) /.121670472927803/
         DATA Y(35) /.115505668053726/
         DATA Y(36) /.107444270115966/
         DATA Y(37) /9.76186521041139E-2/
         DATA Y(38) /8.61901615319533E-2/
         DATA Y(39) /7.334648141108031E-2/
         DATA Y(40) /5.92985849154368E-2/
         DATA Y(41) /4.42774388174198E-2/
         DATA Y(42) /2.85313886289337E-2/
         DATA Y(43) /1.23412297999872E-2/
C
C
C                  CHECK IF P=0 OR P=1
C
         P = 0.0
         TOL = ERR*S1*S2
         RR = R*R
         IF (RR .LE. TOL) RETURN
C
         H2 = H*H + K*K
         D  = AMAX1(S1,S2)
         T  = R - A1*D
         IF (T .LT. 0.0) GO TO 10
              IF (T*T .LT. H2) GO TO 10
                   P = 1.0
                   RETURN
C
   10    H8 = ABS(H)
         K8 = ABS(K)
         IF (R + A*S1 .LE. H8) RETURN
         IF (R + A*S2 .LE. K8) RETURN
         S0 = SQRT(H2)
         IF (S1 .NE. S2) GO TO 20
              H8 = S0
              K8 = 0.0
              IF (R + A*S1 .LE. H8) RETURN
C
   20    IF (S0 .LE. R) GO TO 50
              D = (S0 - R)/D
              IF (RR*EXP(-0.5*D*D) .LE. TOL) RETURN
              IF (S0 .LT. R + A1*AMIN1(S1,S2)) GO TO 50
C
         IF (H8*K8 .EQ. 0.0) GO TO 50
         H9 = H8/S1
         K9 = K8/S2
         D = H9*H9 + K9*K9
         IF (D .LE. B) GO TO 50
C
         Q = (S2/S1)**2
         F = Q*H9*H9 + K9*K9
         Z1 = (F/D)*(R/S2)**2
         Z = D - Z1 - B
         IF (Z .LE. 0.0) GO TO 30
         IF (Z*Z .GE. 4.0*B*Z1) RETURN
C
   30    Z1 = B*S1*S1*(H*H + Q*K*K)/H2
         Z = H2 - RR - Z1
         IF (Z .LE. 0.0) GO TO 50
         IF (Z*Z .GE. 4.0*RR*Z1) RETURN
C
C             FIND THE LIMITS OF INTEGRATION
C
   50    Z  = K8 + A3 * S2
         U  = K8 - A3 * S2
         S0 = S1
         S  = S2
         I  = 0
C
   60    Z  = R - Z
         E0 = 0.0
         IF (Z .GT. 0.0) E0 = SQRT(Z/R)
         IF (U .GE. 0.0) GO TO 61
              E  = 1.0 - E0
              H5 = 0.0
              GO TO 70
   61    E  = SQRT(1.0 - U/R) - E0
         H5 = 1.0
C
   70    IF (I .NE. 0) GO TO 80
              I  = 1
              T  = E0
              F  = E
              H6 = H5
              Z  = H8 + A3 * S1
              U  = H8 - A3 * S1
              GO TO 60
C
   80    IF (F .GE. E) GO TO 100
              E0 = T
              E  = F
              T  = H8
              H8 = K8
              K8 = T
              H5 = H6
              S0 = S2
              S  = S1
C
C        SELECT THE APPROPRIATE INTEGRATION PROCESS
C
  100    E = 0.5*E
         N = E*R*(0.34/S0 + 1.0/(.025*ABS(R-K8) + 5.0*S))
         IF (N .LT. 2.75) GO TO 101
              J  = 31
              N1 = 12
              GO TO 110
  101    IF (N .LT. 1.35) GO TO 102
              J  = 21
              N1 = 10
              GO TO 110
  102    IF (N .LT. 0.75) GO TO 103
              J  = 13
              N1 =  8
              GO TO 110
  103    IF (N .LT. 0.35) GO TO 104
              J  =  7
              N1 =  6
              GO TO 110
  104    IF (N .LT. 0.15) GO TO 105
              J  =  3
              N1 =  4
              GO TO 110
  105    J  = 0
         N1 = 3
C
C              PERFORM GAUSSIAN INTEGRATION
C
  110    M  = N1 + N1
         R8 = R/(RT2*S0)
         H8 = H8/(RT2*S0)
         RR = R/(RT2*S)
         K8 = K8/(RT2*S)
         Q  = RTPINV*E*R8
         I  = -N1
         Z  = 0.0
         SUM = 0.0
         IF (K8 .EQ. 0.0) GO TO 200
C
         DO 170 L = 1,M
         IF (I .EQ. 0) I = 1
              II = J + IABS(I)
              T  = E*(1.0 + SIGN(X(II),FLOAT(I))) + E0
              TT = T*T
              D  = R8*(1.0 - TT)
              U  = H8 - D
              F  = EXP(-U*U)
C
              IF (H8 .NE. 0.0) GO TO 120
                   F = F + F
                   GO TO 130
  120         IF (H5 .NE. 0.0) GO TO 130
                   U = H8 + D
                   F = F + EXP(-U*U)
C
  130         IF (Z .EQ. 0.0) GO TO 140
                   P = C
                   GO TO 161
C
  140         D = RR*T*SQRT(2.0 - TT)
              U = K8 - D
              IF (ABS(U) .LT. A2) GO TO 150
              IF (U .GT. 0.0) GO TO 170
                   P = C
                   Z = 1.0
                   GO TO 161
C
  150         P = ERFC2(U)
              U = K8 + D
              IF (U .LT. A2) GO TO 160
                   P = P - EPS1
                   GO TO 161
C
  160         P = P - ERFC2(U)
  161         SUM = SUM + F*P*T*Y(II)
  170    I = I + 1
         P = Q*SUM
         RETURN
C
  200    DO 250 L = 1,M
         IF (I .EQ. 0) I = 1
              II = J + IABS(I)
              T  = E*(1.0 + SIGN(X(II),FLOAT(I))) + E0
              TT = T*T
              D  = R8*(1.0 - TT)
              U  = H8 - D
              F  = EXP(-U*U)
C
              IF (H8 .NE. 0.0) GO TO 210
                   F = F + F
                   GO TO 220
  210         IF (H5 .NE. 0.0) GO TO 220
                   U = H8 + D
                   F = F + EXP(-U*U)
C
  220         IF (Z .EQ. 0.0) GO TO 230
                   P = C
                   GO TO 241
C
  230         U = RR*T*SQRT(2.0 - TT)
              IF (U .LT. A2) GO TO 240
                   P = C
                   Z = 1.0
                   GO TO 241
C
  240         P = 2.0*(1.0 - ERFC2(U))
  241         SUM = SUM + F*P*T*Y(II)
  250    I = I + 1
         P = Q*SUM
         RETURN
         END
