]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/jimmy/cernmisc/jmdbska.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Dummy codes for CTEQ6 pdfs
[u/mrichter/AliRoot.git] / HERWIG / jimmy / cernmisc / jmdbska.F
CommitLineData
ef94df36 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1 2005/07/29 16:15:46 jmb
6* Include the various CERNLIB functions jimmy needs
7*
8* Revision 1.1.1.1 1996/04/01 15:02:08 mclareni
9* Mathlib gen
10*
11*
12 SUBROUTINE JMDBSKA(X,IA,JA,NL,B)
13 IMPLICIT DOUBLE PRECISION (A-H,J,O-Z)
14 SAVE
15 LOGICAL LEX
16 INTEGER JA
17
18 CHARACTER NAME*(*),ENAM*(*)
19 CHARACTER*80 ERRTXT
20 PARAMETER (NAME = 'BSKA/JMDBSKA', ENAM = 'EBSKA/DEBKA')
21
22 PARAMETER (Z1 = 1, Z2 = 2, Z3 = 3, Z4 = 4)
23 PARAMETER (Z12 = Z1/Z2, Z13 = Z1/Z3, Z14 = Z1/Z4, Z23 = Z2/Z3)
24 PARAMETER (Z34 = Z3/Z4)
25
26 DIMENSION B(0:*)
27
28 PARAMETER (PI = 3.14159 26535 89793D0, PIV = PI/4)
29
30 LEX=.FALSE.
31 GO TO 9
32
33 ENTRY JMDEBKA(X,IA,JA,NL,B)
34 LEX=.TRUE.
35
36 9 MODE=10*IA+JA
37 N=NL-1
38 U=2/X
39 IF(LEX) THEN
40 IF(X .LE. 0) THEN
41 N=0
42 WRITE(ERRTXT,101) X
43 CALL JMMLPT(ENAM,'C341.1',ERRTXT)
44 ELSEIF(NL .LT. 0 .OR. NL .GT. 100) THEN
45 N=0
46 WRITE(ERRTXT,103) NL
47 CALL JMMLPT(ENAM,'C341.3',ERRTXT)
48 ELSEIF(IA .EQ. 0) THEN
49 A=0
50 B(0)=JMDEBSK0(X)
51 B(1)=JMDEBSK1(X)
52 ELSEIF(MODE .EQ. 12) THEN
53 A=Z12
54 B(0)=SQRT(PIV*U)
55 B(1)=B(0)*(1+A*U)
56 ELSEIF(MODE .EQ. 13) THEN
57 A=Z13
58 B(0)=JMDEBKR3(X,1)
59 B(1)=JMDEBKR3(X,2)+A*U*B(0)
60 ELSEIF(MODE .EQ. 14) THEN
61 A=Z14
62 B(0)=JMDEBKR4(X,1)
63 B(1)=JMDEBKR4(X,3)+A*U*B(0)
64 ELSEIF(MODE .EQ. 23) THEN
65 A=Z23
66 B(0)=JMDEBKR3(X,2)
67 B(1)=JMDEBKR3(X,1)+A*U*B(0)
68 ELSEIF(MODE .EQ. 34) THEN
69 A=Z34
70 B(0)=JMDEBKR4(X,3)
71 B(1)=JMDEBKR4(X,1)+A*U*B(0)
72 ELSE
73 N=0
74 WRITE(ERRTXT,102) IA,JA
75 CALL JMMLPT(ENAM,'C341.2',ERRTXT)
76 ENDIF
77 ELSE
78 IF(X .LE. 0) THEN
79 N=0
80 WRITE(ERRTXT,101) X
81 CALL JMMLPT(NAME,'C341.1',ERRTXT)
82 ELSEIF(NL .LT. 0 .OR. NL .GT. 100) THEN
83 N=0
84 WRITE(ERRTXT,103) NL
85 CALL JMMLPT(NAME,'C341.3',ERRTXT)
86 ELSEIF(IA .EQ. 0) THEN
87 A=0
88 B(0)=JMDBESK0(X)
89 B(1)=JMDBESK1(X)
90 ELSEIF(MODE .EQ. 12) THEN
91 A=Z12
92 B(0)=EXP(-X)*SQRT(PIV*U)
93 B(1)=B(0)*(1+A*U)
94 ELSEIF(MODE .EQ. 13) THEN
95 A=Z13
96 B(0)=JMDBSKR3(X,1)
97 B(1)=JMDBSKR3(X,2)+A*U*B(0)
98 ELSEIF(MODE .EQ. 14) THEN
99 A=Z14
100 B(0)=JMDBSKR4(X,1)
101 B(1)=JMDBSKR4(X,3)+A*U*B(0)
102 ELSEIF(MODE .EQ. 23) THEN
103 A=Z23
104 B(0)=JMDBSKR3(X,2)
105 B(1)=JMDBSKR3(X,1)+A*U*B(0)
106 ELSEIF(MODE .EQ. 34) THEN
107 A=Z34
108 B(0)=JMDBSKR4(X,3)
109 B(1)=JMDBSKR4(X,1)+A*U*B(0)
110 ELSE
111 N=0
112 WRITE(ERRTXT,102) IA,JA
113 CALL JMMLPT(NAME,'C341.2',ERRTXT)
114 ENDIF
115 ENDIF
116 DO 1 J = 1,N
117 A=A+1
118 1 B(J+1)=B(J-1)+A*U*B(J)
119 RETURN
120 101 FORMAT('NON-POSITIVE ARGUMENT X = ',E15.6)
121 102 FORMAT('PAIR (IA,JA) = (',I5,I5,') ILLEGAL')
122 103 FORMAT('ILLEGAL NL =',I5)
123 END
124
125
126* $Id$
127*
128* $Log$
129* Revision 1.1 2005/07/29 16:15:46 jmb
130* Include the various CERNLIB functions jimmy needs
131*
132* Revision 1.1.1.1 1996/04/01 15:02:00 mclareni
133* Mathlib gen
134*
135*
136 DOUBLE PRECISION FUNCTION JMDBESI0(X)
137 IMPLICIT DOUBLE PRECISION (A-H,J,O-Z)
138 SAVE
139
140 LOGICAL LEX
141 CHARACTER NAME0*(*),NAME1*(*),NAME0E*(*),NAME1E*(*)
142 CHARACTER*80 ERRTXT
143 DIMENSION CI(0:24,0:1),CK(0:16,0:1)
144
145 PARAMETER (NAME0 = 'JMBESK0/JMDBESK0', NAME0E =
146 $ 'JMEBESK0/JMDBESK0')
147 PARAMETER (NAME1 = 'JMBESK1/JMDBESK1', NAME1E =
148 $ 'JMEBESK1/JMDBESK1')
149
150 PARAMETER (EPS=1D-15)
151 PARAMETER (Z1 = 1, HF = Z1/2)
152 PARAMETER (PI = 3.14159 26535 89793 24D0)
153 PARAMETER (CE = 0.57721 56649 01532 86D0)
154 PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
155
156 DATA CI( 0,0) /+1.00827 92054 58740 032D0/
157 DATA CI( 1,0) /+0.00844 51226 24920 943D0/
158 DATA CI( 2,0) /+0.00017 27006 30777 567D0/
159 DATA CI( 3,0) /+0.00000 72475 91099 959D0/
160 DATA CI( 4,0) /+0.00000 05135 87726 878D0/
161 DATA CI( 5,0) /+0.00000 00568 16965 808D0/
162 DATA CI( 6,0) /+0.00000 00085 13091 223D0/
163 DATA CI( 7,0) /+0.00000 00012 38425 364D0/
164 DATA CI( 8,0) /+0.00000 00000 29801 672D0/
165 DATA CI( 9,0) /-0.00000 00000 78956 698D0/
166 DATA CI(10,0) /-0.00000 00000 33127 128D0/
167 DATA CI(11,0) /-0.00000 00000 04497 339D0/
168 DATA CI(12,0) /+0.00000 00000 01799 790D0/
169 DATA CI(13,0) /+0.00000 00000 00965 748D0/
170 DATA CI(14,0) /+0.00000 00000 00038 604D0/
171 DATA CI(15,0) /-0.00000 00000 00104 039D0/
172 DATA CI(16,0) /-0.00000 00000 00023 950D0/
173 DATA CI(17,0) /+0.00000 00000 00009 554D0/
174 DATA CI(18,0) /+0.00000 00000 00004 443D0/
175 DATA CI(19,0) /-0.00000 00000 00000 859D0/
176 DATA CI(20,0) /-0.00000 00000 00000 709D0/
177 DATA CI(21,0) /+0.00000 00000 00000 087D0/
178 DATA CI(22,0) /+0.00000 00000 00000 112D0/
179 DATA CI(23,0) /-0.00000 00000 00000 012D0/
180 DATA CI(24,0) /-0.00000 00000 00000 018D0/
181
182 DATA CI( 0,1) /+0.97580 06023 26285 926D0/
183 DATA CI( 1,1) /-0.02446 74429 63276 385D0/
184 DATA CI( 2,1) /-0.00027 72053 60763 829D0/
185 DATA CI( 3,1) /-0.00000 97321 46728 020D0/
186 DATA CI( 4,1) /-0.00000 06297 24238 640D0/
187 DATA CI( 5,1) /-0.00000 00659 61142 154D0/
188 DATA CI( 6,1) /-0.00000 00096 13872 919D0/
189 DATA CI( 7,1) /-0.00000 00014 01140 901D0/
190 DATA CI( 8,1) /-0.00000 00000 47563 167D0/
191 DATA CI( 9,1) /+0.00000 00000 81530 681D0/
192 DATA CI(10,1) /+0.00000 00000 35408 148D0/
193 DATA CI(11,1) /+0.00000 00000 05102 564D0/
194 DATA CI(12,1) /-0.00000 00000 01804 409D0/
195 DATA CI(13,1) /-0.00000 00000 01023 594D0/
196 DATA CI(14,1) /-0.00000 00000 00052 678D0/
197 DATA CI(15,1) /+0.00000 00000 00107 094D0/
198 DATA CI(16,1) /+0.00000 00000 00026 120D0/
199 DATA CI(17,1) /-0.00000 00000 00009 561D0/
200 DATA CI(18,1) /-0.00000 00000 00004 713D0/
201 DATA CI(19,1) /+0.00000 00000 00000 829D0/
202 DATA CI(20,1) /+0.00000 00000 00000 743D0/
203 DATA CI(21,1) /-0.00000 00000 00000 080D0/
204 DATA CI(22,1) /-0.00000 00000 00000 117D0/
205 DATA CI(23,1) /+0.00000 00000 00000 011D0/
206 DATA CI(24,1) /+0.00000 00000 00000 019D0/
207
208 DATA CK( 0,0) /+0.98840 81742 30825 800D0/
209 DATA CK( 1,0) /-0.01131 05046 46928 281D0/
210 DATA CK( 2,0) /+0.00026 95326 12762 724D0/
211 DATA CK( 3,0) /-0.00001 11066 85196 665D0/
212 DATA CK( 4,0) /+0.00000 06325 75108 500D0/
213 DATA CK( 5,0) /-0.00000 00450 47337 641D0/
214 DATA CK( 6,0) /+0.00000 00037 92996 456D0/
215 DATA CK( 7,0) /-0.00000 00003 64547 179D0/
216 DATA CK( 8,0) /+0.00000 00000 39043 756D0/
217 DATA CK( 9,0) /-0.00000 00000 04579 936D0/
218 DATA CK(10,0) /+0.00000 00000 00580 811D0/
219 DATA CK(11,0) /-0.00000 00000 00078 832D0/
220 DATA CK(12,0) /+0.00000 00000 00011 360D0/
221 DATA CK(13,0) /-0.00000 00000 00001 727D0/
222 DATA CK(14,0) /+0.00000 00000 00000 275D0/
223 DATA CK(15,0) /-0.00000 00000 00000 046D0/
224 DATA CK(16,0) /+0.00000 00000 00000 008D0/
225
226 DATA CK( 0,1) /+1.03595 08587 72358 331D0/
227 DATA CK( 1,1) /+0.03546 52912 43331 114D0/
228 DATA CK( 2,1) /-0.00046 84750 28166 889D0/
229 DATA CK( 3,1) /+0.00001 61850 63810 053D0/
230 DATA CK( 4,1) /-0.00000 08451 72048 124D0/
231 DATA CK( 5,1) /+0.00000 00571 32218 103D0/
232 DATA CK( 6,1) /-0.00000 00046 45554 607D0/
233 DATA CK( 7,1) /+0.00000 00004 35417 339D0/
234 DATA CK( 8,1) /-0.00000 00000 45757 297D0/
235 DATA CK( 9,1) /+0.00000 00000 05288 133D0/
236 DATA CK(10,1) /-0.00000 00000 00662 613D0/
237 DATA CK(11,1) /+0.00000 00000 00089 048D0/
238 DATA CK(12,1) /-0.00000 00000 00012 726D0/
239 DATA CK(13,1) /+0.00000 00000 00001 921D0/
240 DATA CK(14,1) /-0.00000 00000 00000 305D0/
241 DATA CK(15,1) /+0.00000 00000 00000 050D0/
242 DATA CK(16,1) /-0.00000 00000 00000 009D0/
243
244 NU=0
245 LEX=.FALSE.
246 GO TO 6
247
248 ENTRY JMDEBSI0(X)
249
250 NU=0
251 LEX=.TRUE.
252 GO TO 6
253
254 ENTRY JMDBESI1(X)
255
256 NU=1
257 LEX=.FALSE.
258 GO TO 6
259
260 ENTRY JMDEBSI1(X)
261
262 NU=1
263 LEX=.TRUE.
264
265 6 V=ABS(X)
266 IF(V .LT. 8) THEN
267 Y=(HF*V)**2
268 XL=NU+2
269 A0=1
270 A1=1+2*Y/((XL+1)*(XL-1))
271 A2=1+Y*(4+3*Y/((XL+2)*XL))/((XL+3)*(XL-1))
272 B0=1
273 B1=1-Y/(XL+1)
274 B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
275 W1=3+XL
276 V1=3-XL
277 V3=XL-1
278 V2=V3+V3
279 C=0
280 DO 3 N = 3,30
281 C0=C
282 FN=N
283 W1=W1+2
284 W2=W1-1
285 W3=W2-1
286 W4=W3-1
287 W5=W4-1
288 W6=W5-1
289 V1=V1+1
290 V2=V2+1
291 V3=V3+1
292 U1=FN*W4
293 E=V3/(U1*W3)
294 U2=E*Y
295 F1=1+Y*V1/(U1*W1)
296 F2=(1+Y*V2/(V3*W2*W5))*U2
297 F3=-Y*Y*U2/(W4*W5*W5*W6)
298 A=F1*A2+F2*A1+F3*A0
299 B=F1*B2+F2*B1+F3*B0
300 C=A/B
301 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 4
302 A0=A1
303 A1=A2
304 A2=A
305 B0=B1
306 B1=B2
307 B2=B
308 3 CONTINUE
309 4 H=C
310 IF(NU .EQ. 1) H=HF*X*H
311 IF(LEX) H=EXP(-V)*H
312 ELSE
313 R=1/V
314 H=16*R-1
315 ALFA=H+H
316 B1=0
317 B2=0
318 DO 1 I = 24,0,-1
319 B0=CI(I,NU)+ALFA*B1-B2
320 B2=B1
321 1 B1=B0
322 H=SQRT(RPI2*R)*(B0-H*B2)
323 IF(NU*X .LT. 0) H=-H
324 IF(.NOT.LEX) H=EXP(V)*H
325 ENDIF
326 GO TO 9
327
328 ENTRY JMDBESK0(X)
329
330 NU=0
331 LEX=.FALSE.
332 GO TO 8
333
334 ENTRY JMDEBSK0(X)
335
336 NU=0
337 LEX=.TRUE.
338 GO TO 8
339
340 ENTRY JMDBESK1(X)
341
342 NU=1
343 LEX=.FALSE.
344 GO TO 8
345
346 ENTRY JMDEBSK1(X)
347
348 NU=1
349 LEX=.TRUE.
350
351 8 IF(X .LE. 0) THEN
352 H=0
353 WRITE(ERRTXT,101) X
354 IF(NU .EQ. 0 .AND. .NOT.LEX) CALL JMMLPT(NAME0 ,'C313.1',ERRTXT)
355 IF(NU .EQ. 0 .AND. LEX) CALL JMMLPT(NAME0E,'C313.1',ERRTXT)
356 IF(NU .EQ. 1 .AND. .NOT.LEX) CALL JMMLPT(NAME1 ,'C313.1',ERRTXT)
357 IF(NU .EQ. 1 .AND. LEX) CALL JMMLPT(NAME1E,'C313.1',ERRTXT)
358 ELSEIF(X .LT. 1) THEN
359 B=HF*X
360 BK=-(LOG(B)+CE)
361 F=BK
362 P=HF
363 Q=HF
364 C=1
365 D=B**2
366 BK1=P
367 DO 11 N = 1,15
368 FN=N
369 RFN=1/FN
370 P=P*RFN
371 Q=Q*RFN
372 F=(F+P+Q)*RFN
373 C=C*D*RFN
374 G=C*(P-FN*F)
375 H=C*F
376 BK=BK+H
377 BK1=BK1+G
378 IF(BK1*H+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
379 11 CONTINUE
380 12 H=BK
381 IF(NU .EQ. 1) H=BK1/B
382 IF(LEX) H=EXP(X)*H
383 ELSEIF(X .LE. 5) THEN
384 XN=4*NU**2
385 A=9-XN
386 B=25-XN
387 C=768*X**2
388 C0=48*X
389 A0=1
390 A1=(16*X+7+XN)/A
391 A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
392 B0=1
393 B1=(16*X+9-XN)/A
394 B2=(C+C0*B)/(A*B)+1
395 C=0
396 DO 24 N = 3,30
397 C0=C
398 FN=N
399 FN2=FN+FN
400 FN1=FN2-1
401 FN3=FN1/(FN2-3)
402 FN4=12*FN**2-(1-XN)
403 FN5=16*FN1*X
404 RAN=1/((FN2+1)**2-XN)
405 F1=FN3*(FN4-20*FN)+FN5
406 F2=28*FN-FN4-8+FN5
407 F3=FN3*((FN2-5)**2-XN)
408 A=(F1*A2+F2*A1+F3*A0)*RAN
409 B=(F1*B2+F2*B1+F3*B0)*RAN
410 C=A/B
411 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
412 A0=A1
413 A1=A2
414 A2=A
415 B0=B1
416 B1=B2
417 B2=B
418 24 CONTINUE
419 25 H=C/SQRT(RPIH*X)
420 IF(.NOT.LEX) H=EXP(-X)*H
421 ELSE
422 R=1/X
423 H=10*R-1
424 ALFA=H+H
425 B1=0
426 B2=0
427 DO 23 I = 16,0,-1
428 B0=CK(I,NU)+ALFA*B1-B2
429 B2=B1
430 23 B1=B0
431 H=SQRT(PIH*R)*(B0-H*B2)
432 IF(.NOT.LEX) H=EXP(-X)*H
433 ENDIF
434 9 CONTINUE
435 JMDBESI0=H
436 RETURN
437 101 FORMAT(' NON-POSITIVE ARGUMENT X = ',1P,E15.6)
438 END