]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MINICERN/mathlib/gen/c/bsir364.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 / bsir364.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1996/04/01 15:02:07 mclareni
6* Mathlib gen
7*
8*
9#include "gen/pilot.h"
10#if !defined(CERNLIB_DOUBLE)
11 FUNCTION BSIR3(X,NU)
12#endif
13#if defined(CERNLIB_DOUBLE)
14 FUNCTION DBSIR3(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 = 'BSIR3', NAMEIE = 'EBSIR3')
21 PARAMETER (NAMEK = 'BSKR3', NAMEKE = 'EBSKR3')
22#endif
23#if defined(CERNLIB_DOUBLE)
24 PARAMETER (NAMEI = 'BSIR3/DBSIR3', NAMEIE = 'EBSIR3/DEBIR3')
25 PARAMETER (NAMEK = 'BSKR3/DBSKR3', NAMEKE = 'EBSKR3/DEBKR3')
26#endif
27 LOGICAL LEX
28
29 DIMENSION BC(0:23,2),CC(0:15,2),PP(-2:2),GG(-2:2)
30
31 PARAMETER (EPS = 1D-15)
32 PARAMETER (Z1 = 1, HF = Z1/2)
33 PARAMETER (PI = 3.14159 26535 89793 24D0)
34 PARAMETER (W3 = 1.73205 08075 68877 29D0)
35 PARAMETER (G1 = 2.67893 85347 07747 63D0)
36 PARAMETER (G2 = 1.35411 79394 26400 42D0)
37 PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
38 PARAMETER (C1 = 2*PI/(3*W3))
39 PARAMETER (GM = 3*(1/G2-3/G1)/2, GP = (3/G1+1/G2)/2)
40
41 DATA PP(-2) /-0.66666 66666 66666 67D0/
42 DATA PP(-1) /-0.33333 33333 33333 33D0/
43 DATA PP( 1) / 0.33333 33333 33333 33D0/
44 DATA PP( 2) / 0.66666 66666 66666 67D0/
45
46 DATA GG(-2) / 0.37328 21739 07395 23D0/
47 DATA GG(-1) / 0.73848 81116 21648 31D0/
48 DATA GG( 1) / 1.11984 65217 22185 68D0/
49 DATA GG( 2) / 1.10773 21674 32472 47D0/
50
51 DATA BC( 0,1) / 1.00458 61710 93207 35D0/
52 DATA BC( 1,1) / 0.00467 34791 99873 60D0/
53 DATA BC( 2,1) / 0.00009 08034 04815 04D0/
54 DATA BC( 3,1) / 0.00000 37262 16110 59D0/
55 DATA BC( 4,1) / 0.00000 02520 73237 90D0/
56 DATA BC( 5,1) / 0.00000 00227 82110 77D0/
57 DATA BC( 6,1) / 0.00000 00012 91332 28D0/
58 DATA BC( 7,1) /-0.00000 00006 11915 16D0/
59 DATA BC( 8,1) /-0.00000 00003 75616 85D0/
60 DATA BC( 9,1) /-0.00000 00001 16415 46D0/
61 DATA BC(10,1) /-0.00000 00000 14443 25D0/
62 DATA BC(11,1) / 0.00000 00000 05373 69D0/
63 DATA BC(12,1) / 0.00000 00000 03074 27D0/
64 DATA BC(13,1) / 0.00000 00000 00297 66D0/
65 DATA BC(14,1) /-0.00000 00000 00265 20D0/
66 DATA BC(15,1) /-0.00000 00000 00091 37D0/
67 DATA BC(16,1) / 0.00000 00000 00015 52D0/
68 DATA BC(17,1) / 0.00000 00000 00014 12D0/
69 DATA BC(18,1) /-0.00000 00000 00000 23D0/
70 DATA BC(19,1) /-0.00000 00000 00001 98D0/
71 DATA BC(20,1) /-0.00000 00000 00000 13D0/
72 DATA BC(21,1) / 0.00000 00000 00000 29D0/
73 DATA BC(22,1) / 0.00000 00000 00000 03D0/
74 DATA BC(23,1) /-0.00000 00000 00000 05D0/
75
76 DATA BC( 0,2) / 0.99363 49867 16925 14D0/
77 DATA BC( 1,2) /-0.00646 71526 00616 03D0/
78 DATA BC( 2,2) /-0.00010 60188 22351 55D0/
79 DATA BC( 3,2) /-0.00000 41406 57716 24D0/
80 DATA BC( 4,2) /-0.00000 02916 95418 21D0/
81 DATA BC( 5,2) /-0.00000 00365 71574 33D0/
82 DATA BC( 6,2) /-0.00000 00075 81590 37D0/
83 DATA BC( 7,2) /-0.00000 00019 23008 52D0/
84 DATA BC( 8,2) /-0.00000 00004 20438 80D0/
85 DATA BC( 9,2) /-0.00000 00000 39372 04D0/
86 DATA BC(10,2) / 0.00000 00000 19007 44D0/
87 DATA BC(11,2) / 0.00000 00000 10137 64D0/
88 DATA BC(12,2) / 0.00000 00000 01331 30D0/
89 DATA BC(13,2) /-0.00000 00000 00676 92D0/
90 DATA BC(14,2) /-0.00000 00000 00311 72D0/
91 DATA BC(15,2) / 0.00000 00000 00011 87D0/
92 DATA BC(16,2) / 0.00000 00000 00040 21D0/
93 DATA BC(17,2) / 0.00000 00000 00004 78D0/
94 DATA BC(18,2) /-0.00000 00000 00004 74D0/
95 DATA BC(19,2) /-0.00000 00000 00001 16D0/
96 DATA BC(20,2) / 0.00000 00000 00000 59D0/
97 DATA BC(21,2) / 0.00000 00000 00000 21D0/
98 DATA BC(22,2) /-0.00000 00000 00000 08D0/
99 DATA BC(23,2) /-0.00000 00000 00000 03D0/
100
101 DATA CC( 0,1) / 0.99353 64122 76093 39D0/
102 DATA CC( 1,1) /-0.00631 44392 60798 63D0/
103 DATA CC( 2,1) / 0.00014 30095 80961 13D0/
104 DATA CC( 3,1) /-0.00000 57870 60592 03D0/
105 DATA CC( 4,1) / 0.00000 03265 50333 20D0/
106 DATA CC( 5,1) /-0.00000 00231 23231 95D0/
107 DATA CC( 6,1) / 0.00000 00019 39555 14D0/
108 DATA CC( 7,1) /-0.00000 00001 85897 89D0/
109 DATA CC( 8,1) / 0.00000 00000 19868 42D0/
110 DATA CC( 9,1) /-0.00000 00000 02326 79D0/
111 DATA CC(10,1) / 0.00000 00000 00294 68D0/
112 DATA CC(11,1) /-0.00000 00000 00039 95D0/
113 DATA CC(12,1) / 0.00000 00000 00005 75D0/
114 DATA CC(13,1) /-0.00000 00000 00000 87D0/
115 DATA CC(14,1) / 0.00000 00000 00000 14D0/
116 DATA CC(15,1) /-0.00000 00000 00000 02D0/
117
118 DATA CC( 0,2) / 1.00914 95380 72789 40D0/
119 DATA CC( 1,2) / 0.00897 12068 42483 60D0/
120 DATA CC( 2,2) /-0.00017 13895 98261 54D0/
121 DATA CC( 3,2) / 0.00000 65547 92549 82D0/
122 DATA CC( 4,2) /-0.00000 03595 19190 49D0/
123 DATA CC( 5,2) / 0.00000 00250 24412 19D0/
124 DATA CC( 6,2) /-0.00000 00020 74924 13D0/
125 DATA CC( 7,2) / 0.00000 00001 97223 67D0/
126 DATA CC( 8,2) /-0.00000 00000 20946 47D0/
127 DATA CC( 9,2) / 0.00000 00000 02440 93D0/
128 DATA CC(10,2) /-0.00000 00000 00307 91D0/
129 DATA CC(11,2) / 0.00000 00000 00041 61D0/
130 DATA CC(12,2) /-0.00000 00000 00005 97D0/
131 DATA CC(13,2) / 0.00000 00000 00000 91D0/
132 DATA CC(14,2) /-0.00000 00000 00000 14D0/
133 DATA CC(15,2) / 0.00000 00000 00000 02D0/
134
135 LEX=.FALSE.
136 GO TO 8
137
138#if !defined(CERNLIB_DOUBLE)
139 ENTRY EBSIR3(X,NU)
140#endif
141#if defined(CERNLIB_DOUBLE)
142 ENTRY DEBIR3(X,NU)
143#endif
144 LEX=.TRUE.
145
146 8 MU=ABS(NU)
147 IF(MU .NE. 1 .AND. MU .NE. 2 .OR. NU .LT. 0 .AND. X .LE. 0
148 1 .OR. NU .GT. 0 .AND. X .LT. 0) THEN
149 S=0
150 WRITE(ERRTXT,101) X,NU
151 IF(.NOT.LEX) CALL MTLPRT(NAMEI ,'C340.1',ERRTXT)
152 IF( LEX) CALL MTLPRT(NAMEIE,'C340.1',ERRTXT)
153 ELSEIF(X .EQ. 0) THEN
154 S=0
155 ELSEIF(X .LT. 8) THEN
156 Y=(HF*X)**2
157 XN=PP(NU)
158 XL=XN+2
159 A0=1
160 A1=1+2*Y/((XL+1)*(XN+1))
161 A2=1+Y*(4+3*Y/((XL+2)*(XN+2)))/((XL+3)*(XN+1))
162 B0=1
163 B1=1-Y/(XL+1)
164 B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
165 T1=3+XL
166 V1=3-XL
167 V3=XL-1
168 V2=V3+V3
169 C=0
170 DO 33 N = 3,30
171 C0=C
172 T1=T1+2
173 T2=T1-1
174 T3=T2-1
175 T4=T3-1
176 T5=T4-1
177 T6=T5-1
178 V1=V1+1
179 V2=V2+1
180 V3=V3+1
181 U1=N*T4
182 E=V3/(U1*T3)
183 U2=E*Y
184 F1=1+Y*V1/(U1*T1)
185 F2=(1+Y*V2/(V3*T2*T5))*U2
186 F3=-Y*Y*U2/(T4*T5*T5*T6)
187 A=F1*A2+F2*A1+F3*A0
188 B=F1*B2+F2*B1+F3*B0
189 C=A/B
190 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 34
191 A0=A1
192 A1=A2
193 A2=A
194 B0=B1
195 B1=B2
196 B2=B
197 33 CONTINUE
198 34 S=GG(NU)*(HF*X)**PP(NU)*C
199 IF(LEX) S=EXP(-X)*S
200 ELSE
201 R=1/X
202 W=SQRT(RPI2*R)
203 H=16*R-1
204 ALFA=H+H
205 B1=0
206 B2=0
207 DO 2 I = 23,0,-1
208 B0=BC(I,MU)+ALFA*B1-B2
209 B2=B1
210 2 B1=B0
211 S=W*(B0-H*B2)
212 IF(.NOT.LEX) S=EXP(X)*S
213 T=0
214 IF(NU .LT. 0) THEN
215 H=10*R-1
216 ALFA=H+H
217 B1=0
218 B2=0
219 DO 3 I = 15,0,-1
220 B0=CC(I,MU)+ALFA*B1-B2
221 B2=B1
222 3 B1=B0
223 R=EXP(-X)
224 T=W3*W*R*(B0-H*B2)
225 IF(LEX) T=R*T
226 END IF
227 S=S+T
228 END IF
229 GO TO 99
230
231#if !defined(CERNLIB_DOUBLE)
232 ENTRY BSKR3(X,NU)
233#endif
234#if defined(CERNLIB_DOUBLE)
235 ENTRY DBSKR3(X,NU)
236#endif
237 LEX=.FALSE.
238 GO TO 9
239
240#if !defined(CERNLIB_DOUBLE)
241 ENTRY EBSKR3(X,NU)
242#endif
243#if defined(CERNLIB_DOUBLE)
244 ENTRY DEBKR3(X,NU)
245#endif
246 LEX=.TRUE.
247
248 9 MU=ABS(NU)
249 IF(MU .NE. 1 .AND. MU .NE. 2 .OR. X .LE. 0) THEN
250 S=0
251 WRITE(ERRTXT,101) X,NU
252 IF(.NOT.LEX) CALL MTLPRT(NAMEK ,'C340.1',ERRTXT)
253 IF( LEX) CALL MTLPRT(NAMEKE,'C340.1',ERRTXT)
254 ELSEIF(X .LE. 1) THEN
255 A0=PP(-1)
256 B=HF*X
257 D=-LOG(B)
258 F=A0*D
259 E=EXP(F)
260 G=(GM*A0+GP)*E
261 BK=C1*(HF*GM*(E+1/E)+GP*D*SINH(F)/F)
262 F=BK
263 E=A0**2
264 P=HF*C1*G
265 Q=HF/G
266 C=1
267 D=B**2
268 BK1=P
269 DO 11 N = 1,15
270 FN=N
271 F=(FN*F+P+Q)/(FN**2-E)
272 C=C*D/FN
273 P=P/(FN-A0)
274 Q=Q/(FN+A0)
275 G=C*(P-FN*F)
276 H=C*F
277 BK=BK+H
278 BK1=BK1+G
279 IF(H*BK1+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
280 11 CONTINUE
281 12 S=BK
282 IF(MU .EQ. 2) S=BK1/B
283 IF(LEX) S=EXP(X)*S
284 ELSEIF(X .LE. 5) THEN
285 XN=4*PP(MU)**2
286 A=9-XN
287 B=25-XN
288 C=768*X**2
289 C0=48*X
290 A0=1
291 A1=(16*X+7+XN)/A
292 A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
293 B0=1
294 B1=(16*X+9-XN)/A
295 B2=(C+C0*B)/(A*B)+1
296 C=0
297 DO 24 N = 3,30
298 C0=C
299 FN=N
300 FN2=FN+FN
301 FN1=FN2-1
302 FN3=FN1/(FN2-3)
303 FN4=12*FN**2-(1-XN)
304 FN5=16*FN1*X
305 RAN=1/((FN2+1)**2-XN)
306 F1=FN3*(FN4-20*FN)+FN5
307 F2=28*FN-FN4-8+FN5
308 F3=FN3*((FN2-5)**2-XN)
309 A=(F1*A2+F2*A1+F3*A0)*RAN
310 B=(F1*B2+F2*B1+F3*B0)*RAN
311 C=A/B
312 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
313 A0=A1
314 A1=A2
315 A2=A
316 B0=B1
317 B1=B2
318 B2=B
319 24 CONTINUE
320 25 S=C/SQRT(RPIH*X)
321 IF(.NOT.LEX) S=EXP(-X)*S
322 ELSE
323 R=1/X
324 H=10*R-1
325 ALFA=H+H
326 B1=0
327 B2=0
328 DO 13 I = 15,0,-1
329 B0=CC(I,MU)+ALFA*B1-B2
330 B2=B1
331 13 B1=B0
332 S=SQRT(PIH*R)*(B0-H*B2)
333 IF(.NOT.LEX) S=EXP(-X)*S
334 END IF
335#if !defined(CERNLIB_DOUBLE)
336 99 BSIR3=S
337#endif
338#if defined(CERNLIB_DOUBLE)
339 99 DBSIR3=S
340#endif
341 RETURN
342 101 FORMAT('INCORRECT ARGUMENT OR INDEX, X = ',1P,E15.6,' NU = ',I5)
343 END