Removing obsolete constants.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapacfgpg.f
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
42 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
61 c 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       
80 c
81 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
90 c
91 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
92       entry ACFGPinit(Eorder,Q2fit)
93       return
94 c
95 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
96       entry ACFGPpdf(mem)
97       imem = mem
98       return
99 c
100  1000 format(5e13.5)
101       end
102 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
103       SUBROUTINE ACFGP1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
104 C
105 C     INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
106 C     QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
107 C
108 C     THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
109 C     CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
110 C
111 C     THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
112 C     BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
113 C
114 C     THE PROGRAM WORKS FOR  2. GEV**2 < Q2 <5.5E+5   AND .00137 < X < .9986
115 C
116 C     THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
117 C
118 C     THE VALUE OF LAMBDA-MSB IS 200 MEV
119 C
120 C     THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
121 C                                  X*U=X*U(X,Q2)
122 C                                  X*D= ...
123 C                                  X*S= ...
124 C                                  X*C= ...  (MASSLESS CHARM WITH     C(X,2.)=0)
125 C                                 X*CM= ...  (MASSIVE CHARM WITH MC=1.5 GEV )
126 C                              X*GLU=GLUON(X,Q2)*X
127 C
128 C
129 C                    F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
130 C                    F2C=  "        "         "     WITH MASSIVE CHARM
131 C
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/
141 C----------------------------------------------------------------------
142        DATA ISTART/0/
143        SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
144 C
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
158 C
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
179 C
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)
202 C
203       RETURN
204       END
205 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
206       SUBROUTINE ACFGP2(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
207 C
208 C     INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIVES THE
209 C     QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF X AND Q2
210 C
211 C     THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MASSIVE
212 C     CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABLE.
213 C
214 C     THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS ARE GIVEN
215 C     BY A VDM "ANSATZ" AT Q2=.25 GEV**2.
216 C
217 C     THE PROGRAM WORKS FOR  2. GEV**2 < Q2 <5.5E+5   AND .00137 < X < .9986
218 C
219 C     THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME.
220 C
221 C     THE VALUE OF LAMBDA-MSB IS 200 MEV
222 C
223 C     THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT':
224 C                                  X*U=X*U(X,Q2)
225 C                                  X*D= ...
226 C                                  X*S= ...
227 C                                  X*C= ...  (MASSLESS CHARM WITH     C(X,2.)=0)
228 C                                 X*CM= ...  (MASSIVE CHARM WITH MC=1.5 GEV )
229 C                              X*GLU=GLUON(X,Q2)*X
230 C
231 C
232 C                    F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM
233 C                    F2C=  "        "         "     WITH MASSIVE CHARM
234 C
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/
244 C----------------------------------------------------------------------
245        DATA ISTART/0/
246        SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/
247 C
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
261 C
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
282 C
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
297 C... 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)
308 C
309       RETURN
310       END
311 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
312       SUBROUTINE SFAFG1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL)
313 C
314 C**************************************************************************
315 C                       ( 1st of February 1994)
316 C     This is an interpolation program which reads the files GRPOL and
317 C     GRVDM and gives the quark and gluon distributions in real photon
318 C     as functions of x and Q**2.
319 C
320 C     The Q**2 evolution is a BLL evolution (MSbar scheme) with Nf=4
321 C     and LAMBDA(MSbar)=.200 Gev.
322 C
323 C     A massless charm distribution is generated for Q**2 > 2 Gev**2.
324 C
325 C     The distributions are the sum of a pointlike part (PL) and of a
326 C     Vdm part (VDM):
327 C                     dist=PL + KA*VDM
328 C     KA is a factor which can be adjusted ( the default value is KA=1.0).
329 C     The file GRPOL contains the pointlike part of the distributions.
330 C     The file GRVDM contains the vdm part (A precise definition of this
331 C     latter is given in the paper "PARTON DISTRIBUTIONS IN THE PHOTON",
332 C     Preprint LPTHE Orsay 93-37, by P.Aurenche,M.Fontannaz and J.Ph.Guillet).
333 C
334 C     The output of the program is written in the file GETOUT with the
335 C     following conventions
336 C                              UPLUS=x(u+ubar)
337 C                              DPLUS=x(d+dbar)
338 C                              SPLUS=x(s+sbar)
339 C                              CPLUS=x(c+cbar)
340 C                              SING =UPLUS+DPLUS+SPLUS+CPLUS
341 C                              GLU  =x*g
342 C
343 C      The interpolation is valid for     2. < Q**2 < 5.5E+5 Gev**2,
344 C                             and for   .0015<  x   < .99
345 C
346 C      The program also gives the structure function F2:
347 C                        F2 = q*Cq + g*Cg + Cgam
348 C      Cq and Cg are the Wilson coeficients and Cgam is the direct term.
349 C
350 C      Although the charm quark evolution is massless, the direct term
351 C      Cgam includes the effects due to the charm quark mass. The charm
352 C      quark threshold is therefore correctly described at the lowest
353 C      ordre in alphastrong (Details are given in the preprint).
354 C
355 C
356 C**************************************************************************
357 C
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/
370 C----------------------------------------------------------------------
371        DATA ISTART/0/
372        SAVE ISTART,OWLAM2,Q02,FLAV,KA, /W5051IA/, /W5051IB/
373 C
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
389 C
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
410 C
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
431 C
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)
455 C
456       RETURN
457       END
458 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
459       SUBROUTINE WATE32
460 C  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
511 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
541 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
552 C     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
571 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
582 C     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
601 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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
616 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
617        FUNCTION AFGETFV(X,FVL)
618 C  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
628 C      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)
664 C 70   IF(X.GT..2.AND.X.LT..8) THEN
665 C            CALL AFPOLIN(A,B,N,X,Y,DY)
666 C      ELSE
667 C            CALL AFRATIN(A,B,N,X,Y,DY)
668 C      ENDIF
669   70   CONTINUE
670              CALL AFPOLIN(A,B,N,X,Y,DY)
671   77   AFGETFV=Y
672        RETURN
673        END
674 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
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)
690 11    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
699 C         IF(DEN.EQ.0.)PAUSE
700           DEN=W/DEN
701           D(I)=HP*DEN
702           C(I)=HO*DEN
703 12      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
711 13    CONTINUE
712       RETURN
713       END
714 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc