      SUBROUTINE CIRCV(R,D,I,P,IERR)
C     ******************************************************************
C     CIRCULAR COVERAGE AND CIRCULAR ERROR FUNCTIONS
C     ******************************************************************
      REAL M,N
C     ------------------
C     RT2INV = 1/SQRT(2)
C     RTPINV = 1/SQRT(PI)
C     ------------------
      DATA RT2INV/.70710678118655/
      DATA RTPINV/.56418958354776/
C     ------------------
      DATA EPS/5.E-7/, TOL/5.E-4/
C
      IF (R .LT. 0.0) GO TO 200
      IF (D .LT. 0.0) GO TO 210
      IERR = 0
      P = 0.0
      IF (I .NE. 0) GO TO 30
C
C                    EVALUATION OF V(R,D)
C
      IF (D .GT. 1.0) GO TO 210
      IF (R .EQ. 0.0) RETURN
      IF (R .LT. 5.386773) GO TO 10
      P = 1.0
      RETURN
C
   10 IF (D .NE. 0.0) GO TO 20
      A = RT2INV*R
      P = ERF0(A)
      RETURN
C
   20 T = 0.5/D
      RR = R*(1.0 + D)*T
      DD = R*(1.0 - D)*T
      DIFF = R
      A = R
      GO TO 40
C
C                    EVALUATION OF P(R,D)
C
   30 IF (R .EQ. 0.0) RETURN
      RR = R
      DD = D
      DIFF = R - D
      A = ABS(DIFF)
      IF (A .LT. 5.386773) GO TO 40
      IF (DIFF .GE. 0.0) P = 1.0
      RETURN
C
   40 T = RR*DD
      T3 = 0.5*RR*RR
      B = 0.5*DD*DD
      N = 0.0
      IF (T .GT. 7.0) GO TO 70
C
C     EVALUATION OF P(RR,DD) WHEN ABS(RR-DD) .LT. 5.386773
C                     AND RR*DD .LE. 7.0
C
      T1 = RT2INV*T - 1.0
      T2 = T3*B
      S0 = EXP(-T3 - B)
      S1 = EXP(-B)
      IF (T3 .GT. TOL) GO TO 50
          S1 = S1*T3
          GO TO 51
   50 S1 = S1 - S0
   51 S2 = S0
      T0 = S1
C
   60 N = N + 1.0
      M = 1.0/N
      S0 = T2*M*M*S0
      T0 = B*M*T0 - S0
      S1 = S1 + T0
      S2 = S2 + S0
      IF (T1.GE.N .OR. T0.GT.EPS) GO TO 60
      P = S1
      IF (I .EQ. 0) P = ABS(P + P + S2 - 1.0)
      RETURN
C
C     EVALUATION OF P(RR,DD) WHEN ABS(RR-DD) .LT. 5.386773
C                     AND RR*DD .GT. 7.0
C
   70 A = RT2INV*A
      S1 = 0.5*DIFF*DIFF
      S2 = EXP(-S1)
      E = ERFC0(A,S1,S2)
      T1 = 2.0*ABS(T3 - B)
      T3 = 0.5/T
      T2 = SQRT(T3)
C
      S0 = RTPINV*T2*S2
      T0 = (RR + DD)*RT2INV*T2*E
      T2 = S1*T3
      T3 = 0.5*T3
      S1 = T0
      S2 = S0
C
   80 N = N + 2.0
      M = N - 1.0
      A = M/N
      S0 = A*T3*S0
      T0 = T1*S0 - T2*A*T0
      S0 = M*S0
      S1 = S1 + T0
      S2 = S2 + S0
      IF (T0 .GT. EPS) GO TO 80
C
   90 IF (S0 .LE. EPS) GO TO 100
      N = N + 2.0
      M = N - 1.0
      S0 = M*M*T3*S0/N
      S2 = S2 + S0
      GO TO 90
C
  100 IF (DIFF) 101,102,103
  101 P = 0.5*ABS(S1 - S2)
      RETURN
  102 P = 0.5*ABS(1.0 - S2)
      RETURN
  103 P = 0.5*ABS(2.0 - S2 - S1)
      IF (I .EQ. 0) P = ABS(P + P + S2 - 1.0)
      RETURN
C
C                       ERROR RETURN
C
  200 IERR = 1
      P = -1.0
      RETURN
  210 IERR = 2
      P = -1.0
      RETURN
      END
