]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.3.1/wrapacfgpg.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
end-of-line normalization
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapacfgpg.f
CommitLineData
4e9e3152 1 subroutine ACFGPevolvep(xin,qin,p2in,ip2in,pdf)
2 include 'parmsetup.inc'
3 real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5
4 real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu
5 real*4 par,par2,calc,celc
6 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
7 character*16 name(nmxset)
8 integer nmem(nmxset),ndef(nmxset),mmem
9 common/NAME/name,nmem,ndef,mmem
10 integer nset
11
12 save
13
14 if(imem.eq.1) then
15 call ACFGP1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
16
17 elseif(imem.eq.2) then
18 call ACFGP2(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
19
20 elseif(imem.eq.3.or.imem.eq.0) then
21 call SFAFG1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
22
23 else
24 CONTINUE
25 endif
26
27 pdf(-6)= 0.0d0
28 pdf(6)= 0.0d0
29 pdf(-5)= 0.0d0
30 pdf(5 )= 0.0d0
31 pdf(-4)= chm
32 pdf(4 )= chm
33 pdf(-3)= str
34 pdf(3 )= str
35 pdf(-2)= usea
36 pdf(2 )= upv
37 pdf(-1)= dsea
38 pdf(1 )= dnv
39 pdf(0 )= glu
40
41 return
42ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
43 entry ACFGPread(nset)
44 read(1,*)nmem(nset),ndef(nset)
45 do iset = 1,3
46 read(1,*)(par(j,iset),j=1,4)
47 read(1,*)(par(j,iset),j=5,8)
48 read(1,*)(par(j,iset),j=9,12)
49 read(1,*)(par(j,iset),j=13,16)
50 read(1,*)(par(j,iset),j=17,20)
51 read(1,*)(par(j,iset),j=21,24)
52 read(1,*)(par(j,iset),j=25,28)
53 read(1,*)(par(j,iset),j=29,30)
54 do j=1,32
55 do k=1,20
56 read(1,*)(calc(i,k,j,iset),i=1,4)
57 read(1,*)(calc(i,k,j,iset),i=5,8)
58 enddo
59 enddo
60 enddo
61c last one
62 read(1,*)(par2(j),j=1,4)
63 read(1,*)(par2(j),j=5,8)
64 read(1,*)(par2(j),j=9,12)
65 read(1,*)(par2(j),j=13,16)
66 read(1,*)(par2(j),j=17,20)
67 read(1,*)(par2(j),j=21,24)
68 read(1,*)(par2(j),j=25,28)
69 read(1,*)(par2(j),j=29,30)
70 do j=1,32
71 do k=1,20
72 read(1,*)(celc(i,k,j),i=1,4)
73 read(1,*)(celc(i,k,j),i=5,8)
74 enddo
75 enddo
76
77 return
78
79
80c
81ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
82 entry ACFGPalfa(alfas,qalfa)
83 call getnset(iset)
84 call GetOrderAsM(iset,iord)
85 call Getlam4M(iset,imem,qcdl4)
86 call Getlam5M(iset,imem,qcdl5)
87 call aspdflib(alfas,Qalfa,iord,qcdl5)
88
89 return
90c
91ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
92 entry ACFGPinit(Eorder,Q2fit)
93 return
94c
95ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
96 entry ACFGPpdf(mem)
97 imem = mem
98 return
99c
100 1000 format(5e13.5)
101 end
102ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
103 SUBROUTINE ACFGP1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
104C
105C INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
106C QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
107C
108C THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
109C CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
110C
111C THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
112C BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
113C
114C THE PROGRAM WORKS FOR 2. GEV**2 < Q2 <5.5E+5 AND .00137 < X < .9986
115C
116C THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
117C
118C THE VALUE OF LAMBDA-MSB IS 200 MEV
119C
120C THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
121C X*U=X*U(X,Q2)
122C X*D= ...
123C X*S= ...
124C X*C= ... (MASSLESS CHARM WITH C(X,2.)=0)
125C X*CM= ... (MASSIVE CHARM WITH MC=1.5 GEV )
126C X*GLU=GLUON(X,Q2)*X
127C
128C
129C F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
130C F2C= " " " WITH MASSIVE CHARM
131C
132 double precision
133 + DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
134 REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL
135 REAL Q2
136 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
137 DIMENSION XPDF(7),CALCO(8,20,32)
138 COMMON/W5051I7/CALCO
139 EXTERNAL AFCPLU
140 DATA ZERO/0.0/
141C----------------------------------------------------------------------
142 DATA ISTART/0/
143 SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
144C
145 IF (ISTART.EQ.0) THEN
146 ISTART=1
147 DO 10 K=1,32
148 DO 10 I=1,20
149 DO 10 M=1,8
150 10 CALCO(M,I,K) = CALC(M,I,K,1)
151 OWLAM=PAR(1,1)
152 OWLAM2=OWLAM**2
153 Q02=PAR(30,1)
154 FLAV=PAR(25,1)
155 DELTA=PAR(29,1)
156 CALL WATE32
157 ENDIF
158C
159 X = DX
160 Q = DQ
161 Q2 = Q*Q
162 IDQ2=2
163 SB=0.
164 IF(Q2-Q02) 1,1,2
165 2 IF(IDQ2-1) 1,1,3
166 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
167 1 CONTINUE
168 CALL AURGAM(8,0,X,SB,XPDF(7))
169 CALL AURGAM(7,0,X,SB,SING)
170 CALL AURGAM(4,0,X,SB,DPLUSNS)
171 CALL AURGAM(3,0,X,SB,CPLUSNS)
172 CALL AURGAM(5,0,X,SB,UPLUSNS)
173 CALL AURGAM(6,0,X,SB,SPLUSNS)
174 XPDF(3) = CPLUSNS
175 XPDF(4) = DPLUSNS
176 XPDF(5) = UPLUSNS
177 XPDF(6) = SPLUSNS
178 XPDF(1) = SING
179C
180 ADD = XPDF(1)/FLAV
181 UPLUS=XPDF(5)+ADD
182 DPLUS=-XPDF(4)+ADD
183 SPLUS=-XPDF(6)+ADD
184 CPLUS=-XPDF(3)+ADD
185 UB=UPLUS*0.5
186 UV=UB
187 DB=DPLUS*0.5
188 DV=DB
189 SB=SPLUS*0.5
190 CB=CPLUS*0.5
191 SING=XPDF(1)
192 GLU=XPDF(7)
193 GL=GLU
194*
195 DUV=MAX(ZERO,UV)
196 DDV=MAX(ZERO,DV)
197 DUB=MAX(ZERO,UB)
198 DDB=MAX(ZERO,DB)
199 DSB=MAX(ZERO,SB)
200 DCB=MAX(ZERO,CB)
201 DGL=MAX(ZERO,GL)
202C
203 RETURN
204 END
205cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
206 SUBROUTINE ACFGP2(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
207C
208C INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
209C QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
210C
211C THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
212C CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
213C
214C THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
215C BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
216C
217C THE PROGRAM WORKS FOR 2. GEV**2 < Q2 <5.5E+5 AND .00137 < X < .9986
218C
219C THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
220C
221C THE VALUE OF LAMBDA-MSB IS 200 MEV
222C
223C THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
224C X*U=X*U(X,Q2)
225C X*D= ...
226C X*S= ...
227C X*C= ... (MASSLESS CHARM WITH C(X,2.)=0)
228C X*CM= ... (MASSIVE CHARM WITH MC=1.5 GEV )
229C X*GLU=GLUON(X,Q2)*X
230C
231C
232C F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
233C F2C= " " " WITH MASSIVE CHARM
234C
235 double precision
236 + DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
237 REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL
238 REAL Q2
239 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
240 DIMENSION XPDF(7),CALCO(8,20,32)
241 COMMON/W5051I7/CALCO
242 EXTERNAL AFCPLU
243 DATA ZERO/0.0/
244C----------------------------------------------------------------------
245 DATA ISTART/0/
246 SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
247C
248 IF (ISTART.EQ.0) THEN
249 ISTART=1
250 DO 10 K=1,32
251 DO 10 I=1,20
252 DO 10 M=1,8
253 10 CALCO(M,I,K) = CALC(M,I,K,2)
254 OWLAM=PAR(1,2)
255 OWLAM2=OWLAM**2
256 Q02=PAR(30,2)
257 FLAV=PAR(25,2)
258 DELTA=PAR(29,2)
259 CALL WATE32
260 ENDIF
261C
262 X = DX
263 Q = DQ
264 Q2 = Q*Q
265 IDQ2=2
266 SB=0.
267 IF(Q2-Q02) 1,1,2
268 2 IF(IDQ2-1) 1,1,3
269 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
270 1 CONTINUE
271 CALL AURGAM(8,0,X,SB,XPDF(7))
272 CALL AURGAM(7,0,X,SB,SING)
273 CALL AURGAM(4,0,X,SB,DPLUSNS)
274 CALL AURGAM(3,0,X,SB,CPLUSNS)
275 CALL AURGAM(5,0,X,SB,UPLUSNS)
276 CALL AURGAM(6,0,X,SB,SPLUSNS)
277 XPDF(3) = CPLUSNS
278 XPDF(4) = DPLUSNS
279 XPDF(5) = UPLUSNS
280 XPDF(6) = SPLUSNS
281 XPDF(1) = SING
282C
283 ADD = XPDF(1)/FLAV
284 UPLUS=XPDF(5)+ADD
285 DPLUS=-XPDF(4)+ADD
286 SPLUS=-XPDF(6)+ADD
287 CPLUS=-XPDF(3)+ADD
288 UB=UPLUS*0.5
289 UV=UB
290 DB=DPLUS*0.5
291 DV=DB
292 SB=SPLUS*0.5
293 CB=CPLUS*0.5
294 SING=XPDF(1)
295 GLU=XPDF(7)
296 GL=GLU
297C... get parton density with massive charm
298 CPLUM=AFCPLU(X,Q2)
299 CB=CPLUM*0.5
300*
301 DUV=MAX(ZERO,UV)
302 DDV=MAX(ZERO,DV)
303 DUB=MAX(ZERO,UB)
304 DDB=MAX(ZERO,DB)
305 DSB=MAX(ZERO,SB)
306 DCB=MAX(ZERO,CB)
307 DGL=MAX(ZERO,GL)
308C
309 RETURN
310 END
311ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
312 SUBROUTINE SFAFG1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
313C
314C**************************************************************************
315C ( 1st of February 1994)
316C This is an interpolation program which reads the files GRPOL and
317C GRVDM and gives the quark and gluon distributions in real photon
318C as functions of x and Q**2.
319C
320C The Q**2 evolution is a BLL evolution (MSbar scheme) with Nf=4
321C and LAMBDA(MSbar)=.200 Gev.
322C
323C A massless charm distribution is generated for Q**2 > 2 Gev**2.
324C
325C The distributions are the sum of a pointlike part (PL) and of a
326C Vdm part (VDM):
327C dist=PL + KA*VDM
328C KA is a factor which can be adjusted ( the default value is KA=1.0).
329C The file GRPOL contains the pointlike part of the distributions.
330C The file GRVDM contains the vdm part (A precise definition of this
331C latter is given in the paper "PARTON DISTRIBUTIONS IN THE PHOTON",
332C Preprint LPTHE Orsay 93-37, by P.Aurenche,M.Fontannaz and J.Ph.Guillet).
333C
334C The output of the program is written in the file GETOUT with the
335C following conventions
336C UPLUS=x(u+ubar)
337C DPLUS=x(d+dbar)
338C SPLUS=x(s+sbar)
339C CPLUS=x(c+cbar)
340C SING =UPLUS+DPLUS+SPLUS+CPLUS
341C GLU =x*g
342C
343C The interpolation is valid for 2. < Q**2 < 5.5E+5 Gev**2,
344C and for .0015< x < .99
345C
346C The program also gives the structure function F2:
347C F2 = q*Cq + g*Cg + Cgam
348C Cq and Cg are the Wilson coeficients and Cgam is the direct term.
349C
350C Although the charm quark evolution is massless, the direct term
351C Cgam includes the effects due to the charm quark mass. The charm
352C quark threshold is therefore correctly described at the lowest
353C ordre in alphastrong (Details are given in the preprint).
354C
355C
356C**************************************************************************
357C
358 double precision
359 + DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL
360 REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL
361 REAL Q2
362 DIMENSION XPDF(7)
363 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32)
364 DIMENSION CALCO(8,20,32)
365 DIMENSION CELCO(8,20,32)
366 COMMON/W5051IA/CALCO
367 COMMON/W5051IB/CELCO
368 EXTERNAL AFCPLU
369 DATA ZERO/0.0/
370C----------------------------------------------------------------------
371 DATA ISTART/0/
372 SAVE ISTART,OWLAM2,Q02,FLAV,KA, /W5051IA/, /W5051IB/
373C
374 IF (ISTART.EQ.0) THEN
375 ISTART=1
376 DO 10 K=1,32
377 DO 10 I=1,20
378 DO 10 M=1,8
379 CALCO(M,I,K) = CALC(M,I,K,3)
380 10 CELCO(M,I,K) = CELC(M,I,K)
381 OWLAM=PAR(1,3)
382 OWLAM2=OWLAM**2
383 Q02=PAR(30,3)
384 FLAV=PAR(25,3)
385 DELTA=PAR(29,3)
386 CALL WATE32
387 KA=1.0
388 ENDIF
389C
390 X = DX
391 Q = DQ
392 Q2 = Q*Q
393 IDQ2=2
394 SB=0.
395 IF(Q2-Q02) 1,1,2
396 2 IF(IDQ2-1) 1,1,3
397 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2))
398 1 CONTINUE
399 CALL AFGINT(8,0,X,SB,XPDF(7))
400 CALL AFGINT(7,0,X,SB,SING)
401 CALL AFGINT(4,0,X,SB,DPLUSNS)
402 CALL AFGINT(3,0,X,SB,CPLUSNS)
403 CALL AFGINT(5,0,X,SB,UPLUSNS)
404 CALL AFGINT(6,0,X,SB,SPLUSNS)
405 XPDF(3) = CPLUSNS
406 XPDF(4) = DPLUSNS
407 XPDF(5) = UPLUSNS
408 XPDF(6) = SPLUSNS
409 XPDF(1) = SING
410C
411 ADD = XPDF(1)/FLAV
412 UPLUS= XPDF(5)+ADD
413 DPLUS=-XPDF(4)+ADD
414 SPLUS=-XPDF(6)+ADD
415 CPLUS=-XPDF(3)+ADD
416 SING=XPDF(1)
417 GLU=XPDF(7)
418 GL=GLU
419*
420 CALL AFGIN2(8,0,X,SB,XPDF(7))
421 CALL AFGIN2(7,0,X,SB,SING)
422 CALL AFGIN2(4,0,X,SB,DPLUSNS)
423 CALL AFGIN2(3,0,X,SB,CPLUSNS)
424 CALL AFGIN2(5,0,X,SB,UPLUSNS)
425 CALL AFGIN2(6,0,X,SB,SPLUSNS)
426 XPDF(3) = CPLUSNS
427 XPDF(4) = DPLUSNS
428 XPDF(5) = UPLUSNS
429 XPDF(6) = SPLUSNS
430 XPDF(1) = SING
431C
432 ADD2 = XPDF(1)/FLAV
433 UPLU2= XPDF(5)+ADD2
434 DPLU2=-XPDF(4)+ADD2
435 SPLU2=-XPDF(6)+ADD2
436 CPLU2=-XPDF(3)+ADD2
437 SING2=XPDF(1)
438 GLU2=XPDF(7)
439 UB=UPLUS+UPLU2*KA
440 UV=UB
441 DB=DPLUS+DPLU2*KA
442 DV=DB
443 SB=SPLUS+SPLU2*KA
444 CB=CPLUS+CPLU2*KA
445 SING=SING+SING2*KA
446 GL=GLU+GLU2*KA
447*
448 DUV=MAX(ZERO,UV)
449 DDV=MAX(ZERO,DV)
450 DUB=MAX(ZERO,UB)
451 DDB=MAX(ZERO,DB)
452 DSB=MAX(ZERO,SB)
453 DCB=MAX(ZERO,CB)
454 DGL=MAX(ZERO,GL)
455C
456 RETURN
457 END
458cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
459 SUBROUTINE WATE32
460C 32 POINT GAUSSIAN QUADRATURE ROUTINE
461 double precision
462 + X(16),W(16)
463 double precision
464 + XI(32),WI(32),XX(33)
465 COMMON/W5051I9/XI,WI,XX,NTERMS
466 NTERMS=32
467 X(1)=0.048307665687738316235D0
468 X(2)=0.144471961582796493485D0
469 X(3)=0.239287362252137074545D0
470 X(4)=0.331868602282127649780D0
471 X(5)=0.421351276130635345364D0
472 X(6)=0.506899908932229390024D0
473 X(7)=0.587715757240762329041D0
474 X(8)=0.663044266930215200975D0
475 X(9)=0.732182118740289680387D0
476 X(10)=0.794483795967942406963D0
477 X(11)=0.849367613732569970134D0
478 X(12)=0.896321155766052123965D0
479 X(13)=0.934906075937739689171D0
480 X(14)=0.964762255587506430774D0
481 X(15)=0.985611511545268335400D0
482 X(16)=0.997263861849481563545D0
483 W(1)=0.096540088514727800567D0
484 W(2)=0.095638720079274859419D0
485 W(3)=0.093844399080804565639D0
486 W(4)=0.091173878695763884713D0
487 W(5)=0.087652093004403811143D0
488 W(6)=0.083311924226946755222D0
489 W(7)=0.078193895787070306472D0
490 W(8)=0.072345794108848506225D0
491 W(9)=0.065822222776361846838D0
492 W(10)=0.058684093478535547145D0
493 W(11)=0.050998059262376176196D0
494 W(12)=0.042835898022226680657D0
495 W(13)=0.034273862913021433103D0
496 W(14)=0.025392065309262059456D0
497 W(15)=0.016274394730905670605D0
498 W(16)=0.007018610009470096600D0
499 NTERMH = NTERMS/2
500 DO 1 I=1,NTERMH
501 XI(I)=-X(17-I)
502 WI(I)=W(17-I)
503 XI(I+16)=X(I)
504 WI(I+16)=W(I)
505 1 CONTINUE
506 DO 2 I=1,NTERMS
507 2 XX(I)=0.5D0*(XI(I)+1.0D0)
508 XX(33)=1.0D0
509 RETURN
510 END
511ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
512 SUBROUTINE AURGAM(I,NDRV,X,S,ANS)
513 DIMENSION F1(32),F2(32),F3(32)
514 DIMENSION AF(3),AS(3)
515 DIMENSION CALCO(8,20,32)
516 COMMON/W5051I7/CALCO
517 DATA DELTA/0.8000E-01/
518 ANS=0.
519 IF(X.GT.0.9985) RETURN
520 N=3
521 IS=S/DELTA+1
522 IF(IS.GE.17) IS=17
523 IS1=IS+1
524 IS2=IS1+1
525 DO 1 L=1,32
526 KL=L+32*NDRV
527 F1(L)=CALCO(I,IS,KL)
528 F2(L)=CALCO(I,IS1,KL)
529 F3(L)=CALCO(I,IS2,KL)
530 1 CONTINUE
531 AF(1)=AFGETFV(X,F1)
532 AF(2)=AFGETFV(X,F2)
533 AF(3)=AFGETFV(X,F3)
534 AS(1)=(IS-1)*DELTA
535 AS(2)=AS(1)+DELTA
536 AS(3)=AS(2)+DELTA
537 CALL AFPOLIN(AS,AF,N,S,AANS,DY)
538 ANS=AANS
539 RETURN
540 END
541ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
542 SUBROUTINE AFGINT(I,NDRV,X,S,ANS)
543 DIMENSION F1(32),F2(32),F3(32)
544 DIMENSION AF(3),AS(3)
545 DIMENSION CALCO(8,20,32)
546 COMMON/W5051IA/CALCO
547 DATA DELTA/0.8000E-01/
548 ANS=0.
549 IF(X.GT.0.9985) RETURN
550 N=3
551 IS=S/DELTA+1
552C IF(IS.GE.17) IS=17
553 IS1=IS+1
554 IS2=IS1+1
555 DO 1 L=1,32
556 KL=L+32*NDRV
557 F1(L)=CALCO(I,IS,KL)
558 F2(L)=CALCO(I,IS1,KL)
559 F3(L)=CALCO(I,IS2,KL)
560 1 CONTINUE
561 AF(1)=AFGETFV(X,F1)
562 AF(2)=AFGETFV(X,F2)
563 AF(3)=AFGETFV(X,F3)
564 AS(1)=(IS-1)*DELTA
565 AS(2)=AS(1)+DELTA
566 AS(3)=AS(2)+DELTA
567 CALL AFPOLIN(AS,AF,N,S,AANS,DY)
568 ANS=AANS
569 RETURN
570 END
571ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
572 SUBROUTINE AFGIN2(I,NDRV,X,S,ANS)
573 DIMENSION F1(32),F2(32),F3(32)
574 DIMENSION AF(3),AS(3)
575 DIMENSION CELCO(8,20,32)
576 COMMON/W5051IB/CELCO
577 DATA DELTA/0.8000E-01/
578 ANS=0.
579 IF(X.GT.0.9985) RETURN
580 N=3
581 IS=S/DELTA+1
582C IF(IS.GE.17) IS=17
583 IS1=IS+1
584 IS2=IS1+1
585 DO 1 L=1,32
586 KL=L+32*NDRV
587 F1(L)=CELCO(I,IS,KL)
588 F2(L)=CELCO(I,IS1,KL)
589 F3(L)=CELCO(I,IS2,KL)
590 1 CONTINUE
591 AF(1)=AFGETFV(X,F1)
592 AF(2)=AFGETFV(X,F2)
593 AF(3)=AFGETFV(X,F3)
594 AS(1)=(IS-1)*DELTA
595 AS(2)=AS(1)+DELTA
596 AS(3)=AS(2)+DELTA
597 CALL AFPOLIN(AS,AF,N,S,AANS,DY)
598 ANS=AANS
599 RETURN
600 END
601ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
602 FUNCTION AFCPLU(X,Q2)
603 CMS=1.5**2
604 BETS=1-4.*CMS*X/(1.-X)/Q2
605 IF(BETS.LE..0) THEN
606 AFCPLU=.0
607 RETURN
608 ENDIF
609 BETA=SQRT(BETS)
610 CPLU=(8.*X*(1.-X)-1.-4.*CMS*X*(1.-X)/Q2)*BETA
611 CAU=X**2+(1.-X)**2+4.*CMS*X*(1.-3.*X)/Q2-8.*CMS**2*X**2/Q2**2
612 CPLU=CPLU+CAU* LOG((1.+BETA)/(1.-BETA))
613 AFCPLU=3.*(4./9.)*CPLU*X/(3.1415*137.)
614 1 RETURN
615 END
616ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
617 FUNCTION AFGETFV(X,FVL)
618C NOUVEAU PROGRAMME D'INTERPOLATION UTILISANT UNE ROUTINE DE MATH. RECIPES
619 DIMENSION FVL(32)
620 double precision
621 + XI(32),WI(32),XX(33)
622 COMMON/W5051I9/XI,WI,XX,NTERMS
623 DIMENSION A(4),B(4)
624 N=4
625 EPS=1.E-7
626 XAM=XX(1)-EPS
627 XAP=XX(1)+EPS
628C IF(X.LT.XAM) PRINT*,' X = ',X
629 IF(X.GT.XAM.AND.X.LT.XAP) GO TO 50
630 GO TO 80
631 50 Y=FVL(1)
632 GO TO 77
633 80 IF(X.LT.XX(2)) GO TO 51
634 IF(X.GT.XX(30)) GO TO 61
635 DO 1 I=3,30
636 IF(X.GT.XX(I)) GO TO 1
637 A(1)=XX(I-2)
638 A(2)=XX(I-1)
639 A(3)=XX(I)
640 A(4)=XX(I+1)
641 B(1)=FVL(I-2)
642 B(2)=FVL(I-1)
643 B(3)=FVL(I)
644 B(4)=FVL(I+1)
645 GO TO 70
646 1 CONTINUE
647 61 A(1)=XX(29)
648 A(2)=XX(30)
649 A(3)=XX(31)
650 A(4)=XX(32)
651 B(1)=FVL(29)
652 B(2)=FVL(30)
653 B(3)=FVL(31)
654 B(4)=FVL(32)
655 GO TO 70
656 51 A(1)=XX(1)
657 A(2)=XX(2)
658 A(3)=XX(3)
659 A(4)=XX(4)
660 B(1)=FVL(1)
661 B(2)=FVL(2)
662 B(3)=FVL(3)
663 B(4)=FVL(4)
664C 70 IF(X.GT..2.AND.X.LT..8) THEN
665C CALL AFPOLIN(A,B,N,X,Y,DY)
666C ELSE
667C CALL AFRATIN(A,B,N,X,Y,DY)
668C ENDIF
669 70 CONTINUE
670 CALL AFPOLIN(A,B,N,X,Y,DY)
671 77 AFGETFV=Y
672 RETURN
673 END
674ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
675 SUBROUTINE AFPOLIN(XA,YA,N,X,Y,DY)
676 PARAMETER (NMAX=10)
677 DIMENSION XA(NMAX),YA(NMAX),C(NMAX),D(NMAX)
678 Y=0.
679 IF(N.GT.NMAX) RETURN
680 NS=1
681 DIF=ABS(X-XA(1))
682 DO 11 I=1,N
683 DIFT=ABS(X-XA(I))
684 IF (DIFT.LT.DIF) THEN
685 NS=I
686 DIF=DIFT
687 ENDIF
688 C(I)=YA(I)
689 D(I)=YA(I)
69011 CONTINUE
691 Y=YA(NS)
692 NS=NS-1
693 DO 13 M=1,N-1
694 DO 12 I=1,N-M
695 HO=XA(I)-X
696 HP=XA(I+M)-X
697 W=C(I+1)-D(I)
698 DEN=HO-HP
699C IF(DEN.EQ.0.)PAUSE
700 DEN=W/DEN
701 D(I)=HP*DEN
702 C(I)=HO*DEN
70312 CONTINUE
704 IF (2*NS.LT.N-M)THEN
705 DY=C(NS+1)
706 ELSE
707 DY=D(NS)
708 NS=NS-1
709 ENDIF
710 Y=Y+DY
71113 CONTINUE
712 RETURN
713 END
714ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc