]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MINICERN/mathlib/gen/c/besi064.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing for Sun
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / besi064.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1996/04/01 15:02:00 mclareni
6* Mathlib gen
7*
8*
9#include "gen/pilot.h"
10#if !defined(CERNLIB_DOUBLE)
11 FUNCTION BESI0(X)
12#endif
13#if defined(CERNLIB_DOUBLE)
14 FUNCTION DBESI0(X)
15#include "gen/imp64.inc"
16#endif
17 LOGICAL LEX
18 CHARACTER NAME0*(*),NAME1*(*),NAME0E*(*),NAME1E*(*)
19 CHARACTER*80 ERRTXT
20 DIMENSION CI(0:24,0:1),CK(0:16,0:1)
21
22#if !defined(CERNLIB_DOUBLE)
23 PARAMETER (NAME0 = 'BESK0', NAME0E = 'EBESK0')
24 PARAMETER (NAME1 = 'BESK1', NAME1E = 'EBESK1')
25#endif
26#if defined(CERNLIB_DOUBLE)
27 PARAMETER (NAME0 = 'BESK0/DBESK0', NAME0E = 'EBESK0/DBESK0')
28 PARAMETER (NAME1 = 'BESK1/DBESK1', NAME1E = 'EBESK1/DBESK1')
29#endif
30 PARAMETER (EPS=1D-15)
31 PARAMETER (Z1 = 1, HF = Z1/2)
32 PARAMETER (PI = 3.14159 26535 89793 24D0)
33 PARAMETER (CE = 0.57721 56649 01532 86D0)
34 PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
35
36 DATA CI( 0,0) /+1.00827 92054 58740 032D0/
37 DATA CI( 1,0) /+0.00844 51226 24920 943D0/
38 DATA CI( 2,0) /+0.00017 27006 30777 567D0/
39 DATA CI( 3,0) /+0.00000 72475 91099 959D0/
40 DATA CI( 4,0) /+0.00000 05135 87726 878D0/
41 DATA CI( 5,0) /+0.00000 00568 16965 808D0/
42 DATA CI( 6,0) /+0.00000 00085 13091 223D0/
43 DATA CI( 7,0) /+0.00000 00012 38425 364D0/
44 DATA CI( 8,0) /+0.00000 00000 29801 672D0/
45 DATA CI( 9,0) /-0.00000 00000 78956 698D0/
46 DATA CI(10,0) /-0.00000 00000 33127 128D0/
47 DATA CI(11,0) /-0.00000 00000 04497 339D0/
48 DATA CI(12,0) /+0.00000 00000 01799 790D0/
49 DATA CI(13,0) /+0.00000 00000 00965 748D0/
50 DATA CI(14,0) /+0.00000 00000 00038 604D0/
51 DATA CI(15,0) /-0.00000 00000 00104 039D0/
52 DATA CI(16,0) /-0.00000 00000 00023 950D0/
53 DATA CI(17,0) /+0.00000 00000 00009 554D0/
54 DATA CI(18,0) /+0.00000 00000 00004 443D0/
55 DATA CI(19,0) /-0.00000 00000 00000 859D0/
56 DATA CI(20,0) /-0.00000 00000 00000 709D0/
57 DATA CI(21,0) /+0.00000 00000 00000 087D0/
58 DATA CI(22,0) /+0.00000 00000 00000 112D0/
59 DATA CI(23,0) /-0.00000 00000 00000 012D0/
60 DATA CI(24,0) /-0.00000 00000 00000 018D0/
61
62 DATA CI( 0,1) /+0.97580 06023 26285 926D0/
63 DATA CI( 1,1) /-0.02446 74429 63276 385D0/
64 DATA CI( 2,1) /-0.00027 72053 60763 829D0/
65 DATA CI( 3,1) /-0.00000 97321 46728 020D0/
66 DATA CI( 4,1) /-0.00000 06297 24238 640D0/
67 DATA CI( 5,1) /-0.00000 00659 61142 154D0/
68 DATA CI( 6,1) /-0.00000 00096 13872 919D0/
69 DATA CI( 7,1) /-0.00000 00014 01140 901D0/
70 DATA CI( 8,1) /-0.00000 00000 47563 167D0/
71 DATA CI( 9,1) /+0.00000 00000 81530 681D0/
72 DATA CI(10,1) /+0.00000 00000 35408 148D0/
73 DATA CI(11,1) /+0.00000 00000 05102 564D0/
74 DATA CI(12,1) /-0.00000 00000 01804 409D0/
75 DATA CI(13,1) /-0.00000 00000 01023 594D0/
76 DATA CI(14,1) /-0.00000 00000 00052 678D0/
77 DATA CI(15,1) /+0.00000 00000 00107 094D0/
78 DATA CI(16,1) /+0.00000 00000 00026 120D0/
79 DATA CI(17,1) /-0.00000 00000 00009 561D0/
80 DATA CI(18,1) /-0.00000 00000 00004 713D0/
81 DATA CI(19,1) /+0.00000 00000 00000 829D0/
82 DATA CI(20,1) /+0.00000 00000 00000 743D0/
83 DATA CI(21,1) /-0.00000 00000 00000 080D0/
84 DATA CI(22,1) /-0.00000 00000 00000 117D0/
85 DATA CI(23,1) /+0.00000 00000 00000 011D0/
86 DATA CI(24,1) /+0.00000 00000 00000 019D0/
87
88 DATA CK( 0,0) /+0.98840 81742 30825 800D0/
89 DATA CK( 1,0) /-0.01131 05046 46928 281D0/
90 DATA CK( 2,0) /+0.00026 95326 12762 724D0/
91 DATA CK( 3,0) /-0.00001 11066 85196 665D0/
92 DATA CK( 4,0) /+0.00000 06325 75108 500D0/
93 DATA CK( 5,0) /-0.00000 00450 47337 641D0/
94 DATA CK( 6,0) /+0.00000 00037 92996 456D0/
95 DATA CK( 7,0) /-0.00000 00003 64547 179D0/
96 DATA CK( 8,0) /+0.00000 00000 39043 756D0/
97 DATA CK( 9,0) /-0.00000 00000 04579 936D0/
98 DATA CK(10,0) /+0.00000 00000 00580 811D0/
99 DATA CK(11,0) /-0.00000 00000 00078 832D0/
100 DATA CK(12,0) /+0.00000 00000 00011 360D0/
101 DATA CK(13,0) /-0.00000 00000 00001 727D0/
102 DATA CK(14,0) /+0.00000 00000 00000 275D0/
103 DATA CK(15,0) /-0.00000 00000 00000 046D0/
104 DATA CK(16,0) /+0.00000 00000 00000 008D0/
105
106 DATA CK( 0,1) /+1.03595 08587 72358 331D0/
107 DATA CK( 1,1) /+0.03546 52912 43331 114D0/
108 DATA CK( 2,1) /-0.00046 84750 28166 889D0/
109 DATA CK( 3,1) /+0.00001 61850 63810 053D0/
110 DATA CK( 4,1) /-0.00000 08451 72048 124D0/
111 DATA CK( 5,1) /+0.00000 00571 32218 103D0/
112 DATA CK( 6,1) /-0.00000 00046 45554 607D0/
113 DATA CK( 7,1) /+0.00000 00004 35417 339D0/
114 DATA CK( 8,1) /-0.00000 00000 45757 297D0/
115 DATA CK( 9,1) /+0.00000 00000 05288 133D0/
116 DATA CK(10,1) /-0.00000 00000 00662 613D0/
117 DATA CK(11,1) /+0.00000 00000 00089 048D0/
118 DATA CK(12,1) /-0.00000 00000 00012 726D0/
119 DATA CK(13,1) /+0.00000 00000 00001 921D0/
120 DATA CK(14,1) /-0.00000 00000 00000 305D0/
121 DATA CK(15,1) /+0.00000 00000 00000 050D0/
122 DATA CK(16,1) /-0.00000 00000 00000 009D0/
123
124 NU=0
125 LEX=.FALSE.
126 GO TO 6
127
128#if !defined(CERNLIB_DOUBLE)
129 ENTRY EBESI0(X)
130#endif
131#if defined(CERNLIB_DOUBLE)
132 ENTRY DEBSI0(X)
133#endif
134 NU=0
135 LEX=.TRUE.
136 GO TO 6
137
138#if !defined(CERNLIB_DOUBLE)
139 ENTRY BESI1(X)
140#endif
141#if defined(CERNLIB_DOUBLE)
142 ENTRY DBESI1(X)
143#endif
144 NU=1
145 LEX=.FALSE.
146 GO TO 6
147
148#if !defined(CERNLIB_DOUBLE)
149 ENTRY EBESI1(X)
150#endif
151#if defined(CERNLIB_DOUBLE)
152 ENTRY DEBSI1(X)
153#endif
154 NU=1
155 LEX=.TRUE.
156
157 6 V=ABS(X)
158 IF(V .LT. 8) THEN
159 Y=(HF*V)**2
160 XL=NU+2
161 A0=1
162 A1=1+2*Y/((XL+1)*(XL-1))
163 A2=1+Y*(4+3*Y/((XL+2)*XL))/((XL+3)*(XL-1))
164 B0=1
165 B1=1-Y/(XL+1)
166 B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
167 W1=3+XL
168 V1=3-XL
169 V3=XL-1
170 V2=V3+V3
171 C=0
172 DO 3 N = 3,30
173 C0=C
174 FN=N
175 W1=W1+2
176 W2=W1-1
177 W3=W2-1
178 W4=W3-1
179 W5=W4-1
180 W6=W5-1
181 V1=V1+1
182 V2=V2+1
183 V3=V3+1
184 U1=FN*W4
185 E=V3/(U1*W3)
186 U2=E*Y
187 F1=1+Y*V1/(U1*W1)
188 F2=(1+Y*V2/(V3*W2*W5))*U2
189 F3=-Y*Y*U2/(W4*W5*W5*W6)
190 A=F1*A2+F2*A1+F3*A0
191 B=F1*B2+F2*B1+F3*B0
192 C=A/B
193 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 4
194 A0=A1
195 A1=A2
196 A2=A
197 B0=B1
198 B1=B2
199 B2=B
200 3 CONTINUE
201 4 H=C
202 IF(NU .EQ. 1) H=HF*X*H
203 IF(LEX) H=EXP(-V)*H
204 ELSE
205 R=1/V
206 H=16*R-1
207 ALFA=H+H
208 B1=0
209 B2=0
210 DO 1 I = 24,0,-1
211 B0=CI(I,NU)+ALFA*B1-B2
212 B2=B1
213 1 B1=B0
214 H=SQRT(RPI2*R)*(B0-H*B2)
215 IF(NU*X .LT. 0) H=-H
216 IF(.NOT.LEX) H=EXP(V)*H
217 ENDIF
218 GO TO 9
219
220#if !defined(CERNLIB_DOUBLE)
221 ENTRY BESK0(X)
222#endif
223#if defined(CERNLIB_DOUBLE)
224 ENTRY DBESK0(X)
225#endif
226 NU=0
227 LEX=.FALSE.
228 GO TO 8
229
230#if !defined(CERNLIB_DOUBLE)
231 ENTRY EBESK0(X)
232#endif
233#if defined(CERNLIB_DOUBLE)
234 ENTRY DEBSK0(X)
235#endif
236 NU=0
237 LEX=.TRUE.
238 GO TO 8
239
240#if !defined(CERNLIB_DOUBLE)
241 ENTRY BESK1(X)
242#endif
243#if defined(CERNLIB_DOUBLE)
244 ENTRY DBESK1(X)
245#endif
246 NU=1
247 LEX=.FALSE.
248 GO TO 8
249
250#if !defined(CERNLIB_DOUBLE)
251 ENTRY EBESK1(X)
252#endif
253#if defined(CERNLIB_DOUBLE)
254 ENTRY DEBSK1(X)
255#endif
256 NU=1
257 LEX=.TRUE.
258
259 8 IF(X .LE. 0) THEN
260 H=0
261 WRITE(ERRTXT,101) X
262 IF(NU .EQ. 0 .AND. .NOT.LEX) CALL MTLPRT(NAME0 ,'C313.1',ERRTXT)
263 IF(NU .EQ. 0 .AND. LEX) CALL MTLPRT(NAME0E,'C313.1',ERRTXT)
264 IF(NU .EQ. 1 .AND. .NOT.LEX) CALL MTLPRT(NAME1 ,'C313.1',ERRTXT)
265 IF(NU .EQ. 1 .AND. LEX) CALL MTLPRT(NAME1E,'C313.1',ERRTXT)
266 ELSEIF(X .LT. 1) THEN
267 B=HF*X
268 BK=-(LOG(B)+CE)
269 F=BK
270 P=HF
271 Q=HF
272 C=1
273 D=B**2
274 BK1=P
275 DO 11 N = 1,15
276 FN=N
277 RFN=1/FN
278 P=P*RFN
279 Q=Q*RFN
280 F=(F+P+Q)*RFN
281 C=C*D*RFN
282 G=C*(P-FN*F)
283 H=C*F
284 BK=BK+H
285 BK1=BK1+G
286 IF(BK1*H+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
287 11 CONTINUE
288 12 H=BK
289 IF(NU .EQ. 1) H=BK1/B
290 IF(LEX) H=EXP(X)*H
291 ELSEIF(X .LE. 5) THEN
292 XN=4*NU**2
293 A=9-XN
294 B=25-XN
295 C=768*X**2
296 C0=48*X
297 A0=1
298 A1=(16*X+7+XN)/A
299 A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
300 B0=1
301 B1=(16*X+9-XN)/A
302 B2=(C+C0*B)/(A*B)+1
303 C=0
304 DO 24 N = 3,30
305 C0=C
306 FN=N
307 FN2=FN+FN
308 FN1=FN2-1
309 FN3=FN1/(FN2-3)
310 FN4=12*FN**2-(1-XN)
311 FN5=16*FN1*X
312 RAN=1/((FN2+1)**2-XN)
313 F1=FN3*(FN4-20*FN)+FN5
314 F2=28*FN-FN4-8+FN5
315 F3=FN3*((FN2-5)**2-XN)
316 A=(F1*A2+F2*A1+F3*A0)*RAN
317 B=(F1*B2+F2*B1+F3*B0)*RAN
318 C=A/B
319 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
320 A0=A1
321 A1=A2
322 A2=A
323 B0=B1
324 B1=B2
325 B2=B
326 24 CONTINUE
327 25 H=C/SQRT(RPIH*X)
328 IF(.NOT.LEX) H=EXP(-X)*H
329 ELSE
330 R=1/X
331 H=10*R-1
332 ALFA=H+H
333 B1=0
334 B2=0
335 DO 23 I = 16,0,-1
336 B0=CK(I,NU)+ALFA*B1-B2
337 B2=B1
338 23 B1=B0
339 H=SQRT(PIH*R)*(B0-H*B2)
340 IF(.NOT.LEX) H=EXP(-X)*H
341 ENDIF
342 9 CONTINUE
343#if !defined(CERNLIB_DOUBLE)
344 BESI0=H
345#endif
346#if defined(CERNLIB_DOUBLE)
347 DBESI0=H
348#endif
349 RETURN
350 101 FORMAT(' NON-POSITIVE ARGUMENT X = ',1P,E15.6)
351 END