1 SUBROUTINE CtLhALFSET (QS, ALFS)
2 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
5 DATA ALAM, BLAM, ERR / 0.01, 10.0, 0.02 /
8 CALL CtLhParQcd (2, 'ORDR', ORDR, IR1)
11 EFLLN = CtLhQZBRNT (CtLhRTALF, ALAM, BLAM, ERR, IR2)
12 EFFLAM = QS / EXP (EFLLN)
13 CALL CtLhSETL1 (NEFF, EFFLAM)
15 FUNCTION CtLhALPI (AMU)
16 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
17 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
18 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
20 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
22 IF(.NOT.SET) CALL CtLhLAMCWZ
25 CtLhALPI = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT)
27 CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI', 'AMU', AMU,
29 ELSEIF (IRT .EQ. 2) THEN
30 CALL CtLhWARNR (IW2, 'CtLhALPI > 3; Be aware!', 'CtLhALPI',
31 > CtLhALPI, D0, D1, 0)
35 FUNCTION CtLhALPQCD (IRDR, NF, RML, IRT)
36 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
37 PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15)
38 PARAMETER (CG = 3.0d0, TR = 0.5d0, CF = 4.0d0/3.0d0)
40 IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN
42 > 'Order out of range in CtLhALPQCD: IRDR = ', IRDR
45 B0 = (11.d0*CG - 2.* NF) / 3.d0
46 B1 = (34.d0*CG**2 - 10.d0*CG*NF - 6.d0*CF*NF) / 3.d0
55 IF (IRDR .GE. 2) AL = AL * (1.d0-B1*LOG(ALN) / ALN / B0**2)
62 FUNCTION CtLhAMHATF(I)
63 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
64 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
65 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
67 IF (.NOT.SET) CALL CtLhLAMCWZ
68 IF ((I.LE.0).OR.(I.GT.9)) THEN
69 print *,'warning I OUT OF RANGE IN CtLhAMHATF'
77 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
78 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
79 DATA HUGE, IWRN / 1.E20, 0 /
86 CALL CtLhWARNR(IWRN, 'CtLhDXDZ singular in CtLhDXDZ; set=HUGE',
93 SUBROUTINE CtLhEVLPAR (IACT, NAME, VALUE, IRET)
94 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
98 WRITE ( NINT(VALUE) , 101)
99 101 FORMAT (/ ' Initiation parameters: Qini, Ipd0, Ihdn ' /
100 > ' Maximum Q, Order of Alpha: Qmax, IKNL ' /
101 > ' X- mesh parameters : Xmin, Xcr, Nx ' /
102 > ' LnQ-mesh parameters : Nt, Jt ' /
103 > ' # of parton flavors : NfMx ' /)
105 ElseIF (IACT .EQ. 1) THEN
106 CALL CtLhEVLSET (NAME, VALUE, IRET)
108 print *,'fatal evlpar'
113 SUBROUTINE CtLhEVLSET (NAME, VALUE, IRET)
114 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
117 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
118 PARAMETER (MXPN = MXF * 2 + 2)
119 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
120 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
121 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
122 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
123 > / PdfSwh / Iset, IpdMod, Iptn0, NuIni
125 IF (NAME .EQ. 'QINI') THEN
126 IF (VALUE .LE. 0) GOTO 12
128 ElseIF (NAME .EQ. 'IPD0') THEN
130 IF (Item .Eq. 10 .or. Item .Eq. 11) GOTO 12
132 ElseIF (NAME .EQ. 'IHDN') THEN
134 IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12
136 ElseIF (NAME .EQ. 'QMAX') THEN
137 IF (VALUE .LE. QINI) GOTO 12
139 ElseIF (NAME .EQ. 'IKNL') THEN
142 IF (ITEM.NE.1.AND.ITEM.NE.2) GOTO 12
144 ElseIF (NAME .EQ. 'XCR') THEN
145 IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12
148 ElseIF (NAME .EQ. 'XMIN') THEN
149 IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12
152 ElseIF (NAME .EQ. 'NX') THEN
154 IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12
157 ElseIF (NAME .EQ. 'NT') THEN
159 IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12
161 ElseIF (NAME .EQ. 'JT') THEN
163 IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12
165 ElseIF (NAME .EQ. 'NFMX') THEN
167 IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12
169 ElseIF (NAME .EQ. 'IPDMOD') THEN
171 IF (Abs(Item) .Gt. 1) GOTO 12
173 ElseIF (NAME .EQ. 'IPTN0') THEN
175 IF (ABS(ITEM) .GT. MXF) GOTO 12
177 ElseIF (NAME .EQ. 'NUINI') THEN
179 IF (ITEM .LE. 0) GOTO 12
188 SUBROUTINE CtLhEVOLVE (FINI, IRET)
189 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
190 include 'parmsetup.inc'
192 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
193 PARAMETER (MXPN = MXF * 2 + 2)
194 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
195 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
196 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
197 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
198 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
199 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
200 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
201 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
202 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
204 DIMENSION JI(-MXF : MXF+1)
205 EXTERNAL LhCtNSRHSP, LhCtNSRHSM, FINI
207 save nxsave, ntsave, jtsave, ngsave,
208 & xcrsave, xminsave, qinisave, qmaxsave, ientry, ishow
210 dimension xvsave(0:MXX), qvsave(0:MXQ), tvsave(0:MXQ)
215 if(ientry .lt. 100) then
216 666 format(1x,'enter evolve',i3)
217 elseif(ientry .eq. 100) then
218 667 format(1x,'enter evolve',i3,' further suppressed')
220 ishow = 0 !set to no display
221 if(ientry .eq. 1) then
222 ishow = 1 !turn display on
225 IF (IHDN .LE. 4) THEN
227 ElseIF (IHDN .LE. 6) THEN
230 IF (.NOT. LSTX) CALL CtLhXARRAY
233 CALL CtLhPARPDF (2, 'ALAM', AL, IR)
234 CALL CtLhQARRAY (NINI)
237 Nelmt = KF * (Nt+1) * (Nx+1)
238 DO 101 IFLV = -NFMX, NFMX+1
240 JI(IFLV) = JFL * (NT+1) * (NX+1)
243 UPD(JI(0)+IZ+1,iset) = FINI (0, XV(IZ))
244 UPD(JI(NFSN)+IZ+1,iset) = 0
245 IF (NFMX .EQ. 0) GOTO 31
246 DO 331 IFLV = 1, NINI
247 A = FINI ( IFLV, XV(IZ))
248 B = FINI (-IFLV, XV(IZ))
250 UPD(JI(NFSN)+IZ+1,iset) =
251 > UPD(JI(NFSN)+IZ+1,iset) + QRKP (IFLV)
252 UPD(JI(-IFLV)+IZ+1,iset) = A - B
254 DO 332 IFLV = 1, NINI
255 UPD(JI( IFLV)+IZ+1,iset) =
256 > QRKP(IFLV) - UPD(JI(NFSN)+IZ+1,iset)/NINI
259 DO 21 NEFF = NINI, NFMX
260 IF (IKNL .EQ. 2) CALL CtLhSTUPKL (NEFF)
261 ICNT = NEFF - NINI + 1
262 IF (NTN(ICNT) .EQ. 0) GOTO 21
266 CALL CtLhSNEVL (IKNL, NX, NITR, JT, DT, TIN, NEFF,
267 > UPD(JI(NFSN)+2,iset), UPD(JI(0)+2,iset),
268 > UPD(JI(NFSN)+1,iset), UPD(JI(0)+1,iset))
269 IF (NEFF .EQ. 0) GOTO 88
270 5 DO 333 IFLV = 1, NEFF
271 CALL CtLhNSEVL (LhCtNSRHSP, IKNL, NX, NITR, JT, DT, TIN,
272 > NEFF, UPD(JI( IFLV)+2,iset), UPD(JI( IFLV)+1,iset))
274 > CALL CtLhNSEVL (LhCtNSRHSM, IKNL, NX, NITR, JT, DT, TIN,
275 > NEFF, UPD(JI(-IFLV)+2,iset), UPD(JI(-IFLV)+1,iset))
278 TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV),iset)
279 TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN),iset) / NEFF
281 IF (IKNL .GT. 0) TP = MAX (TP, DZER)
282 IF (IFLV .LE. MXVAL) THEN
283 TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV),iset)
284 IF (IKNL .GT. 0) THEN
291 UPD (JI( IFLV) + IS*(NX+1) + IX + 1,iset) = (TP + TM)/2.
292 UPD (JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = (TP - TM)/2.
296 DO 334 IFLV = NEFF + 1, NFMX
299 UPD(JI( IFLV) + IS*(NX+1) + IX + 1,iset) = 0
300 UPD(JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = 0
305 IF (NFMX .EQ. NEFF) GOTO 21
306 DO 335 IFLV = -NFMX, NFMX+1
307 JI(IFLV) = JI(IFLV) + NITR * (NX+1)
309 CALL CtLhHQRK (NX, TT, NEFF+1, UPD(JI(0)+2,iset),
310 > UPD(JI(NEFF+1)+2,iset))
312 QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1,iset)
313 UPD (JI(NFSN)+IZ+1,iset) = UPD (JI(NFSN)+IZ+1,iset)
315 VS00 = UPD (JI(NFSN)+IZ+1,iset) / (NEFF+1)
316 UPD(JI( NEFF+1) + IZ + 1,iset) = QRKP(NEFF+1) - VS00
318 A = UPD(JI( IFL)+IZ+1,iset)
319 B = UPD(JI(-IFL)+IZ+1,iset)
321 UPD(JI( IFL)+IZ+1,iset) = QRKP(IFL) - VS00
322 IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1,iset) = A - B
326 if(ientry .eq. 1) then
336 if((nx .ne. nxsave) .or.
337 & (nt .ne. ntsave) .or.
338 & (jt .ne. jtsave) .or.
339 & (ng .ne. ngsave) .or.
340 & (xcr .ne. xcrsave) .or.
341 & (xmin .ne. xminsave) .or.
342 & (qini .ne. qinisave) .or.
343 & (qmax .ne. qmaxsave)) then
344 write(6,669) nx, nt, jt, ng, xcr, xmin,
346 669 format(1x,'evolve.f: nx,nt,jt,ng=',4i4,
347 & ' xcr,xmin=',2f9.6,
348 & ' qini, qmax',f7.4,1x,e12.5,' ientry=',i6)
360 if(xvsave(i) .ne. xv(i)) then
368 if(qvsave(i) .ne. qv(i)) then
372 if(tvsave(i) .ne. tv(i)) then
379 FUNCTION CtLhFINTRP (FF, X0, DX, NX, XV, ERR, IR)
380 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
381 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
383 DIMENSION FF (0:NX), XX(MX)
384 DATA SML, XX / 1.D-5, 0., 1.0, 2.0 /
385 DATA IW1, IW3, IW5 / 3 * 0 /
392 CALL CtLhWARNI(IW1, 'Nx < 1, error in CtLhFINTRP.',
393 > 'NX', NX, 1, 256, 1)
400 CALL CtLhWARNR(IW3, 'DX < 0, error in CtLhFINTRP.',
401 > 'DX', DX, D0, D1, 1)
406 IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN
408 > 'X out of range in CtLhFINTRP, Extrapolation used.',
415 ELSEIF (TX .GE. ANX-1.) THEN
421 CALL CtLhRATINT (XX, FF(IX), MNX, DDX, TEM, ERR)
425 FUNCTION CtLhGausInt(F,XL,XR,AERR,RERR,ERR,IRT)
426 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
427 DIMENSION XLIMS(100), R(93), W(93)
428 INTEGER PTR(4),NORD(4)
430 DATA PTR,NORD/4,10,22,46, 6,12,24,48/
431 DATA R/.2386191860,.6612093865,.9324695142,
432 1 .1252334085,.3678314990,.5873179543,.7699026742,.9041172563,
433 1 .9815606342,.0640568929,.1911188675,.3150426797,.4337935076,
434 1 .5454214714,.6480936519,.7401241916,.8200019860,.8864155270,
435 1 .9382745520,.9747285560,.9951872200,.0323801710,.0970046992,
436 1 .1612223561,.2247637903,.2873624873,.3487558863,.4086864820,
437 1 .4669029048,.5231609747,.5772247261,.6288673968,.6778723796,
438 1 .7240341309,.7671590325,.8070662040,.8435882616,.8765720203,
439 1 .9058791367,.9313866907,.9529877032,.9705915925,.9841245837,
440 1 .9935301723,.9987710073,.0162767488,.0488129851,.0812974955,
441 1 .1136958501,.1459737146,.1780968824,.2100313105,.2417431561,
442 1 .2731988126,.3043649444,.3352085229,.3656968614,.3957976498,
443 1 .4254789884,.4547094222,.4834579739,.5116941772,.5393881083,
444 1 .5665104186,.5930323648,.6189258401,.6441634037,.6687183100,
445 1 .6925645366,.7156768123,.7380306437,.7596023411,.7803690438,
446 1 .8003087441,.8194003107,.8376235112,.8549590334,.8713885059,
447 1 .8868945174,.9014606353,.9150714231,.9277124567,.9393703398,
448 1 .9500327178,.9596882914,.9683268285,.9759391746,.9825172636,
449 1 .9880541263,.9925439003,.9959818430,.9983643759,.9996895039/
450 DATA W/.4679139346,.3607615730,.1713244924,
451 1 .2491470458,.2334925365,.2031674267,.1600783285,.1069393260,
452 1 .0471753364,.1279381953,.1258374563,.1216704729,.1155056681,
453 1 .1074442701,.0976186521,.0861901615,.0733464814,.0592985849,
454 1 .0442774388,.0285313886,.0123412298,.0647376968,.0644661644,
455 1 .0639242386,.0631141923,.0620394232,.0607044392,.0591148397,
456 1 .0572772921,.0551995037,.0528901894,.0503590356,.0476166585,
457 1 .0446745609,.0415450829,.0382413511,.0347772226,.0311672278,
458 1 .0274265097,.0235707608,.0196161605,.0155793157,.0114772346,
459 1 .0073275539,.0031533461,.0325506145,.0325161187,.0324471637,
460 1 .0323438226,.0322062048,.0320344562,.0318287589,.0315893308,
461 1 .0313164256,.0310103326,.0306713761,.0302999154,.0298963441,
462 1 .0294610900,.0289946142,.0284974111,.0279700076,.0274129627,
463 1 .0268268667,.0262123407,.0255700360,.0249006332,.0242048418,
464 1 .0234833991,.0227370697,.0219666444,.0211729399,.0203567972,
465 1 .0195190811,.0186606796,.0177825023,.0168854799,.0159705629,
466 1 .0150387210,.0140909418,.0131282296,.0121516047,.0111621020,
467 1 .0101607705,.0091486712,.0081268769,.0070964708,.0060585455,
468 1 .0050142027,.0039645543,.0029107318,.0018539608,.0007967921/
469 DATA TOLABS,TOLREL,NMAX/1.E-35,5.E-4,100/
477 10 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0
478 BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0
481 15 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
485 DO 20 I=PTR(J),PTR(J)-1+NORD(J)
486 20 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I)))
488 TOL=MAX(TOLABS,TOLREL*ABS(VAL))
489 IF (ABS(TVAL-VAL).LT.TOL) THEN
490 CtLhGausInt=CtLhGausInt+VAL
492 IF (NLIMS.NE.0) GO TO 10
500 IF (NLIMS.GT.(NMAX-2)) THEN
501 write(*,50) CtLhGausInt,NMAX,BB-AA,BB+AA
509 50 FORMAT (' CtLhGausInt FAILS, CtLhGausInt,NMAX,XL,XR=',
512 SUBROUTINE CtLhHINTEG (NX, F, H)
513 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
514 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
515 PARAMETER (MXPN = MXF * 2 + 2)
516 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
517 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
518 COMMON / LhCtHINTEC / GH(NDG, MXX)
519 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
520 DIMENSION F(NX), H(NX), G(MXX)
524 TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2)
527 W = XA(I,1) / XA(IY,1)
528 G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W)
530 HTEM = CtLhSMPSNA (NP-2, DZ, G(3), ERR)
532 H(I) = TEM + HTEM + TEM1
534 H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0))
538 SUBROUTINE CtLhHQRK (NX, TT, NQRK, Y, F)
539 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
540 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
541 DIMENSION Y(NX), F(NX)
542 IF (NX .GT. 1) GOTO 11
552 SUBROUTINE CtLhINTEGR (NX, M, F, G, IR)
553 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
555 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
556 PARAMETER (MXPN = MXF * 2 + 2)
557 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
558 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
559 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
560 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
561 DIMENSION F(NX), G(NX)
562 DATA IWRN1, IWRN2 / 0, 0 /
564 IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN
565 MSG = 'NX out of range in CtLhINTEGR call'
566 CALL CtLhWARNI (IWRN1, MSG, 'NX', NX, 0, MXX, 0)
569 IF (M .LT. M1 .OR. M .GT. M2) THEN
570 MSG ='Exponent M out of range in CtLhINTEGR'
571 CALL CtLhWARNI (IWRN2, MSG, 'M', M, M1, M2, 1)
575 TEM = H(1, NX-1, -M) * F(NX-2) + H(2, NX-1, -M) * F(NX-1)
576 > + H(3, NX-1, -M) * F(NX)
580 G(NX-1) = TEM * XA(NX-1, M)
582 DO 10 I = NX-2, 2, -1
583 TEM = TEM + H(1,I,-M)*F(I-1) + H(2,I,-M)*F(I)
584 > + H(3,I,-M)*F(I+1) + H(4,I,-M)*F(I+2)
588 G(I) = TEM * XA(I, M)
591 TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3)
595 G(1) = TEM * XA(1, M)
600 SUBROUTINE CtLhKERNEL
601 >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT)
602 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
603 PARAMETER (PI = 3.141592653589793d0, PI2 = PI ** 2)
604 PARAMETER (D0 = 0.0, D1 = 1.0)
605 DATA CF, CG, TR, IWRN / 1.33333333333333d0, 3.0d0, 0.5d0, 0 /
609 IF (X .LE. 0. .OR. X .GE. 1.) THEN
610 CALL CtLhWARNR(IWRN, 'X out of range in CtLhKERNEL', 'X', X,
622 SPEN2 = CtLhSPENC2 (X)
623 FFP = (1.+ X2) * XM1I
624 FGP = (2.- 2.* X + X2) / X
625 GFP = 1. - 2.* X + 2.* X2
626 GGP = XM1I + XI - 2. + X - X2
627 FFM = (1.+ X2) * XP1I
628 FGM = - (2.+ 2.* X + X2) / X
629 GFM = 1. + 2.* X + 2.* X2
630 GGM = XP1I - XI - 2. - X - X2
631 FF1 = CF * FFP * (1.- X)
634 GG1 = 2.* CG * GGP * X * (1.-X)
635 PCF2 = -2.* FFP *XLN*XLN1M - (3.*XM1I + 2.*X)*XLN
636 > - (1.+X)/2.*XLN2 - 5.*(1.-X)
637 PCFG = FFP * (XLN2 + 11.*XLN/3.+ 67./9.- PI**2 / 3.)
638 > + 2.*(1.+X) * XLN + 40.* (1.-X) / 3.
639 PCFT = (FFP * (- XLN - 5./3.) - 2.*(1.-X)) * 2./ 3.
640 PQQB = 2.* FFM * SPEN2 + 2.*(1.+X)*XLN + 4.*(1.-X)
641 PQQB = (CF**2-CF*CG/2.) * PQQB
642 PQQ2 = CF**2 * PCF2 + CF*CG * PCFG / 2. + CF*TRNF * PCFT
643 PNSP = (PQQ2 + PQQB) * (1.-X)
644 PNSM = (PQQ2 - PQQB) * (1.-X)
645 FFCF2 = - 1. + X + (1.- 3.*X) * XLN / 2. - (1.+ X) * XLN2 / 2.
646 > - FFP * (3.* XLN / 2. + 2.* XLN * XLN1M)
648 FFCFG = 14./3.* (1.-X)
649 > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.)
651 FFCFT = - 16./3. + 40./3.* X + (10.* X + 16./3.* X2 + 2.) * XLN
652 > - 112./9.* X2 + 40./9./X - 2.* (1.+ X) * XLN2
653 > - FFP * (10./9. + 2./3. * XLN)
654 FGCF2 = - 5./2.- 7./2.* X + (2.+ 7./2.* X) * XLN + (X/2.-1.)*XLN2
656 > - FGP * (3.* XLN1M + XLN1M ** 2)
657 FGCFG = 28./9. + 65./18.* X + 44./9. * X2 - (12.+ 5.*X + 8./3.*X2)
658 > * XLN + (4.+ X) * XLN2 + 2.* X * XLN1M
659 > + FGP * (-2.*XLN*XLN1M + XLN2/2. + 11./3.*XLN1M + XLN1M**2
662 FGCFT = -4./3.* X - FGP * (20./9.+ 4./3.*XLN1M)
663 GFCFT = 4.- 9.*X + (-1.+ 4.*X)*XLN + (-1.+ 2.*X)*XLN2 + 4.*XLN1M
664 > + GFP * (-4.*XLN*XLN1M + 4.*XLN + 2.*XLN2 - 4.*XLN1M
665 > + 2.*XLN1M**2 - 2./3.* PI2 + 10.)
666 GFCGT = 182./9.+ 14./9.*X + 40./9./X + (136./3.*X - 38./3.)*XLN
667 > - 4.*XLN1M - (2.+ 8.*X)*XLN2
668 > + GFP * (-XLN2 + 44./3.*XLN - 2.*XLN1M**2 + 4.*XLN1M
669 > + PI2/3. - 218./9.)
671 GGCFT = -16.+ 8.*X + 20./3.*X2 + 4./3./X + (-6.-10.*X)*XLN
672 > - 2.* (1.+ X) * XLN2
673 GGCGT = 2.- 2.*X + 26./9.*X2 - 26./9./X - 4./3.*(1.+X)*XLN
675 GGCG2 = 27./2.*(1.-X) + 67./9.*(X2-XI) + 4.*(1.+X)*XLN2
676 > + (-25.+ 11.*X - 44.*X2)/3.*XLN
677 > + GGP * (67./9.- 4.*XLN*XLN1M + XLN2 - PI2/3.)
679 FF2 = CF * TRNF * FFCFT + CF ** 2 * FFCF2 + CF * CG * FFCFG
680 FG2 = CF * TRNF * FGCFT + CF ** 2 * FGCF2 + CF * CG * FGCFG
681 GF2 = CF * TRNF * GFCFT + CG * TRNF * GFCGT
682 GG2 = CF * TRNF * GGCFT + CG ** 2 * GGCG2 + CG * TRNF * GGCGT
683 XLG = (LOG(1./(1.-X)) + 1.)
685 FF2 = FF2 * X * (1.- X)
688 GG2 = GG2 * X * (1.- X)
691 SUBROUTINE CtLhLAMCWZ
692 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
693 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
695 CALL CtLhSETL1 (NF, AL)
697 FUNCTION LhCtNAMQCD(NNAME)
698 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
699 CHARACTER NNAME*(*), NAME*8
700 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
707 IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR.
708 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') )
710 IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR.
711 1 (NAME .EQ. '# FL') )
714 IF (NAME .EQ. 'M'//CHAR(I+IASC0))
718 IF (NAME .EQ. 'LAM'//CHAR(I+IASC0))
721 IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') LhCtNAMQCD = 24
724 FUNCTION LhCtNFL(AMU)
725 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
726 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
727 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
729 IF (.NOT. SET) CALL CtLhLAMCWZ
731 IF ((LhCtNFL .EQ. NF) .OR. (AMU .LE. AMN)) GOTO 20
732 DO 10 I = NF - NHQ + 1, NF
733 IF (AMU .GE. AMHAT(I)) THEN
741 SUBROUTINE CtLhNSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN)
742 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
743 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
744 PARAMETER (MXPN = MXF * 2 + 2)
745 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
746 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
747 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
748 DIMENSION U0(NX), UN(0:NX, 0:NT)
749 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
752 IF (NX .GT. MXX) THEN
753 WRITE (*,*) 'Nx =', NX, ' greater than Max pts in CtLhNSEVL.'
754 STOP 'Program stopped in CtLhNSEVL'
756 TMD = TIN + DT * NT / 2.
758 TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
763 UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1)
770 IRND = (IS-1) * JT + JS
771 IF (IRND .EQ. 1) THEN
772 CALL RHS (TT, Neff, Y0, F0)
774 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
777 CALL RHS (TT, NEFF, Y0, F1)
779 Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
782 CALL RHS (TT, NEFF, Y1, F1)
784 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
787 CALL RHS (TT, NEFF, YP, FP)
789 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
797 UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
801 SUBROUTINE LhCtNSRHSM (TT, NEFF, FI, FO)
802 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
804 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
805 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
806 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
807 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
808 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
809 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
810 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
811 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
812 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
813 DIMENSION G1(MXX), FI(NX), FO(NX)
814 DIMENSION W0(MXX), W1(MXX), WH(MXX)
818 CPL2= CPL ** 2 / 2. * S
820 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
821 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
822 CALL CtLhHINTEG (NX, FI, WH)
824 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
825 FO(IZ) = CPL * FO(IZ)
827 IF (IKNL .EQ. 2) THEN
837 G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX))
839 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
840 TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2
841 FO(IX) = FO(IX) + TMP2
846 SUBROUTINE LhCtNSRHSP (TT, NEFF, FI, FO)
847 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
849 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
850 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
851 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
852 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
853 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
854 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
855 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
856 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
857 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
858 DIMENSION G1(MXX), FI(NX), FO(NX)
859 DIMENSION W0(MXX), W1(MXX), WH(MXX)
863 CPL2= CPL ** 2 / 2. * S
865 CALL CtLhINTEGR (NX, 0, FI, W0, IR1)
866 CALL CtLhINTEGR (NX, 1, FI, W1, IR2)
867 CALL CtLhHINTEG (NX, FI, WH)
869 FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ))
870 FO(IZ) = CPL * FO(IZ)
872 IF (IKNL .EQ. 2) THEN
879 XY = ZZ (NX-IX+1, NX-IY+1)
880 G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX))
882 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
883 TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + ZQQB)) * CPL2
884 FO(IX) = FO(IX) + TMP2
889 FUNCTION CtLhPARDIS (IPRTN, XX, QQ)
890 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
891 include 'parmsetup.inc'
894 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
895 PARAMETER (MXPN = MXF * 2 + 2)
896 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
897 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
898 PARAMETER (Smll = 1D-9)
900 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
901 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
902 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
903 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
904 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
905 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
906 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
907 COMMON / LhCtCOMQMS / VALQMS(9)
908 dimension fvec(4), fij(4)
909 dimension xvpow(0:mxx)
910 Data Iwrn1, Iwrn2, Iwrn3, OneP / 3*0, 1.00001 /
911 data xpow / 0.3d0 / !**** choice of interpolation variable
915 save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
916 save const1 , const2, const3, const4, const5, const6
917 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
921 if(nx .ne. nxsave) then
924 xvpow(i) = xv(i)**xpow
932 c enforce threshold early to improve speed...
935 if(QQ .lt. VALQMS(ii) ) then
941 c force pardis = 0.0d0 at exactly =1.0d0 - added mrw 10/May/06
942 if(xx .eq. 1.0d0) then
947 c skip the initialization in x if same as in the previous call.
948 if(x .eq. xlast) goto 100
953 11 If (JU-JLx .GT. 1) Then
955 If (X .Ge. XV(JM)) Then
962 If (JLx .LE. -1) Then
963 Print '(A,1pE12.4)','Severe error: x <= 0 in ParDis x=', x
965 ElseIf (JLx .Eq. 0) Then
967 Msg = '0 < X < Xmin in ParDis; extrapolation used!'
968 CALL CtLhWARNR (IWRN1, Msg, 'X', X, Xmin, 1D0, 1)
969 Elseif (JLx .LE. Nx-2) Then
971 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
974 Print '(A,1pE12.4)','Severe error: x > 1 in ParDis x=', x
978 If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
996 sdet = s12*s34 - s1213*s2434
998 const5 = (s34*sy2-s2434*sy3)*tmp/s12
999 const6 = (s1213*sy2-s12*sy3)*tmp/s34
1004 c skip the initialization in q if same as in the previous call.
1005 if(q .eq. qlast) goto 110
1012 12 If (JU-JLq .GT. 1) Then
1014 If (Q .GE. QV(JM)) Then
1021 If (JLq .LE. 0) Then
1023 If (JLq .LT. 0) Then
1024 Msg = 'Q < Q0 in ParDis; extrapolation used!'
1025 CALL CtLhWARNR (IWRN2, Msg, 'Q', Q, Qini, 1D0, 1)
1027 Elseif (JLq .LE. Nt-2) Then
1031 If (JLq .GE. Nt) Then
1032 Msg = 'Q > Qmax in ParDis; extrapolation used!'
1033 CALL CtLhWARNR (IWRN3, Msg, 'Q', Q, Qmax, 1D0, 1)
1037 If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
1051 tdet = t12*t34 - tmp1*tmp2
1056 jtmp = ((IPRTN + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
1058 J1 = jtmp + it*(NX+1)
1061 fij(2) = Upd(J1+1,iset) * Xa(1,2)
1062 fij(3) = Upd(J1+2,iset) * Xa(2,2)
1063 fij(4) = Upd(J1+3,iset) * Xa(3,2)
1064 Call CtLhPolint4 (XVpow(0), Fij(1), 4, ss, Fx, Dfx)
1066 If (x .GT. 0D0) Fvec(it) = Fx / x**2
1067 ElseIf (JLx .Eq. Nx-1) Then
1068 Call CtLhPolint4 (XVpow(Nx-3), Upd(J1,iset), 4, ss, Fx, Dfx)
1071 sf2 = Upd(J1+1,iset)
1072 sf3 = Upd(J1+2,iset)
1073 Fvec(it) = (const5*(Upd(J1,iset)
1074 & - sf2*const1 + sf3*const2)
1075 & + const6*(Upd(J1+3,iset)
1076 & + sf2*const3 - sf3*const4)
1077 & + sf2*sy3 - sf3*sy2) / s23
1080 If (JLq .LE. 0) Then
1081 Call CtLhPolint4 (TV(0), Fvec(1), 4, tt, ff, Dfq)
1082 ElseIf (JLq .GE. Nt-1) Then
1083 Call CtLhPolint4 (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
1087 g1 = ( tf2*t13 - tf3*t12) / t23
1088 g4 = (-tf2*t34 + tf3*t24) / t23
1089 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
1090 & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
1091 ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
1097 SUBROUTINE CtLhPARPDF (IACT, NAME, VALUE, IRET)
1098 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1099 CHARACTER NAME*(*), Uname*10
1101 DATA ILEVEL, LRET / 1, 1 /
1103 CALL CtLhUPC (NAME, Ln, Uname)
1104 IF (IACT .EQ. 0 .OR. IACT .EQ. 4) then
1105 c > IVALUE = NINT (VALUE) !tentatively remove this since it seems not to be used
1106 print *,'Fatal error: iact=',iact
1109 START1 = (IACT .NE. 1) .AND. (IACT .NE. 2)
1110 c prepare to remove this stuff, since I think IACT=1 or 2 always
1112 print *,'Fatal error: start1=',start1
1115 IF (START1) ILEVEL = 1
1119 CALL CtLhParQcd (IACT, Uname(1:Ln), VALUE, JRET)
1120 IF (JRET .EQ. 1) GOTO 11
1121 IF (JRET .EQ. 2) GOTO 12
1122 IF (JRET .EQ. 3) GOTO 13
1123 IF (JRET .GT. 4) GOTO 15
1125 2 CALL CtLhEVLPAR (IACT, Uname(1:Ln), VALUE, JRET)
1126 IF (JRET .EQ. 1) GOTO 11
1127 IF (JRET .EQ. 2) GOTO 12
1128 IF (JRET .EQ. 3) GOTO 13
1129 IF (JRET .GT. 4) GOTO 15
1131 IF (.NOT. START1) GOTO 1
1132 IF (JRET .EQ. 0) GOTO 10
1141 IF (JRET .NE. 4) LRET = JRET
1142 IF (LRET.EQ.0 .OR. LRET.EQ.2 .OR. LRET.EQ.3) THEN
1143 PRINT *, 'Error in CtLhPARPDF: IRET, IACT, NAME, VALUE =',
1144 > LRET, IACT, NAME, VALUE
1145 PRINT *, 'fatal error in CtLhparpdf'
1152 SUBROUTINE CtLhParQcd(IACT,NAME,VALUE,IRET)
1153 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1158 WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ',
1159 > '(i in 1 to 9), LAMi (i in 1 to NFL)'
1161 ELSEIF (IACT.EQ.1) THEN
1162 CALL CtLhQCDSET (NAME,VALUE,IRET)
1163 ELSEIF (IACT.EQ.2) THEN
1164 CALL CtLhQCDGET (NAME,VALUE,IRET)
1170 FUNCTION CtLhPFF1 (XX)
1171 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1172 LOGICAL LA, LB, LSTX
1173 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1174 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1175 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1177 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1178 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1179 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1180 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1181 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1182 > PNSP(0:MXX), PNSM(0:MXX)
1184 DATA LA, LB / 2 * .FALSE. /
1188 2 IF (LA .AND. .NOT.LB) THEN
1197 ElseIF (X .GE. XMIN) THEN
1199 TEM = CtLhFINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT)
1201 CALL CtLhPOLIN1 (XL, FF1(1), MX, X, TEM, ERR)
1205 CtLhPFF1 = TEM / (1.-X)
1208 CtLhTFF1 = TEM / (1.-X) * CtLhDXDZ(Z)
1216 RFF1 = TEM * X / (1.-X)
1225 ElseIF (X .GE. XMIN) THEN
1227 TEM = CtLhFINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT)
1229 CALL CtLhPOLIN1 (XL, PNSP(1), MX, X, TEM, ERR)
1231 CtLhFNSP = TEM / (1.- X)
1238 ElseIF (X .GE. XMIN) THEN
1240 TEM = CtLhFINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT)
1242 CALL CtLhPOLIN1 (XL, PNSM(1), MX, X, TEM, ERR)
1244 CtLhFNSM = TEM / (1.- X)
1251 ElseIF (X .GE. XMIN) THEN
1253 TEM = CtLhFINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT)
1255 CALL CtLhPOLIN1 (XL, GG1(1), MX, X, TEM, ERR)
1258 PGG1 = TEM / X / (1.-X)
1265 CtLhRGG1 = TEM / (1.-X)
1274 ElseIF (X .GE. XMIN) THEN
1276 TEM = CtLhFINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT)
1278 CALL CtLhPOLIN1 (XL, FF2(1), MX, X, TEM, ERR)
1281 PFF2 = TEM / X / (1.-X)
1288 CtLhRFF2 = TEM / (1.-X)
1297 ElseIF (X .GE. XMIN) THEN
1299 TEM = CtLhFINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT)
1301 CALL CtLhPOLIN1 (XL, GG2(1), MX, X, TEM, ERR)
1304 PGG2 = TEM / X / (1.-X)
1311 CtLhRGG2 = TEM / (1.-X)
1316 SUBROUTINE CtLhPOLIN1 (XA,YA,N,X,Y,DY)
1317 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1319 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1324 IF (DIFT.LT.DIF) THEN
1343 IF (2*NS.LT.N-M)THEN
1353 SUBROUTINE CtLhQARRAY (NINI)
1354 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1355 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1356 PARAMETER (MXPN = MXF * 2 + 2)
1357 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
1358 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1359 COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF)
1360 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1362 IF (NT .GE. mxq) NT = mxq - 1
1367 1 DT0 = (TMAX - TINI) / float(NT)
1368 NINI = LhCtNFL(QINI)
1369 NFMX = LhCtNFL(QMAX)
1370 Call CtLhParQcd (2, 'ORDER', Ord, Ir)
1371 Call CtLhParQcd (2, 'ALAM', Al0, Ir)
1372 Call CtLhParQcd (2, 'NFL', Afl0, Ir)
1374 Call CtLhParQcd (1, 'NFL', AFL, Ir)
1377 Call CtLhSetLam (Ifl0, Al0, Iordr)
1378 NG = NFMX - NINI + 1
1387 DO 20 NEFF = NINI, NFMX
1388 ICNT = NEFF - NINI + 1
1389 IF (NEFF .LT. NFMX) THEN
1390 THRN = CtLhAMHATF (NEFF + 1)
1391 QOUN = MIN (QMAX, THRN)
1395 IF (QOUN-QOUT .LE. 0.0001) THEN
1403 NITR = INT (TEM / DT0) + 1
1409 NTL (ICNT+1) = NTL(ICNT) + NITR
1410 IF (NITR .NE. 0) THEN
1412 TV (NTL(ICNT)+I) = TIN + DT * I
1413 S = EXP (TV(NTL(ICNT)+I))
1414 QV (NTL(ICNT)+I) = AL * EXP (S)
1423 IF (NTP .GE. MXQ) THEN
1424 NT = MXQ - ND - NCNT
1430 SUBROUTINE CtLhQCDGET(NAME,VALUE,IRET)
1431 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1433 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1434 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1435 COMMON / LhCtCOMQMS / VALQMS(9)
1437 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1438 ICODE = LhCtNAMQCD(NAME)
1440 IF (ICODE .EQ. 1) THEN
1442 ELSEIF (ICODE .EQ. 2) THEN
1444 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 12)) THEN
1445 VALUE = VALQMS(ICODE - 2)
1446 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1447 VALUE = ALAM(ICODE - 13)
1448 ELSEIF (ICODE .EQ. 24) THEN
1454 SUBROUTINE CtLhQCDSET (NAME,VALUE,IRET)
1455 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1457 COMMON / LhCtCOMQMS / VALQMS(9)
1458 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1460 PARAMETER (PI=3.141592653589793d0, EULER=0.57721566)
1461 IVALUE = NINT(VALUE)
1462 ICODE = LhCtNAMQCD(NAME)
1463 IF (ICODE .EQ. 0) THEN
1465 c print *,'warning empty CtLhQCDSET call: NAME=',
1466 c & NAME,' VALUE=',VALUE
1470 IF (ICODE .EQ. 1) THEN
1471 IF (VALUE.LE.0) GOTO 12
1473 ELSEIF (ICODE .EQ. 2) THEN
1474 IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12
1476 ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 11)) THEN
1477 IF (VALUE .LT. 0) GOTO 12
1478 Scle = Min (Value , VALQMS(ICODE - 2))
1479 AlfScle = CtLhALPI(Scle) * Pi
1480 VALQMS(ICODE - 2) = VALUE
1481 Call CtLhAlfSet (Scle, AlfScle)
1482 ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN
1483 IF (VALUE .LE. 0) GOTO 12
1484 CALL CtLhSETL1 (ICODE-13, VALUE)
1485 ELSEIF (ICODE .EQ. 24) THEN
1486 IF ((IVALUE .LT. 1) .OR. (IVALUE .GT. 2)) GOTO 12
1489 IF (.NOT. SET) CALL CtLhLAMCWZ
1495 FUNCTION CtLhQZBRNT(FUNC, X1, X2, TOLIN, IRT)
1496 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1497 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
1504 IF(FB*FA.GT.0.) THEN
1505 WRITE (*, *) 'Root must be bracketed for CtLhQZBRNT.'
1510 IF(FB*FC.GT.0.) THEN
1516 IF(ABS(FC).LT.ABS(FB)) THEN
1524 TOL1=2.*EPS*ABS(B)+0.5*TOL
1526 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
1530 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
1538 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
1539 Q=(Q-1.)*(R-1.)*(S-1.)
1543 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
1556 IF(ABS(D) .GT. TOL1) THEN
1563 WRITE (*, *) 'CtLhQZBRNT exceeding maximum iterations.'
1568 SUBROUTINE CtLhRATINT(XA,YA,N,X,Y,DY)
1569 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1570 PARAMETER (NMAX=10,TINY=1.E-25)
1571 DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
1580 ELSE IF (H.LT.HH) THEN
1599 IF (2*NS.LT.N-M)THEN
1609 FUNCTION CtLhRTALF (EFLLN)
1610 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1611 include 'parmsetup.inc'
1612 PARAMETER (PI = 3.141592653589793d0)
1613 COMMON / CtLhRTALFC / ALFST, JORD, NEFF
1614 EFMULM = EXP (EFLLN)
1616 TEM2 = 1. / CtLhALPQCD (JORD, NEFF, EFMULM, I)
1617 CtLhRTALF = TEM1 - TEM2
1619 Subroutine CtLhbldat1
1620 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1621 include 'parmsetup.inc'
1623 PARAMETER (MXX = 105, MXQ = 25, MxF = 6)
1624 PARAMETER (MxPN = MxF * 2 + 2)
1625 PARAMETER (MxQX= MXQ * MXX, MxPQX = MxQX * MxPN)
1626 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1627 COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG
1628 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1629 COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt
1630 PARAMETER (NF0 = 4, Nshp = 8,NEX = Nshp+2)
1636 SUBROUTINE CtLhSETL1 (NEF, VLAM)
1637 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1639 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
1640 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1641 COMMON / LhCtCOMQMS / VALQMS(9)
1642 IF (NEF .LT. 0 .OR. NEF .GT. NF) THEN
1643 WRITE(*,*)'NEF out of range in CtLhSETL1: NEF NF =',NEF,NF
1648 AMHAT(N) = VALQMS(N)
1651 DO 10 N = NEF, 1, -1
1652 CALL CtLhTRNLAM(NORDER, N, -1, IR1)
1655 CALL CtLhTRNLAM(NORDER, N, 1, IR1)
1657 DO 30, N = NF, 1, -1
1658 IF ((ALAM(N) .GE. 0.7 * AMHAT(N))
1659 > .OR. (ALAM(N-1) .GE. 0.7 * AMHAT(N)))THEN
1666 DO 50, N = NF-NHQ, 1, -1
1672 IF (ALAM(N) .GT. AMN) AMN = ALAM(N)
1679 SUBROUTINE CtLhSETLAM (NEF, WLAM, IRDR)
1680 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1681 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1683 IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN
1684 WRITE(*,*)'NEF out of range in CtLhSETLAM: NEF NF=',NEF,NF
1688 IF (IRDR .NE. NORDER) then
1689 PRINT *,'fatal error: wanted cnvl1'
1692 CALL CtLhSETL1 (NEF, VLAM)
1694 Subroutine CtLhbldat2
1695 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1696 COMMON / LhCtCOMQMS / VALQMS(9)
1697 COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET
1714 FUNCTION CtLhSMPNOL (NX, DX, FN, ERR)
1715 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1718 IF (NX .LE. 1 .OR. NX .GT. 1000) THEN
1719 PRINT *, 'NX =', NX, ' OUT OF RANGE IN CtLhSMPNOL!'
1721 ELSEIF (NX .EQ. 2) THEN
1723 ELSEIF (NX .EQ. 3) THEN
1724 TEM = DX * FN(2) * 2.
1727 TEM = DX * (23.* FN(2) - 16.* FN(3) + 5.* FN(4)) / 12.
1728 TMP = DX * (3.* FN(2) - FN(3)) / 2.
1729 ERR = ABS(TEM - TMP)
1730 TEM = TEM + CtLhSMPSNA (NX-1, DX, FN(2), ER1)
1731 ERR = ABS(ER1) + ERR
1733 TEM = DX * (8.* FN(2) - 4.* FN(3) + 8.* FN(4)) / 3.
1734 TMP = DX * (3.* FN(2) + 2.* FN(3) + 3.* FN(4)) / 2.
1735 ERR = ABS(TEM - TMP)
1736 TEM = TEM + CtLhSMPSNA (NX-4, DX, FN(5), ER1)
1737 ERR = ABS(ER1) + ERR
1743 FUNCTION CtLhSMPSNA (NX, DX, F, ERR)
1744 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1745 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1746 PARAMETER (MAXX = 1000)
1748 DATA IW1, IW2, TINY / 2*0, 1.E-35 /
1749 IF (DX .LE. 0.) THEN
1750 CALL CtLhWARNR(IW2,'DX cannot be < 0. in CtLhSMPSNA', 'DX',
1755 IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN
1756 CALL CtLhWARNI(IW1, 'NX out of range in CtLhSMPSNA', 'NX', NX,
1759 ELSEIF (NX .EQ. 1) THEN
1761 ELSEIF (NX .EQ. 2) THEN
1762 SIMP = (F(1) + F(2)) / 2.
1763 ERRD = (F(1) - F(2)) / 2.
1767 ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24.
1774 SIMP = (F(1) + 4.* F(2) + F(3)) / 3.
1775 TRPZ = (F(1) + 2.* F(2) + F(3)) / 2.
1785 SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3.
1786 TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2.
1791 CtLhSMPSNA = SIMP * DX
1792 IF (ABS(SIMP) .GT. TINY) THEN
1799 SUBROUTINE CtLhSNEVL(IKNL,NX,NT,JT,DT,TIN,NEFF,UI,GI,US,GS)
1800 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1801 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1802 PARAMETER (MXQX= MXQ * MXX)
1803 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1804 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1805 DIMENSION UI(NX), US(0:NX, 0:NT)
1806 DIMENSION GI(NX), GS(0:NX, 0:NT)
1807 DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX)
1808 DIMENSION Z0(MXX), Z1(MXX), ZP(MXX), G0(MXX), G1(MXX), GP(MXX)
1812 IF (NX .GT. MXX) THEN
1813 WRITE (*,*) 'Nx =', NX, ' too many pts in CtLhSNEVL'
1814 STOP 'Program stopped in CtLhSNEVL'
1816 TMD = TIN + DT * NT / 2.
1818 TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU)
1824 US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3)
1825 GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3)
1833 IRND = (IS-1) * JTT + JS
1834 IF (IRND .EQ. 1) THEN
1835 CALL CtLhSNRHS (TT, NEFF, Y0,Z0, F0,G0)
1837 Y0(IZ) = Y0(IZ) + DDT * F0(IZ)
1838 Z0(IZ) = Z0(IZ) + DDT * G0(IZ)
1841 CALL CtLhSNRHS (TT, NEFF, Y0, Z0, F1, G1)
1843 Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0
1844 Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0
1847 CALL CtLhSNRHS (TT, NEFF, Y1, Z1, F1, G1)
1849 YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0
1850 ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0
1853 CALL CtLhSNRHS (TT, NEFF, YP, ZP, FP, GP)
1855 Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0
1856 Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0
1863 IF (IKNL .GT. 0) THEN
1864 US (IX, IS) = MAX(Y1(IX), D0)
1865 GS (IX, IS) = MAX(Z1(IX), D0)
1867 US (IX, IS) = Y1(IX)
1868 GS (IX, IS) = Z1(IX)
1871 US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3)
1872 GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3)
1876 SUBROUTINE CtLhSNRHS (TT, NEFF, FI, GI, FO, GO)
1877 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1879 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1880 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1881 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1882 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1883 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1884 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1885 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1886 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1887 COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx
1888 DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX)
1889 DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), WM(MXX)
1890 DIMENSION R0(MXX), R1(MXX), R2(MXX), RH(MXX), RM(MXX)
1894 CPL2= CPL ** 2 / 2. * S
1896 CALL CtLhINTEGR (NX,-1, FI, WM, IR1)
1897 CALL CtLhINTEGR (NX, 0, FI, W0, IR2)
1898 CALL CtLhINTEGR (NX, 1, FI, W1, IR3)
1899 CALL CtLhINTEGR (NX,-1, GI, RM, IR4)
1900 CALL CtLhINTEGR (NX, 0, GI, R0, IR5)
1901 CALL CtLhINTEGR (NX, 1, GI, R1, IR6)
1902 CALL CtLhINTEGR (NX, 2, GI, R2, IR7)
1903 CALL CtLhHINTEG (NX, FI, WH)
1904 CALL CtLhHINTEG (NX, GI, RH)
1905 IF (IKNL .GT. 0) THEN
1907 FO(IZ) = ( 2D0 * FI(IZ)
1908 > + 4D0 / 3D0 * ( 2D0 * WH(IZ) - W0(IZ) - W1(IZ) ))
1909 > + NEFF * ( R0(IZ) - 2D0 * R1(IZ) + 2D0 * R2(IZ) )
1910 FO(IZ) = FO(IZ) * CPL
1911 GO(IZ) = 4D0 / 3D0 * ( 2D0 * WM(IZ) - 2D0 * W0(IZ) + W1(IZ) )
1912 > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ)
1913 > + 6D0 * (RH(IZ) + RM(IZ) - 2D0 * R0(IZ) + R1(IZ) - R2(IZ))
1914 GO(IZ) = GO(IZ) * CPL
1918 FO(IZ) = NEFF * (-R0(IZ) + 2.* R1(IZ) )
1919 > + 2.* FI(IZ) + 4./ 3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ) )
1920 FO(IZ) = FO(IZ) * CPL
1921 GO(IZ) = 4./ 3.* ( 2.* W0(IZ) - W1(IZ) )
1922 >+ (33.- 2.* NEFF) / 6.* GI(IZ) + 6.*(RH(IZ) + R0(IZ) - 2.* R1(IZ))
1923 GO(IZ) = GO(IZ) * CPL
1926 IF (IKNL .EQ. 2) THEN
1938 G1(KZ) = FFG(I, IY) * (FI(IY) - XY**2 *FI(I))
1939 G4(KZ) = GGF(I, IY) * (GI(IY) - XY**2 *GI(I))
1940 G2(KZ) = FFG(IS,IT) * (GI(IY) - xy*GI(I)) !FG
1941 G3(KZ) = GGF(IS,IT) * (FI(IY) - XY*FI(I)) !GF (usual notations)
1943 TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR)
1944 TEM2 = CtLhSMPSNA (NP, DZ, G2, ERR)
1945 TEM3 = CtLhSMPSNA (NP, DZ, G3, ERR)
1946 TEM4 = CtLhSMPNOL (NP, DZ, G4, ERR)
1947 TEM1 = TEM1 - FI(I) * (AFF2(I) + ZGF2)
1948 TEM4 = TEM4 - GI(I) * (AGG2(I) + ZFG2)
1949 tem2 = tem2 + GI(I)*AFG2(I)
1950 tem3= tem3 + FI(I)*AGF2(I)
1953 FO(I) = FO(I) + TMF * CPL2
1954 GO(I) = GO(I) + TMG * CPL2
1959 FUNCTION CtLhSPENC2 (X)
1960 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1962 COMMON / LhCtSPENCC / XX
1963 DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 /
1965 TEM = CtLhGausInt(CtLhSPN2IN, XX, U1, AERR, RERR, ERR, IRT)
1966 CtLhSPENC2 = TEM + LOG (XX) ** 2 / 2.
1969 FUNCTION CtLhSPN2IN (ZZ)
1970 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1971 COMMON / LhCtSPENCC / X
1973 TEM = LOG (1.+ X - Z) / Z
1977 SUBROUTINE CtLhSTUPKL (NFL)
1978 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1980 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
1981 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
1983 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
1984 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
1985 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
1986 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
1987 COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX),
1988 > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX),
1989 > PNSP(0:MXX), PNSM(0:MXX)
1990 COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX)
1991 COMMON / LhCtKRNL00 / DZ, XL(MX), NNX
1992 COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX),
1993 > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB
1994 EXTERNAL CtLhPFF1, CtLhRGG1, CtLhRFF2, CtLhRGG2
1995 EXTERNAL CtLhFNSP, CtLhFNSM
1996 dimension aff1(mxx),agg1(mxx)
1997 PARAMETER (PI = 3.141592653589793d0, PI2 = PI**2)
1998 DATA CF, CG, TR / 1.33333333333333d0, 3.0, 0.5 /
1999 data zeta3/1.20205690315959d0/ ! zeta(3.0)
2001 DATA AERR, RERR / 0.0, 0.02 /
2009 CALL CtLhKERNEL (XZ, FF1(I), GF1(I), FG1(I), GG1(I), PNSP(I),
2010 > PNSM(I), FF2(I), GF2(I), FG2(I), GG2(I), NFL, IRT)
2012 FF1(0) = FF1(1) * 3. - FF1(2) * 3. + FF1(3)
2013 FG1(0) = FG1(1) * 3. - FG1(2) * 3. + FG1(3)
2014 GF1(0) = GF1(1) * 3. - GF1(2) * 3. + GF1(3)
2015 GG1(0) = GG1(1) * 3. - GG1(2) * 3. + GG1(3)
2016 PNSP(0) = PNSP(1) * 3. - PNSP(2) * 3. + PNSP(3)
2017 PNSM(0) = PNSM(1) * 3. - PNSM(2) * 3. + PNSM(3)
2018 FF2(0) = FF2(1) * 3. - FF2(2) * 3. + FF2(3)
2019 FG2(0) = FG2(1) * 3. - FG2(2) * 3. + FG2(3)
2020 GF2(0) = GF2(1) * 3. - GF2(2) * 3. + GF2(3)
2021 GG2(0) = GG2(1) * 3. - GG2(2) * 3. + GG2(3)
2022 FF1(NX) = FF1(NX-1) * 3. - FF1(NX-2) * 3. + FF1(NX-3)
2023 FG1(NX) = FG1(NX-1) * 3. - FG1(NX-2) * 3. + FG1(NX-3)
2024 GF1(NX) = GF1(NX-1) * 3. - GF1(NX-2) * 3. + GF1(NX-3)
2025 GG1(NX) = GG1(NX-1) * 3. - GG1(NX-2) * 3. + GG1(NX-3)
2026 PNSM(NX) = PNSM(NX-1) * 3. - PNSM(NX-2) * 3. + PNSM(NX-3)
2027 PNSP(NX) = PNSP(NX-1) * 3. - PNSP(NX-2) * 3. + PNSP(NX-3)
2028 FF2(NX) = FF2(NX-1) * 3. - FF2(NX-2) * 3. + FF2(NX-3)
2029 FG2(NX) = FG2(NX-1) * 3. - FG2(NX-2) * 3. + FG2(NX-3)
2030 GF2(NX) = GF2(NX-1) * 3. - GF2(NX-2) * 3. + GF2(NX-3)
2031 GG2(NX) = GG2(NX-1) * 3. - GG2(NX-2) * 3. + GG2(NX-3)
2033 AFF1(1) = CtLhGausInt(CtLhPFF1,D0,XV(1),AERR,RERR,ER1,IRT)
2035 TMPG = CtLhGausInt(CtLhRGG1,D0,XV(1),AERR,RERR,ER3,IRT)
2036 AGG1(1) = TMPG + DGG1
2037 ANSM(1) = CtLhGausInt(CtLhFNSM,D0,XV(1),AERR,RER,ER2,IRT)
2038 ANSP(1) = CtLhGausInt(CtLhFNSP,D0,XV(1),AERR,RER,ER2,IRT)
2040 AFF2(1) = CtLhGausInt(CtLhRFF2, D0, XV(1), AER, RER, ER2, IRT)
2042 AGG2(1) = CtLhGausInt(CtLhRGG2, D0, XV(1), AER, RER, ER4, IRT)
2044 TEM =CtLhGausInt(CtLhPFF1,XV(I2-1),XV(I2),AERR,RERR,ER1,IRT)
2045 AFF1(I2) = TEM + AFF1(I2-1)
2046 AER = ABS(TEM * RER)
2047 AFF2(I2)=CtLhGausInt(CtLhRFF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT)
2049 TEM = CtLhGausInt(CtLhRGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT)
2051 AGG1(I2) = TMPG + DGG1
2052 AER = ABS(TEM * RER)
2053 AGG2(I2)=CtLhGausInt(CtLhRGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT)
2055 ANSP(I2)=CtLhGausInt(CtLhFNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2057 ANSM(I2)=CtLhGausInt(CtLhFNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT)
2060 ANSP(NX)=CtLhGausInt(CtLhFNSP,XV(NX-1),D1,AERR,RER,ERR,
2062 ANSM(NX)=CtLhGausInt(CtLhFNSM,XV(NX-1),D1,AERR,RER,ERR,
2065 do i2=1,nx-1 !loop over x
2067 XI = 1./ X !auxiliary definitions
2072 XLN1M = DLOG (1.- X)
2076 xLi31m=CtLhxLi(3,1d0-x)
2085 > (9 + 4*Pi2 - 22*x + 13*x2 + 6*(3 - 4*x + x2)*xln1m +
2086 > 40*xln - 24*xLi2)/9.
2088 > (6*(2*(-9 + Pi2) + 3*x*(5 + x)) +4*(3 +2*Pi2+3*x*(-3 + 2*x))*
2089 > xln1m + 6*x3m*x1m*xln1m2 - 6*(x*(8 + 3*x) + 4*xln1m2)*
2090 > xln - 3*(-4 + x)*x*xln2)/12 - 2*(3 + 2*xln1m)*xLi2 - 4*xLi31m
2092 > (3637-186*Pi2-x*(3198+72*Pi2+x*(231 + 208*x)))/108.- xln +
2093 > (3*xln1m*(-33 - 4*Pi2 + (50 - 17*x)*x - 3*x3m*x1m*xln1m) +
2094 > 2*(x*(198 + x*(27+8*x))+9*xln1m*(3 - 4*x + x2 + 2*xln1m))*
2095 > xln - 9*x*(4 + x)*xln2)/18- x1p*x3p*xln*xln1p-
2096 > (x1p*x3p - 4*xln)*xLi2m + (31d0/3d0 +4*xln1m- 4*xln)*xLi2 +
2097 > 4*xLi31m + 12*xLi3 - 2*xLi32 - 10*zeta3
2099 > (18 - 81*x + 6*Pi2*x + 123*x2 - 6*Pi2*x2 - 60*x3 +
2100 > 4*Pi2*x3 - 6*(-2 + 3*x - 3*x2 + 2*x3)*xln1m2 -33*x*xln +
2101 > 15*x2*xln - 24*x3*xln - 9*x*xln2 + 9*x2*xln2 -
2102 > 12*x3*xln2 - 12*x1m*xln1m*(-1 + 2*x2 + 2*xln - x*xln +
2103 > 2*x2*xln) - 24*xLi2)/9.
2105 > (2*(-67 + 2*Pi2 + x*(64 + x*(-91 + 3*Pi2 + 94*x)) +
2106 > x1m*(7+x*(-5+16*x))*xln1m -3*x1m*(2+ x*(-1+2*x))*xln1m2 -
2107 > 20*xln - 3*x*xln*(13 + 16*x*x1p - 3*x1p*xln) +
2108 > 6*x1p*(2+x+2*x2)*xln*xln1p+6*x1p*(2+x+2*x2)*xLi2m))/9.
2109 AGF2(I2) = CF*TRNF*WGFCFT+CF**2* WGFCF2+CF*CG*WGFCFG
2110 AFG2(I2) = CF*TRNF*WFGCFT +CG*TRNF*WFGCGT
2114 ZGF2=-28./27.*Cf**2+94./27.*Cf*Cg -52./27.*Cf*TrNf
2115 ZFG2= 37./27.*Cf*TrNf + 35./54.*Cg*TrNf
2116 ZQQB=1.43862321154902*(Cf**2-0.5*Cf*Cg)
2121 XG2 = (LOG(1./(1.-X)) + 1.) ** 2
2122 FFG (IS, IS) = FG2(NX) * DXTZ(I) * XG2
2123 GGF (IS, IS) = GF2(NX) * DXTZ(I) * XG2
2124 PNS (IS, IS) =PNSM(NX) * DXTZ(I)
2130 XG2 = (LOG(1./XM1) + 1.) ** 2
2137 TEM = (FF2(IZ) * (1.- DT) + FF2(IZ+1) * DT) / XM1 / XY
2138 FFG (IX, IY) = TEM * DXTZ(IY)
2139 TEM = (FG2(IZ) * (1.- DT) + FG2(IZ+1) * DT) * XG2 / XY
2140 FFG (IS, IT) = TEM * DXTZ(IY)
2141 TEM = (GF2(IZ) * (1.- DT) + GF2(IZ+1) * DT) * XG2 / XY
2142 GGF (IS, IT) = TEM * DXTZ(IY)
2143 TEM = (GG2(IZ) * (1.- DT) + GG2(IZ+1) * DT) / XM1 / XY
2144 GGF (IX, IY) = TEM * DXTZ(IY)
2145 TEM = (PNSP(IZ) * (1.- DT) + PNSP(IZ+1) * DT) / XM1
2146 PNS (IX, IY) = TEM * DXTZ(IY)
2147 TEM = (PNSM(IZ) * (1.- DT) + PNSM(IZ+1) * DT) / XM1
2148 PNS (IS, IT) = TEM * DXTZ(IY)
2153 SUBROUTINE CtLhTRNLAM (IRDR, NF, IACT, IRT)
2154 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2155 COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ
2156 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2158 DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 /
2159 DATA IR1, SML / 0, 1.E-5 /
2165 IF (JACT .GT. 0) THEN
2168 ALM = LOG (THMS/VLAM)
2174 THMS = MAX (THMS, SML)
2175 BLM = LOG (THMS/VLAM)
2177 IF (VLAM .GE. 0.7 * THMS) THEN
2178 IF (JACT . EQ. 1) THEN
2187 IF (ALM .GE. BLM) THEN
2188 WRITE (*, *) 'CtLhTRNLAM has ALM >= BLM: ', ALM, BLM
2189 WRITE (*, *) 'I do not know how to continue'
2193 ERR = RERR * LOG (VMULM)
2194 WLLN = CtLhQZBRNT (CtLhZBRLAM, ALM, BLM, ERR, IR1)
2195 ALAM(N1) = THMS / EXP (WLLN)
2196 IF (IR1 .NE. 0) THEN
2197 WRITE (*, *) 'CtLhQZBRNT failed in CtLhTRNLAM; ',
2198 > 'NF, VLAM =', NF, VLAM
2199 WRITE (*, *) 'I do not know how to continue'
2204 SUBROUTINE CtLhUPC (A, La, UpA)
2205 CHARACTER A*(*), UpA*(*), C*(1)
2209 If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!'
2210 Ld = ICHAR('A')-ICHAR('a')
2214 IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN
2215 UpA (I:I) = CHAR(Ichar(c) + ld)
2226 SUBROUTINE CtLhWARNI (IWRN, MSG, NMVAR, IVAB,
2228 CHARACTER*(*) MSG, NMVAR
2235 PRINT '(1X,A/1X, 2A,I10 /A,I4)', MSG, NMVAR, ' = ', IV
2236 IF (IACT .EQ. 1) THEN
2237 PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX
2240 If (Iw .LT. Nmax) Then
2241 PRINT '(1X,A/1X,I10,A, I10)', MSG, NMVAR, ' = ', IV
2242 Elseif (Iw .Eq. Nmax) Then
2243 Print '(/A/)', 'CtLhWARNI Severe Warning: Too many errors'
2248 SUBROUTINE CtLhWARNR (IWRN, MSG, NMVAR, VARIAB,
2250 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2251 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2252 CHARACTER*(*) MSG, NMVAR
2258 PRINT '(1X, A/1X,2A,1PD16.7/A,I4)', MSG, NMVAR, ' = ', VR
2259 IF (IACT .EQ. 1) THEN
2260 PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX
2263 If (Iw .LT. Nmax) Then
2264 PRINT '(I5, 2A/1X,2A,I10,1PD16.7)', IW, ' ', MSG,
2266 Elseif (Iw .Eq. Nmax) Then
2267 Print '(/A/)', 'CtLhWARNR Severe Warning: Too many errors'
2272 SUBROUTINE CtLhXARRAY
2273 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2275 PARAMETER (D0 = 0.0, D10=10.0)
2276 PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
2277 PARAMETER (MXPN = MXF * 2 + 2)
2278 PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN)
2279 PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2)
2281 COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX)
2282 COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2)
2283 COMMON / LhCtHINTEC / GH(NDG, MXX)
2284 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2285 COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX)
2286 DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG)
2287 DATA F12, F22, F32 / 1D0, 1D0, 1D0 /
2288 DATA (G1(I,NDH), G2(I,1), I=1,NDG) / 0.0,0.0,0.0,0.0,0.0,0.0 /
2296 DXTZ(I) = CtLhDXDZ(Z) / X
2301 IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L
2307 DXTZ(NX) = CtLhDXDZ(1.D0)
2313 ELY(I) = LOG(1D0 - XV(I))
2315 ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3)
2319 XY = XV(IX) / XV(IY)
2320 ZZ (IX, IY) = CtLhZFRMX (XY)
2321 ZZ (NX-IX+1, NX-IY+1) = XY
2324 IF (I .NE. NX-1) THEN
2331 DET = F11*F22*F33 + F21*F32*F13 + F31*F12*F23
2332 > - F31*F22*F13 - F21*F12*F33 - F11*F32*F23
2333 IF (ABS(DET) .LT. PUNY) THEN
2334 Msg='Determinant close to zero; will be arbitrarily set to:'
2335 CALL CtLhWARNR(IWRN, Msg, 'DET', PUNY, D0, D0, 0)
2338 G2(1,2) = (F22*F33 - F23*F32) / DET
2339 G2(1,3) = (F32*F13 - F33*F12) / DET
2340 G2(1,4) = (F12*F23 - F13*F22) / DET
2341 G2(2,2) = (F23*F31 - F21*F33) / DET
2342 G2(2,3) = (F33*F11 - F31*F13) / DET
2343 G2(2,4) = (F13*F21 - F11*F23) / DET
2344 G2(3,2) = (F21*F32 - F22*F31) / DET
2345 G2(3,3) = (F31*F12 - F32*F11) / DET
2346 G2(3,4) = (F11*F22 - F12*F21) / DET
2347 B2 = LOG (XV(I+2)/XV(I))
2348 B3 = XV(I) * (B2 - 1.) + XV(I+2)
2349 GH (1,I) = B2 * G2 (2,2) + B3 * G2 (3,2)
2350 GH (2,I) = B2 * G2 (2,3) + B3 * G2 (3,3)
2351 GH (3,I) = B2 * G2 (2,4) + B3 * G2 (3,4)
2357 ElseIF (I .EQ. NX-1) THEN
2360 GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0
2368 A(K) = XA(I+1, 0) - XA(I, 0)
2370 A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK)
2376 TEM = TEM + A(L) * GB(L,J,I)
2390 FUNCTION CtLhXFRMZ (Z)
2391 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2393 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2394 PARAMETER (MXX = 105)
2395 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2396 COMMON / LhCtINVERT / ZA
2398 DATA TEM, RER / D1, 1E-3 /
2399 DATA ZLOW, ZHIGH, IWRN2 / -10.0, 1.00002, 0 /
2402 IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN
2403 XLA = LOG (XMIN) * 1.5
2405 TEM = CtLhZBRNT (CtLhZFXL, XLA, XLB, EPS, IRT)
2407 CALL CtLhWARNR (IWRN2, 'Z out of range in CtLhXFRMZ, X set=0.',
2408 > 'Z', Z, ZLOW, ZHIGH, 1)
2411 CtLhXFRMZ = EXP(TEM)
2414 FUNCTION CtLhxLi(n,x)
2416 integer NCUT, i,n,m3
2417 real*8 CtLhxLi,Out,x,pi2by6,zeta3,c1,c2
2420 dimension c1(2:m3),c2(2:m3)
2422 data c1/0.75,-0.5833333333333333d0,0.454861111111111d0,
2423 > -0.3680555555555555d0,0.3073611111111111d0,
2424 > -0.2630555555555555d0,0.2294880243764172d0/
2425 data c2/-0.5d0,0.5d0,-0.4583333333333333d0,0.416666666666666d0,
2426 > -0.3805555555555555d0,0.35d0,-0.3241071428571428d0/
2427 data zeta3,pi2by6 /1.20205690315959d0,1.64493406684823d0/
2431 if (abs(x).gt.r) then
2432 PRINT *,'Li: x out of range (-1,1) , x=',x
2436 PRINT *,'Polylogarithm Li undefined for n=',n
2438 elseif (n.eq.0) then
2440 elseif (n.eq.1) then
2442 elseif (n.eq.2) then
2443 !Calculate dilogarithm
2444 !separately for x<0.5 and x>0.5
2445 if (x.ge.(-0.5).and.x.le.0.5) then
2452 elseif (x.eq.0) then
2454 elseif(x.gt.0.5) then !n.eq.2,x>0.5
2456 L = pi2by6 - dlog(x)*dlog(xt)
2463 elseif (x.lt.(-0.5)) then
2466 do while (i.le.NCUT)
2473 elseif (n.eq.3.and.x.ge.0.8) then !use the expansion of Li3 near x=1
2474 L=zeta3+pi2by6*dlog(x)
2478 L=L+(c1(i)+c2(i)*xln1m)*xt**i
2481 else !n>3 or x=3,x<0.8
2485 L=L+r/dble(i)**dble(n)
2491 FUNCTION CtLhZBRLAM (WLLN)
2492 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2493 COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1
2495 TEM1 = 1./ CtLhALPQCD(JRDR, N1, WMULM, I)
2496 TEM2 = 1./ CtLhALPQCD(JRDR, N, VMULM, I)
2497 CtLhZBRLAM = TEM1 - TEM2
2499 FUNCTION CtLhZBRNT(FUNC, X1, X2, TOL, IRT)
2500 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2501 PARAMETER (ITMAX = 1000, EPS = 3.E-12)
2509 IF(FB*FA.GT.0.) THEN
2510 PRINT *, 'Root must be bracketed for CtLhZBRNT. Set = 0'
2517 IF(FB*FC.GT.0.) THEN
2523 IF(ABS(FC).LT.ABS(FB)) THEN
2531 TOL1=2.*EPS*ABS(B)+0.5*TOL
2533 IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN
2537 IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
2545 P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.))
2546 Q=(Q-1.)*(R-1.)*(S-1.)
2550 IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
2563 IF(ABS(D) .GT. TOL1) THEN
2570 PRINT *, 'CtLhZBRNT exceeding maximum iterations.'
2575 FUNCTION CtLhZFRMX (XX)
2576 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2578 PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1)
2579 PARAMETER (MXX = 105)
2580 COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX
2581 DATA IWRN1, HUGE, TINY / 0, 1.E35, 1.E-35 /
2582 F(X) = (XCR-XMIN) * LOG (X/XMIN) + LOG (XCR/XMIN) * (X-XMIN)
2583 D(X) = (XCR-XMIN) / X + LOG (XCR/XMIN)
2585 IF (X .GE. XMIN) THEN
2587 ElseIF (X .GE. D0) THEN
2591 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhZFRMX'
2592 > , 'X', X, TINY, HUGE, 1)
2600 IF (X .GE. XMIN) THEN
2602 ElseIF (X .GE. D0) THEN
2606 CALL CtLhWARNR(IWRN1, 'X out of range in CtLhDZDX '
2607 > , 'X', X, TINY, HUGE, 1)
2614 FUNCTION CtLhZFXL (XL)
2615 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2616 COMMON / LhCtINVERT / ZA
2618 TT = CtLhZFRMX (X) - ZA