]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MINICERN/mathlib/gen/c/bsir464.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / bsir464.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1996/04/01 15:02:03 mclareni
6* Mathlib gen
7*
8*
9#include "gen/pilot.h"
10#if !defined(CERNLIB_DOUBLE)
11 FUNCTION BSIR4(X,NU)
12#endif
13#if defined(CERNLIB_DOUBLE)
14 FUNCTION DBSIR4(X,NU)
15#include "gen/imp64.inc"
16#endif
17 CHARACTER*80 ERRTXT
18 CHARACTER NAMEI*(*),NAMEK*(*),NAMEIE*(*),NAMEKE*(*)
19#if !defined(CERNLIB_DOUBLE)
20 PARAMETER (NAMEI = 'BSIR4/DBSIR4', NAMEIE = 'EBSIR4/DEBIR4')
21 PARAMETER (NAMEK = 'BSKR4/DBSKR4', NAMEKE = 'EBSKR4/DEBKR4')
22#endif
23#if defined(CERNLIB_DOUBLE)
24 PARAMETER (NAMEI = 'BSIR4/DBSIR4', NAMEIE = 'EBSIR4/DEBIR4')
25 PARAMETER (NAMEK = 'BSKR4/DBSKR4', NAMEKE = 'EBSKR4/DEBKR4')
26#endif
27 LOGICAL LEX
28
29 DIMENSION BC(0:23,2),CC(0:15,2),PP(-3:3),GG(-3:3)
30
31 PARAMETER (EPS = 1D-14)
32 PARAMETER (Z1 = 1, HF =Z1/2)
33 PARAMETER (PI = 3.14159 26535 89793 24D0)
34 PARAMETER (W2 = 1.41421 35623 73095 05D0)
35 PARAMETER (G1 = 3.62560 99082 21908 31D0)
36 PARAMETER (G3 = 1.22541 67024 65177 65D0)
37 PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI = 1/PI, RW2 = 1/W2)
38 PARAMETER (C1 = PI/(2*W2))
39 PARAMETER (GM = 2*(1/G3-4/G1), GP = (4/G1+1/G3)/2)
40
41 DATA GG(-3) /0.27581 56628 30209 31D0/, PP(-3) /-0.75D0/
42 DATA GG(-1) /0.81604 89390 98262 98D0/, PP(-1) /-0.25D0/
43 DATA GG( 1) /1.10326 26513 20837 26D0/, PP( 1) / 0.25D0/
44 DATA GG( 3) /1.08806 52521 31017 31D0/, PP( 3) / 0.75D0/
45
46 DATA BC( 0,1) / 1.00619 92270 14122 57D0/
47 DATA BC( 1,1) / 0.00631 99620 31140 72D0/
48 DATA BC( 2,1) / 0.00012 56131 27965 64D0/
49 DATA BC( 3,1) / 0.00000 52052 40761 57D0/
50 DATA BC( 4,1) / 0.00000 03591 84411 39D0/
51 DATA BC( 5,1) / 0.00000 00355 85362 89D0/
52 DATA BC( 6,1) / 0.00000 00036 05011 66D0/
53 DATA BC( 7,1) /-0.00000 00001 26294 10D0/
54 DATA BC( 8,1) /-0.00000 00002 96595 12D0/
55 DATA BC( 9,1) /-0.00000 00001 18337 70D0/
56 DATA BC(10,1) /-0.00000 00000 21655 68D0/
57 DATA BC(11,1) / 0.00000 00000 03032 04D0/
58 DATA BC(12,1) / 0.00000 00000 03041 10D0/
59 DATA BC(13,1) / 0.00000 00000 00530 77D0/
60 DATA BC(14,1) /-0.00000 00000 00204 53D0/
61 DATA BC(15,1) /-0.00000 00000 00105 49D0/
62 DATA BC(16,1) / 0.00000 00000 00005 50D0/
63 DATA BC(17,1) / 0.00000 00000 00014 36D0/
64 DATA BC(18,1) / 0.00000 00000 00001 14D0/
65 DATA BC(19,1) /-0.00000 00000 00001 87D0/
66 DATA BC(20,1) /-0.00000 00000 00000 32D0/
67 DATA BC(21,1) / 0.00000 00000 00000 26D0/
68 DATA BC(22,1) / 0.00000 00000 00000 06D0/
69 DATA BC(23,1) /-0.00000 00000 00000 04D0/
70
71 DATA BC( 0,2) / 0.98980 19115 24008 91D0/
72 DATA BC( 1,2) /-0.01035 09365 14827 02D0/
73 DATA BC( 2,2) /-0.00015 85263 84973 08D0/
74 DATA BC( 3,2) /-0.00000 60527 21962 69D0/
75 DATA BC( 4,2) /-0.00000 04158 38597 31D0/
76 DATA BC( 5,2) /-0.00000 00487 99346 57D0/
77 DATA BC( 6,2) /-0.00000 00089 86835 44D0/
78 DATA BC( 7,2) /-0.00000 00019 83283 58D0/
79 DATA BC( 8,2) /-0.00000 00003 58969 60D0/
80 DATA BC( 9,2) /-0.00000 00000 08766 62D0/
81 DATA BC(10,2) / 0.00000 00000 25819 45D0/
82 DATA BC(11,2) / 0.00000 00000 09780 24D0/
83 DATA BC(12,2) / 0.00000 00000 00565 05D0/
84 DATA BC(13,2) /-0.00000 00000 00851 66D0/
85 DATA BC(14,2) /-0.00000 00000 00270 25D0/
86 DATA BC(15,2) / 0.00000 00000 00040 96D0/
87 DATA BC(16,2) / 0.00000 00000 00040 50D0/
88 DATA BC(17,2) / 0.00000 00000 00001 11D0/
89 DATA BC(18,2) /-0.00000 00000 00005 25D0/
90 DATA BC(19,2) /-0.00000 00000 00000 70D0/
91 DATA BC(20,2) / 0.00000 00000 00000 70D0/
92 DATA BC(21,2) / 0.00000 00000 00000 14D0/
93 DATA BC(22,2) /-0.00000 00000 00000 10D0/
94 DATA BC(23,2) /-0.00000 00000 00000 02D0/
95
96 DATA CC( 0,1) / 0.99128 81656 75147 07D0/
97 DATA CC( 1,1) /-0.00850 62567 20022 24D0/
98 DATA CC( 2,1) / 0.00019 70491 57408 35D0/
99 DATA CC( 3,1) /-0.00000 80377 10166 54D0/
100 DATA CC( 4,1) / 0.00000 04554 01498 43D0/
101 DATA CC( 5,1) /-0.00000 00323 27352 82D0/
102 DATA CC( 6,1) / 0.00000 00027 16130 28D0/
103 DATA CC( 7,1) /-0.00000 00002 60644 07D0/
104 DATA CC( 8,1) / 0.00000 00000 27882 69D0/
105 DATA CC( 9,1) /-0.00000 00000 03267 69D0/
106 DATA CC(10,1) / 0.00000 00000 00414 09D0/
107 DATA CC(11,1) /-0.00000 00000 00056 17D0/
108 DATA CC(12,1) / 0.00000 00000 00008 09D0/
109 DATA CC(13,1) /-0.00000 00000 00001 23D0/
110 DATA CC(14,1) / 0.00000 00000 00000 20D0/
111 DATA CC(15,1) /-0.00000 00000 00000 03D0/
112
113 DATA CC( 0,2) / 1.01476 24350 64637 87D0/
114 DATA CC( 1,2) / 0.01449 34617 87809 66D0/
115 DATA CC( 2,2) /-0.00025 87162 07241 80D0/
116 DATA CC( 3,2) / 0.00000 96912 18911 49D0/
117 DATA CC( 4,2) /-0.00000 05261 29313 99D0/
118 DATA CC( 5,2) / 0.00000 00363 96854 29D0/
119 DATA CC( 6,2) /-0.00000 00030 05472 76D0/
120 DATA CC( 7,2) / 0.00000 00002 84827 80D0/
121 DATA CC( 8,2) /-0.00000 00000 30182 91D0/
122 DATA CC( 9,2) / 0.00000 00000 03511 10D0/
123 DATA CC(10,2) /-0.00000 00000 00442 27D0/
124 DATA CC(11,2) / 0.00000 00000 00059 70D0/
125 DATA CC(12,2) /-0.00000 00000 00008 56D0/
126 DATA CC(13,2) / 0.00000 00000 00001 30D0/
127 DATA CC(14,2) /-0.00000 00000 00000 21D0/
128 DATA CC(15,2) / 0.00000 00000 00000 03D0/
129
130 LEX=.FALSE.
131 GO TO 8
132
133#if !defined(CERNLIB_DOUBLE)
134 ENTRY EBSIR4(X,NU)
135#endif
136#if defined(CERNLIB_DOUBLE)
137 ENTRY DEBIR4(X,NU)
138#endif
139 LEX=.TRUE.
140
141 8 MU=ABS(NU)
142 IF(MU .NE. 1 .AND. MU .NE. 2 .AND. MU .NE. 3 .OR.
143 1 NU .LT. 0 .AND. X .LE. 0 .OR. NU .GT. 0 .AND. X .LT. 0) THEN
144 S=0
145 WRITE(ERRTXT,101) X,NU
146 IF(.NOT.LEX) CALL MTLPRT(NAMEI ,'C327.1',ERRTXT)
147 IF( LEX) CALL MTLPRT(NAMEIE,'C327.1',ERRTXT)
148 ELSEIF(X .EQ. 0) THEN
149 S=0
150 ELSEIF(NU .EQ. -2) THEN
151 IF(LEX) THEN
152 S=HF*(1+EXP(-X-X))/SQRT(PIH*X)
153 ELSE
154 S=COSH(X)/SQRT(PIH*X)
155 ENDIF
156 ELSEIF(NU .EQ. 2) THEN
157 IF(LEX) THEN
158 IF(X .LT. HF) THEN
159 S=SINH(X)*EXP(-X)/SQRT(PIH*X)
160 ELSE
161 S=HF*(1-EXP(-X-X))/SQRT(PIH*X)
162 ENDIF
163 ELSE
164 S=SINH(X)/SQRT(PIH*X)
165 ENDIF
166 ELSEIF(X .LT. 8) THEN
167 Y=(HF*X)**2
168 XN=PP(NU)
169 XL=XN+2
170 A0=1
171 A1=1+2*Y/((XL+1)*(XN+1))
172 A2=1+Y*(4+3*Y/((XL+2)*(XN+2)))/((XL+3)*(XN+1))
173 B0=1
174 B1=1-Y/(XL+1)
175 B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
176 T1=3+XL
177 V1=3-XL
178 V3=XL-1
179 V2=V3+V3
180 C=0
181 DO 33 N = 3,30
182 C0=C
183 T1=T1+2
184 T2=T1-1
185 T3=T2-1
186 T4=T3-1
187 T5=T4-1
188 T6=T5-1
189 V1=V1+1
190 V2=V2+1
191 V3=V3+1
192 U1=N*T4
193 E=V3/(U1*T3)
194 U2=E*Y
195 F1=1+Y*V1/(U1*T1)
196 F2=(1+Y*V2/(V3*T2*T5))*U2
197 F3=-Y*Y*U2/(T4*T5*T5*T6)
198 A=F1*A2+F2*A1+F3*A0
199 B=F1*B2+F2*B1+F3*B0
200 C=A/B
201 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 34
202 A0=A1
203 A1=A2
204 A2=A
205 B0=B1
206 B1=B2
207 B2=B
208 33 CONTINUE
209 34 S=GG(NU)*(HF*X)**PP(NU)*C
210 IF(LEX) S=EXP(-X)*S
211 ELSE
212 K=(MU+1)/2
213 R=1/X
214 W=SQRT(RPI*R)
215 H=16*R-1
216 ALFA=H+H
217 B1=0
218 B2=0
219 DO 2 I = 23,0,-1
220 B0=BC(I,K)+ALFA*B1-B2
221 B2=B1
222 2 B1=B0
223 S=RW2*W*(B0-H*B2)
224 IF(.NOT.LEX) S=EXP(X)*S
225 T=0
226 IF(NU .LT. 0) THEN
227 H=10*R-1
228 ALFA=H+H
229 B1=0
230 B2=0
231 DO 3 I = 15,0,-1
232 B0=CC(I,K)+ALFA*B1-B2
233 B2=B1
234 3 B1=B0
235 R=EXP(-X)
236 T=W*R*(B0-H*B2)
237 IF(LEX) T=R*T
238 ENDIF
239 S=S+T
240 ENDIF
241 GO TO 99
242
243#if !defined(CERNLIB_DOUBLE)
244 ENTRY BSKR4(X,NU)
245#endif
246#if defined(CERNLIB_DOUBLE)
247 ENTRY DBSKR4(X,NU)
248#endif
249 LEX=.FALSE.
250 GO TO 9
251
252#if !defined(CERNLIB_DOUBLE)
253 ENTRY EBSKR4(X,NU)
254#endif
255#if defined(CERNLIB_DOUBLE)
256 ENTRY DEBKR4(X,NU)
257#endif
258 LEX=.TRUE.
259
260 9 MU=ABS(NU)
261 IF(MU .NE. 1 .AND. MU .NE. 2 .AND. MU .NE. 3 .OR. X .LE. 0) THEN
262 S=0
263 WRITE(ERRTXT,101) X,NU
264 IF(.NOT.LEX) CALL MTLPRT(NAMEK ,'C327.1',ERRTXT)
265 IF( LEX) CALL MTLPRT(NAMEKE,'C327.1',ERRTXT)
266 ELSEIF(MU .EQ. 2) THEN
267 S=SQRT(PIH/X)
268 IF(.NOT.LEX) S=EXP(-X)*S
269 ELSEIF(X .LE. 1) THEN
270 A0=PP(-1)
271 B=HF*X
272 D=-LOG(B)
273 F=A0*D
274 E=EXP(F)
275 G=(GM*A0+GP)*E
276 BK=C1*(HF*GM*(E+1/E)+GP*D*SINH(F)/F)
277 F=BK
278 E=A0**2
279 P=HF*C1*G
280 Q=HF/G
281 C=1
282 D=B**2
283 BK1=P
284 DO 11 N = 1,15
285 FN=N
286 F=(FN*F+P+Q)/(FN**2-E)
287 C=C*D/FN
288 P=P/(FN-A0)
289 Q=Q/(FN+A0)
290 G=C*(P-FN*F)
291 H=C*F
292 BK=BK+H
293 BK1=BK1+G
294 IF(H*BK1+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
295 11 CONTINUE
296 12 S=BK
297 IF(MU .EQ. 3) S=BK1/B
298 IF(LEX) S=EXP(X)*S
299 ELSEIF(X .LE. 5) THEN
300 XN=4*PP(MU)**2
301 A=9-XN
302 B=25-XN
303 C=768*X**2
304 C0=48*X
305 A0=1
306 A1=(16*X+7+XN)/A
307 A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
308 B0=1
309 B1=(16*X+9-XN)/A
310 B2=(C+C0*B)/(A*B)+1
311 C=0
312 DO 24 N = 3,30
313 C0=C
314 FN=N
315 FN2=FN+FN
316 FN1=FN2-1
317 FN3=FN1/(FN2-3)
318 FN4=12*FN**2-(1-XN)
319 FN5=16*FN1*X
320 RAN=1/((FN2+1)**2-XN)
321 F1=FN3*(FN4-20*FN)+FN5
322 F2=28*FN-FN4-8+FN5
323 F3=FN3*((FN2-5)**2-XN)
324 A=(F1*A2+F2*A1+F3*A0)*RAN
325 B=(F1*B2+F2*B1+F3*B0)*RAN
326 C=A/B
327 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
328 A0=A1
329 A1=A2
330 A2=A
331 B0=B1
332 B1=B2
333 B2=B
334 24 CONTINUE
335 25 S=C/SQRT(RPIH*X)
336 IF(.NOT.LEX) S=EXP(-X)*S
337 ELSE
338 K=(MU+1)/2
339 R=1/X
340 Y=5*R
341 H=10*R-1
342 ALFA=H+H
343 B1=0
344 B2=0
345 DO 13 I = 15,0,-1
346 B0=CC(I,K)+ALFA*B1-B2
347 B2=B1
348 13 B1=B0
349 S=SQRT(PIH*R)*(B0-H*B2)
350 IF(.NOT.LEX) S=EXP(-X)*S
351 ENDIF
352#if !defined(CERNLIB_DOUBLE)
353 99 BSIR4=S
354#endif
355#if defined(CERNLIB_DOUBLE)
356 99 DBSIR4=S
357#endif
358 RETURN
359 101 FORMAT('INCORRECT ARGUMENT OR INDEX, X = ',1P,E15.6,5X,'NU = ',I5)
360 END