1      SUBROUTINE DRCVAL (X, Y, RC, IERR)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>>   2001-07-16 DRCVAL Krogh  Change -1.0 to -1.d0.
6c>>   1996-03-30 DRCVAL Krogh  Added external statement.
7C>>   1995-11-17 DRCVAL Krogh  Converted SFTRAN to Fortran 77.
8C>>   1994-10-19 DRCVAL Krogh  Changes to use M77CON
9C>>   1991-10-31 DRCVAL WV Snyder JPL Incorporate changes from Carlson
10C>>   1990-12-20 DRCVAL WV Snyder JPL Convert from NSWC for Math 77.
11C
12C     THIS SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL
13C
14C     RC(X,Y) = INTEGRAL FROM ZERO TO INFINITY OF
15C
16C                         -1/2     -1
17C               (1/2)(T+X)    (T+Y)  DT,
18C
19C     WHERE X IS NONNEGATIVE AND Y IS NONZERO.  IF Y IS NEGATIVE,
20C     THE CAUCHY PRINCIPAL VALUE IS COMPUTED BY USING A PRELIMI-
21C     NARY TRANSFORMATION TO MAKE Y POSITIVE; SEE EQUATION (2.12)
22C     OF THE SECOND REFERENCE BELOW.  WHEN Y IS POSITIVE, THE
23C     DUPULICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
24C     NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
25C     SERIES TO FIFTH ORDER.  LOGARITHMIC, INVERSE CIRCULAR, AND
26C     INVERSE HYPERBOLIC FUNCTIONS ARE EXPRESSED IN TERMS OF RC
27C     BY EQUATIONS (4.9)-(4.13) OF THE SECOND REFERENCE BELOW.
28C
29C     REFERENCES: B. C. CARLSON AND E. M. NOTIS, ALGORITHMS FOR
30C     INCOMPLETE ELLIPTIC INTEGRALS, ACM TRANSACTIONS ON MATHEMA-
31C     TICAL SOFTWARE, 7 (1981), 398-403;
32C     B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
33C     NUMER. MATH. 33 (1979), 1-16.
34C     AUTHORS: B. C. CARLSON AND ELAINE M. NOTIS, AMES LABORATORY-
35C     DOE, IOWA STATE UNIVERSITY, AMES, IA 50011, AND R. L. PEXTON,
36C     LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE, CA 94550.
37C     AUG. 1, 1979, REVISED SEPT. 1, 1987.
38C
39C     CHECK VALUES: RC(0,1/4) = RC(1/16,1/8) = PI,
40C                   RC(9/4,2) = LN(2),
41C                   RC(1/4,-2) = LN(2)/3.
42C     CHECK BY ADDITION THEOREM: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z),
43C     WHERE X, Y, AND Z ARE POSITIVE AND  X * Y = Z * Z.
44C
45C    *****     Formal Arguments     ***********************************
46C
47C          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
48C     LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
49C     UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
50C
51C     INPUT ...
52C
53C     X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
54C
55C     OUTPUT ...
56C
57C     RC IS THE VALUE OF THE INTEGRAL.
58C
59C     IERR IS THE RETURN ERROR CODE.
60C          0  FOR NORMAL COMPLETION OF THE SUBROUTINE.
61C          1  X IS NEGATIVE, OR Y = 0.
62C          2  X+ABS(Y) IS TOO SMALL.
63C          3  X OR ABS(Y) IS TOO LARGE, OR X + ABS(Y) IS TOO LARGE.
64C          4  Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25
65C
66      DOUBLE PRECISION X, Y, RC
67      INTEGER IERR
68C
69C--D replaces "?": ?ERM1, ?ERV1, ?RCVAL
70C
71C     *****     External References     ********************************
72C
73      EXTERNAL D1MACH
74      DOUBLE PRECISION D1MACH
75C
76C     *****     Local Variables     ************************************
77C
78      DOUBLE PRECISION C1,C2,C3,C4,ERRTOL,LAMDA,LOLIM
79      DOUBLE PRECISION MU,SN,UPLIM,XLU225,XN,Y2236L,YN,W
80      PARAMETER (C1 = 3.0d0 / 10.0d0)
81      PARAMETER (C2 = 1.0d0 / 7.0d0)
82      PARAMETER (C3 = 3.0d0 / 8.0d0)
83      PARAMETER (C4 = 9.0d0 / 22.0d0)
84      SAVE ERRTOL, LOLIM, UPLIM, XLU225, Y2236L
85      DATA LOLIM /-1.0d0/
86C
87C     ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
88C     RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
89C     16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
90C
91C     SAMPLE CHOICES   ERRTOL   RELATIVE TRUNCATION
92C                               ERROR LESS THAN
93C                      1.D-3    2.D-17
94C                      3.D-3    2.D-14
95C                      1.D-2    2.D-11
96C                      3.D-2    2.D-8
97C                      1.D-1    2.D-5
98C
99C     We could put in a Newton iteration to solve for ERRTOL in
100C     terms of D1MACH(4), but it seems good enough to put
101C     ERRTOL = (D1MACH(4)/16.0)**(1.0/6.0)
102C
103c ----------------------------------------------------------------------
104C          WARNING. CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
105C          EXPENSE OF ROBUSTNESS.
106c ----------------------------------------------------------------------
107C
108      IF (X .LT. 0.0d0 .OR. Y .EQ. 0.0d0) THEN
109         RC = 0.0d0
110         CALL DERM1 ('DRCVAL',1,0,'X < 0 or Y = 0','X',X,',')
111         CALL DERV1 ('Y',Y,'.')
112         IERR = 1
113         RETURN
114      END IF
115      IF (LOLIM .LT. 0.0d0) THEN
116         LOLIM = 5.0d0 * D1MACH(1)
117         UPLIM = D1MACH(2) / 5.0d0
118         XLU225 = (LOLIM*UPLIM)**2 / 25.0d0
119         Y2236L = -2.236d0/SQRT(LOLIM)
120         ERRTOL = (D1MACH(4)/16.0d0) ** (1.0d0/6.0d0)
121      END IF
122      YN = ABS(Y)
123      IF (X .GT. UPLIM .OR. YN .GT. UPLIM) THEN
124         RC = 0.0d0
125       CALL DERM1 ('DRCVAL',3,0,'X > UPLIM or ABS(Y) > UPLIM','X',X,',')
126         CALL DERV1 ('Y',Y,',')
127         CALL DERV1 ('UPLIM',UPLIM,'.')
128         IERR = 3
129         RETURN
130      END IF
131      IF (X + YN .LT. LOLIM) THEN
132         RC = 0.0d0
133         CALL DERM1 ('DRCVAL',2,0,'X + ABS(Y) < LOLIM','X',X,',')
134         CALL DERV1 ('Y',Y,',')
135         CALL DERV1 ('LOLIM',LOLIM,'.')
136         IERR = 2
137         RETURN
138      END IF
139      IF (X + YN .GT. UPLIM) THEN
140         RC = 0.0d0
141         CALL DERM1 ('DRCVAL',3,0,'X + ABS(Y) > UPLIM','X',X,',')
142         CALL DERV1 ('Y',Y,',')
143         CALL DERV1 ('UPLIM',UPLIM,'.')
144         IERR = 3
145         RETURN
146      END IF
147      IF (Y.LT.Y2236L .AND. X.GT.0.0d0 .AND. X.LT.XLU225) THEN
148         RC = 0.0d0
149         CALL DERM1 ('DRCVAL',4,0,
150     1      'Y < -2.236/SQRT(LOLIM) AND 0 < X < (LOLIM*UPLIM)**2/25',
151     2      'X',X,',')
152         CALL DERV1 ('Y',Y,',')
153         CALL DERV1 ('LOLIM',LOLIM,',')
154         CALL DERV1 ('UPLIM',UPLIM,'.')
155         IERR = 4
156         RETURN
157      END IF
158C
159      IERR = 0
160      IF (Y.GT.0.0d0) THEN
161         XN = X
162         W = 1.0d0
163      ELSE
164C        TRANSFORM TO POSITIVE Y
165         XN = X - Y
166         W = SQRT(X) / SQRT(XN)
167      END IF
168C
169   20 CONTINUE
170         MU = (XN + YN + YN) / 3.0d0
171         SN = (YN + MU) / MU - 2.0d0
172         IF (ABS(SN) .LT. ERRTOL) THEN
173            RC = W * (1.0d0+SN*SN*(C1+SN*(C2+SN*(C3+SN*C4)))) / SQRT(MU)
174            RETURN
175         END IF
176         LAMDA = 2.0 * SQRT(XN) * SQRT(YN) + YN
177         XN = (XN + LAMDA) * 0.25d0
178         YN = (YN + LAMDA) * 0.25d0
179      GO TO 20
180      END
181