]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:11 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | #if defined(CERNLIB_DOUBLE) | |
11 | FUNCTION DELI3(X,AKP,P) | |
12 | C | |
13 | #include "gen/imp64.inc" | |
14 | C | |
15 | CHARACTER*(*) NAME | |
16 | PARAMETER(NAME='RELI3/DELI3') | |
17 | #endif | |
18 | #if !defined(CERNLIB_DOUBLE) | |
19 | FUNCTION RELI3(X,AKP,P) | |
20 | C | |
21 | CHARACTER*(*) NAME | |
22 | PARAMETER(NAME='RELI3') | |
23 | #endif | |
24 | C | |
25 | C Translation of Algol procedure el3(x,kc,p) in | |
26 | C R. BULIRSCH Numerical Calculation of Elliptic Integrals and | |
27 | C Elliptic Functions III., Numer. Math. 13 (1969) 305-315 | |
28 | C | |
29 | LOGICAL LBO,LBK | |
30 | ||
31 | PARAMETER (ID = 16, IB = 4) | |
32 | PARAMETER (PI = 3.14159 26535 89793 24D0) | |
33 | PARAMETER (AL2 = 0.69314 71805 59945 31D0) | |
34 | PARAMETER (ALB = IB*AL2, RLB = 1/ALB) | |
35 | PARAMETER (Z1 = 1, Z10 = 10, HF = Z1/2, C1 = Z1/10) | |
36 | PARAMETER (ND = ID-2, CA = Z10**(-ID/2), CB = Z10**(-(ID+2))) | |
37 | PARAMETER (ZD = HF/(ND+1)) | |
38 | ||
39 | CHARACTER*80 ERRTXT | |
40 | ||
41 | DIMENSION RA(2:ND),RB(2:ND),RR(2:ND) | |
42 | ||
43 | DATA (RB(K),RA(K),K=2,ND) | |
44 | 2/2.50000000000000000D-01, 7.50000000000000000D-01, | |
45 | 3 1.66666666666666667D-01, 8.33333333333333333D-01, | |
46 | 4 1.25000000000000000D-01, 8.75000000000000000D-01, | |
47 | 5 1.00000000000000000D-01, 9.00000000000000000D-01, | |
48 | 6 8.33333333333333333D-02, 9.16666666666666667D-01, | |
49 | 7 7.14285714285714286D-02, 9.28571428571428571D-01, | |
50 | 8 6.25000000000000000D-02, 9.37500000000000000D-01, | |
51 | 9 5.55555555555555556D-02, 9.44444444444444445D-01, | |
52 | A 5.00000000000000000D-02, 9.50000000000000000D-01, | |
53 | B 4.54545454545454545D-02, 9.54545454545454545D-01, | |
54 | C 4.16666666666666667D-02, 9.58333333333333333D-01, | |
55 | D 3.84615384615384615D-02, 9.61538461538461538D-01, | |
56 | E 3.57142857142857143D-02, 9.64285714285714286D-01/ | |
57 | ||
58 | IF(X .EQ. 0) THEN | |
59 | H=0 | |
60 | GO TO 9 | |
61 | ENDIF | |
62 | HH=X**2 | |
63 | F=P*HH | |
64 | S=AKP | |
65 | IF(S .EQ. 0) S=CA/(1+ABS(X)) | |
66 | T=S**2 | |
67 | PM=HF*T | |
68 | E=HH*T | |
69 | Z=ABS(F) | |
70 | R=ABS(P) | |
71 | H1=HH+1 | |
72 | IF(E .LT. C1 .AND. Z .LT. C1 .AND. T .LT. 1 .AND. R .LT. 1) THEN | |
73 | S=P+PM | |
74 | DO 1 K = 2,ND | |
75 | RR(K)=S | |
76 | PM=PM*T*RA(K) | |
77 | 1 S=S*P+PM | |
78 | U=S*ZD | |
79 | S=U | |
80 | LBO=.FALSE. | |
81 | DO 2 K = ND,2,-1 | |
82 | U=U+(RR(K)-U)*RB(K) | |
83 | LBO=.NOT.LBO | |
84 | V=U | |
85 | IF(LBO) V=-U | |
86 | 2 S=S*HH+V | |
87 | IF(LBO) S=-S | |
88 | U=HF*(U+1) | |
89 | #if defined(CERNLIB_DOUBLE) | |
90 | H=(U-S*H1)*SQRT(H1)*X+U*DASINH(X) | |
91 | #endif | |
92 | #if !defined(CERNLIB_DOUBLE) | |
93 | H=(U-S*H1)*SQRT(H1)*X+U*ASINH(X) | |
94 | #endif | |
95 | ELSE | |
96 | W=1+F | |
97 | IF(W .EQ. 0) THEN | |
98 | H=0 | |
99 | WRITE(ERRTXT,103) X,P | |
100 | CALL MTLPRT(NAME,'C346.3',ERRTXT) | |
101 | GO TO 9 | |
102 | ENDIF | |
103 | PP=P | |
104 | IF(PP .EQ. 0) PP=CB/HH | |
105 | S=ABS(S) | |
106 | Y=ABS(X) | |
107 | G=PP-1 | |
108 | IF(G .EQ. 0) G=CB | |
109 | F=PP-T | |
110 | IF(F .EQ. 0) F=CB*T | |
111 | AM=1-T | |
112 | AP=1+E | |
113 | R=PP*H1 | |
114 | FA=G/(F*PP) | |
115 | LBO=FA .GT. 0 | |
116 | FA=ABS(FA) | |
117 | PZ=ABS(G*F) | |
118 | DE=SQRT(PZ) | |
119 | Q=SQRT(ABS(PP)) | |
120 | PM=PP-MIN(PM,HF) | |
121 | IF(PM .GE. 0) THEN | |
122 | U=SQRT(R*AP) | |
123 | V=Y*DE | |
124 | IF(G .LT. 0) V=-V | |
125 | D=1/Q | |
126 | C=1 | |
127 | ELSE | |
128 | U=SQRT(H1*AP*PZ) | |
129 | YE=Y*Q | |
130 | V=AM*YE | |
131 | Q=-DE/G | |
132 | D=-AM/DE | |
133 | C=0 | |
134 | PZ=AP-R | |
135 | ENDIF | |
136 | IF(LBO) THEN | |
137 | R=V/U | |
138 | Z=1 | |
139 | K=1 | |
140 | IF(PM .LT. 0) THEN | |
141 | H1=Y*SQRT(H1/(AP*FA)) | |
142 | H1=1/H1-H1 | |
143 | Z=H1-2*R | |
144 | R=2+R*H1 | |
145 | IF(R .EQ. 0) R=CB | |
146 | IF(Z .EQ. 0) Z=H1*CB | |
147 | R=R/Z | |
148 | Z=R | |
149 | W=PZ | |
150 | ENDIF | |
151 | U=U/W | |
152 | V=V/W | |
153 | ELSE | |
154 | T=U+ABS(V) | |
155 | LBK=.TRUE. | |
156 | IF(PP .LT. 0) THEN | |
157 | DE=V/PZ | |
158 | YE=2*U*YE | |
159 | U=T/PZ | |
160 | V=-(F+G*E)/T | |
161 | T=PZ*ABS(W) | |
162 | Z=(HH*R*F-G*AP+YE)/T | |
163 | YE=YE/T | |
164 | ELSE | |
165 | DE=V/W | |
166 | YE=0 | |
167 | U=(E+PP)/T | |
168 | V=T/W | |
169 | Z=1 | |
170 | ENDIF | |
171 | IF(S .GT. 1) THEN | |
172 | H1=U | |
173 | U=V | |
174 | V=H1 | |
175 | ENDIF | |
176 | ENDIF | |
177 | Y=1/Y | |
178 | E=S | |
179 | T=1 | |
180 | N=1 | |
181 | L=0 | |
182 | M=0 | |
183 | 3 Y=Y-E/Y | |
184 | IF(Y .EQ. 0) Y =SQRT(E)*CB | |
185 | F=C | |
186 | C=D/Q+C | |
187 | G=E/Q | |
188 | D=2*(F*G+D) | |
189 | Q=G+Q | |
190 | G=T | |
191 | T=S+T | |
192 | N=2*N | |
193 | M=2*M | |
194 | IF(LBO) THEN | |
195 | IF(Z .LT. 0) M=K+M | |
196 | K=0 | |
197 | IF(R .GT. 0) K=1 | |
198 | IF(R .LT. 0) K=-1 | |
199 | H1=E/(U**2+V**2) | |
200 | U=U*(1+H1) | |
201 | V=V*(1-H1) | |
202 | ELSE | |
203 | R=U/V | |
204 | H1=Z*R | |
205 | Z=H1*Z | |
206 | HH=E/V | |
207 | IF(LBK) THEN | |
208 | DE=DE/U | |
209 | YE=YE*(H1+1/H1)+DE*(1+R) | |
210 | DE=DE*(U-HH) | |
211 | LBK=ABS(YE) .LT. 1 | |
212 | ELSE | |
213 | A=LOG(Z) | |
214 | K=RLB*A | |
215 | IF(Z .GE. 1) K=K+1 | |
216 | Z=EXP(A-K*ALB) | |
217 | M=M+K | |
218 | ENDIF | |
219 | ENDIF | |
220 | IF(ABS(G-S) .GT. CA*G) THEN | |
221 | IF(LBO) THEN | |
222 | G=HF*(1/R-R) | |
223 | HH=U+V*G | |
224 | H1=G*U-V | |
225 | IF(HH .EQ. 0) HH=U*CB | |
226 | IF(H1 .EQ. 0) H1=V*CB | |
227 | Z=R*H1 | |
228 | R=HH/H1 | |
229 | ELSE | |
230 | U=U+E/U | |
231 | V=V+HH | |
232 | ENDIF | |
233 | S=2*SQRT(E) | |
234 | E=S*T | |
235 | L=2*L | |
236 | IF(Y .LT. 0) L=L+1 | |
237 | GO TO 3 | |
238 | ENDIF | |
239 | IF(Y .LT. 0) L=L+1 | |
240 | E=(ATAN(T/Y)+PI*L)*(C*T+D)/(T*(T+Q)) | |
241 | IF(LBO) THEN | |
242 | H1=V/(T+U) | |
243 | Z=1-R*H1 | |
244 | H1=R+H1 | |
245 | IF(Z .EQ. 0) Z=CB | |
246 | IF(Z .LT. 0) THEN | |
247 | IF(H1 .GT. 0) M=M+1 | |
248 | IF(H1 .LT. 0) M=M-1 | |
249 | ENDIF | |
250 | S=ATAN(H1/Z)+M*PI | |
251 | ELSE | |
252 | IF(LBK) THEN | |
253 | #if defined(CERNLIB_DOUBLE) | |
254 | S=HF*DASINH(YE) | |
255 | #endif | |
256 | #if !defined(CERNLIB_DOUBLE) | |
257 | S=HF*ASINH(YE) | |
258 | #endif | |
259 | ELSE | |
260 | S=HF*(LOG(Z)+M*ALB) | |
261 | ENDIF | |
262 | ENDIF | |
263 | H=(E+SQRT(FA)*S)/N | |
264 | IF(X .LT. 0) H=-H | |
265 | ENDIF | |
266 | #if defined(CERNLIB_DOUBLE) | |
267 | 9 DELI3=H | |
268 | #endif | |
269 | #if !defined(CERNLIB_DOUBLE) | |
270 | 9 RELI3=H | |
271 | #endif | |
272 | RETURN | |
273 | 103 FORMAT('FUNCTION SINGULAR FOR ',1P, | |
274 | 1 'X = ',D15.8,3X,'P = ',D15.8,4X,'(P*X**2 = -1)') | |
275 | END |