]>
Commit | Line | Data |
---|---|---|
4e9e3152 | 1 | SUBROUTINE CtLhALFSET (QS, ALFS) |
2 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
3 | EXTERNAL CtLhRTALF | |
4 | COMMON / CtLhRTALFC / ALFST, JORD, NEFF | |
5 | DATA ALAM, BLAM, ERR / 0.01, 10.0, 0.02 / | |
6 | ALFST = ALFS | |
7 | CALL CtLhParQcd (2, 'ORDR', ORDR, IR1) | |
8 | JORD = Int(ORDR) | |
9 | NEFF = LhCtNFL(QS) | |
10 | EFLLN = CtLhQZBRNT (CtLhRTALF, ALAM, BLAM, ERR, IR2) | |
11 | EFFLAM = QS / EXP (EFLLN) | |
12 | CALL CtLhSETL1 (NEFF, EFFLAM) | |
13 | END | |
14 | FUNCTION CtLhALPI (AMU) | |
15 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
16 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
17 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
18 | LOGICAL SET | |
19 | PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15) | |
20 | DATA IW1, IW2 / 2*0 / | |
21 | IF(.NOT.SET) CALL CtLhLAMCWZ | |
22 | NEFF = LhCtNFL(AMU) | |
23 | ALM = ALAM(NEFF) | |
24 | CtLhALPI = CtLhALPQCD (NORDER, NEFF, AMU/ALM, IRT) | |
25 | IF (IRT .EQ. 1) THEN | |
26 | CALL CtLhWARNR (IW1, 'AMU < ALAM in CtLhALPI', 'AMU', AMU, | |
27 | > ALM, BIG, 1) | |
28 | ELSEIF (IRT .EQ. 2) THEN | |
29 | CALL CtLhWARNR (IW2, 'CtLhALPI > 3; Be aware!', 'CtLhALPI', | |
30 | > CtLhALPI, D0, D1, 0) | |
31 | ENDIF | |
32 | RETURN | |
33 | END | |
34 | FUNCTION CtLhALPQCD (IRDR, NF, RML, IRT) | |
35 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
36 | PARAMETER (D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15) | |
37 | PARAMETER (CG = 3.0d0, TR = 0.5d0, CF = 4.0d0/3.0d0) | |
38 | IRT = 0 | |
39 | IF (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN | |
40 | print *, | |
41 | > 'Order out of range in CtLhALPQCD: IRDR = ', IRDR | |
42 | STOP | |
43 | ENDIF | |
44 | B0 = (11.d0*CG - 2.* NF) / 3.d0 | |
45 | B1 = (34.d0*CG**2 - 10.d0*CG*NF - 6.d0*CF*NF) / 3.d0 | |
46 | RM2 = RML**2 | |
47 | IF (RM2 .LE. 1.) THEN | |
48 | IRT = 1 | |
49 | CtLhALPQCD = 99. | |
50 | RETURN | |
51 | ENDIF | |
52 | ALN = LOG (RM2) | |
53 | AL = 4.d0/ B0 / ALN | |
54 | IF (IRDR .GE. 2) AL = AL * (1.d0-B1*LOG(ALN) / ALN / B0**2) | |
55 | IF (AL .GE. 3.) THEN | |
56 | IRT = 2 | |
57 | ENDIF | |
58 | CtLhALPQCD = AL | |
59 | RETURN | |
60 | END | |
61 | FUNCTION CtLhAMHATF(I) | |
62 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
63 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
64 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
65 | LOGICAL SET | |
66 | IF (.NOT.SET) CALL CtLhLAMCWZ | |
67 | IF ((I.LE.0).OR.(I.GT.9)) THEN | |
68 | print *,'warning I OUT OF RANGE IN CtLhAMHATF' | |
69 | CtLhAMHATF = 0 | |
70 | ELSE | |
71 | CtLhAMHATF = AMHAT(I) | |
72 | ENDIF | |
73 | RETURN | |
74 | END | |
75 | FUNCTION CtLhDXDZ (Z) | |
76 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
77 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
78 | DATA HUGE, IWRN / 1.E20, 0 / | |
79 | ZZ = Z | |
80 | X = CtLhXFRMZ (ZZ) | |
81 | TEM = CtLhDZDX (X) | |
82 | IF (TEM .NE. D0) THEN | |
83 | TMP = D1 / TEM | |
84 | Else | |
85 | CALL CtLhWARNR(IWRN, 'CtLhDXDZ singular in CtLhDXDZ; set=HUGE', | |
86 | > 'Z', Z, D0, D1, 0) | |
87 | TMP = HUGE | |
88 | EndIf | |
89 | CtLhDXDZ = TMP | |
90 | RETURN | |
91 | END | |
92 | SUBROUTINE CtLhEVLPAR (IACT, NAME, VALUE, IRET) | |
93 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
94 | CHARACTER*(*) NAME | |
95 | IRET = 1 | |
96 | IF (IACT .EQ. 0) THEN | |
97 | WRITE ( NINT(VALUE) , 101) | |
98 | 101 FORMAT (/ ' Initiation parameters: Qini, Ipd0, Ihdn ' / | |
99 | > ' Maximum Q, Order of Alpha: Qmax, IKNL ' / | |
100 | > ' X- mesh parameters : Xmin, Xcr, Nx ' / | |
101 | > ' LnQ-mesh parameters : Nt, Jt ' / | |
102 | > ' # of parton flavors : NfMx ' /) | |
103 | IRET = 4 | |
104 | ElseIF (IACT .EQ. 1) THEN | |
105 | CALL CtLhEVLSET (NAME, VALUE, IRET) | |
106 | Else | |
107 | print *,'fatal evlpar' | |
108 | stop | |
109 | EndIf | |
110 | RETURN | |
111 | END | |
112 | SUBROUTINE CtLhEVLSET (NAME, VALUE, IRET) | |
113 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
114 | LOGICAL LSTX | |
115 | CHARACTER*(*) NAME | |
116 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
117 | PARAMETER (MXPN = MXF * 2 + 2) | |
118 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
119 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
120 | COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG | |
121 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
122 | > / PdfSwh / Iset, IpdMod, Iptn0, NuIni | |
123 | IRET = 1 | |
124 | IF (NAME .EQ. 'QINI') THEN | |
125 | IF (VALUE .LE. 0) GOTO 12 | |
126 | QINI = VALUE | |
127 | ElseIF (NAME .EQ. 'IPD0') THEN | |
128 | ITEM = NINT(VALUE) | |
129 | IF (Item .Eq. 10 .or. Item .Eq. 11) GOTO 12 | |
130 | IPD0 = ITEM | |
131 | ElseIF (NAME .EQ. 'IHDN') THEN | |
132 | ITEM = NINT(VALUE) | |
133 | IF (ITEM .LT. -1 .OR. ITEM .GT. 5) GOTO 12 | |
134 | IHDN = ITEM | |
135 | ElseIF (NAME .EQ. 'QMAX') THEN | |
136 | IF (VALUE .LE. QINI) GOTO 12 | |
137 | QMAX = VALUE | |
138 | ElseIF (NAME .EQ. 'IKNL') THEN | |
139 | ITMP = NINT(VALUE) | |
140 | ITEM = ABS(ITMP) | |
141 | IF (ITEM.NE.1.AND.ITEM.NE.2) GOTO 12 | |
142 | IKNL = ITMP | |
143 | ElseIF (NAME .EQ. 'XCR') THEN | |
144 | IF (VALUE .LT. XMIN .OR. VALUE .GT. 10.) GOTO 12 | |
145 | XCR = VALUE | |
146 | LSTX = .FALSE. | |
147 | ElseIF (NAME .EQ. 'XMIN') THEN | |
148 | IF (VALUE .LT. 1D-7 .OR. VALUE .GT. 1D0) GOTO 12 | |
149 | XMIN = VALUE | |
150 | LSTX = .FALSE. | |
151 | ElseIF (NAME .EQ. 'NX') THEN | |
152 | ITEM = NINT(VALUE) | |
153 | IF (ITEM .LT. 10 .OR. ITEM .GT. MXX-1) GOTO 12 | |
154 | NX = ITEM | |
155 | LSTX = .FALSE. | |
156 | ElseIF (NAME .EQ. 'NT') THEN | |
157 | ITEM = NINT(VALUE) | |
158 | IF (ITEM .LT. 2 .OR. ITEM .GT. MXQ) GOTO 12 | |
159 | NT = ITEM | |
160 | ElseIF (NAME .EQ. 'JT') THEN | |
161 | ITEM = NINT(VALUE) | |
162 | IF (ITEM .LT. 1 .OR. ITEM .GT. 5) GOTO 12 | |
163 | JT = ITEM | |
164 | ElseIF (NAME .EQ. 'NFMX') THEN | |
165 | ITEM = NINT(VALUE) | |
166 | IF (ITEM .LT. 1 .OR. ITEM .GT. MXPN) GOTO 12 | |
167 | NfMx = ITEM | |
168 | ElseIF (NAME .EQ. 'IPDMOD') THEN | |
169 | ITEM = NINT(VALUE) | |
170 | IF (Abs(Item) .Gt. 1) GOTO 12 | |
171 | IpdMod = ITEM | |
172 | ElseIF (NAME .EQ. 'IPTN0') THEN | |
173 | ITEM = NINT(VALUE) | |
174 | IF (ABS(ITEM) .GT. MXF) GOTO 12 | |
175 | IPTN0 = ITEM | |
176 | ElseIF (NAME .EQ. 'NUINI') THEN | |
177 | ITEM = NINT(VALUE) | |
178 | IF (ITEM .LE. 0) GOTO 12 | |
179 | NuIni = ITEM | |
180 | Else | |
181 | IRET = 0 | |
182 | EndIf | |
183 | RETURN | |
184 | 12 IRET = 2 | |
185 | RETURN | |
186 | END | |
187 | SUBROUTINE CtLhEVOLVE (FINI, IRET) | |
188 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
189 | include 'parmsetup.inc' | |
190 | LOGICAL LSTX | |
191 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
192 | PARAMETER (MXPN = MXF * 2 + 2) | |
193 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
194 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
195 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
196 | COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG | |
197 | COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) | |
198 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
199 | COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt | |
200 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
201 | COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) | |
202 | DIMENSION QRKP(MXF) | |
203 | DIMENSION JI(-MXF : MXF+1) | |
204 | EXTERNAL LhCtNSRHSP, LhCtNSRHSM, FINI | |
205 | DATA DZER / 0.0 / | |
206 | save nxsave, ntsave, jtsave, ngsave, | |
207 | & xcrsave, xminsave, qinisave, qmaxsave, ientry | |
208 | data ientry / 0 / | |
209 | c | |
210 | call getnset(iset) | |
211 | c | |
212 | ientry = ientry + 1 | |
213 | IRET = 0 | |
214 | IF (IHDN .LE. 4) THEN | |
215 | MXVAL = 2 | |
216 | ElseIF (IHDN .LE. 6) THEN | |
217 | MXVAL = 3 | |
218 | EndIf | |
219 | IF (.NOT. LSTX) CALL CtLhXARRAY | |
220 | CALL CtLhPARPDF (2, 'ALAM', AL, IR) | |
221 | CALL CtLhQARRAY (NINI) | |
222 | NFSN = NFMX + 1 | |
223 | KF = 2 * NFMX + 2 | |
224 | Nelmt = KF * (Nt+1) * (Nx+1) | |
225 | DO 101 IFLV = -NFMX, NFMX+1 | |
226 | JFL = NFMX + IFLV | |
227 | JI(IFLV) = JFL * (NT+1) * (NX+1) | |
228 | 101 CONTINUE | |
229 | DO 31 IZ = 1, NX | |
230 | UPD(JI(0)+IZ+1,iset) = FINI (0, XV(IZ)) | |
231 | UPD(JI(NFSN)+IZ+1,iset) = 0 | |
232 | IF (NFMX .EQ. 0) GOTO 31 | |
233 | DO 331 IFLV = 1, NINI | |
234 | A = FINI ( IFLV, XV(IZ)) | |
235 | B = FINI (-IFLV, XV(IZ)) | |
236 | QRKP (IFLV) = A + B | |
237 | UPD(JI(NFSN)+IZ+1,iset) = | |
238 | > UPD(JI(NFSN)+IZ+1,iset) + QRKP (IFLV) | |
239 | UPD(JI(-IFLV)+IZ+1,iset) = A - B | |
240 | 331 CONTINUE | |
241 | DO 332 IFLV = 1, NINI | |
242 | UPD(JI( IFLV)+IZ+1,iset) = | |
243 | > QRKP(IFLV) - UPD(JI(NFSN)+IZ+1,iset)/NINI | |
244 | 332 CONTINUE | |
245 | 31 CONTINUE | |
246 | DO 21 NEFF = NINI, NFMX | |
247 | IF (IKNL .EQ. 2) CALL CtLhSTUPKL (NEFF) | |
248 | ICNT = NEFF - NINI + 1 | |
249 | IF (NTN(ICNT) .EQ. 0) GOTO 21 | |
250 | NITR = NTN (ICNT) | |
251 | DT = DTN (ICNT) | |
252 | TIN = TLN (ICNT) | |
253 | CALL CtLhSNEVL (IKNL, NX, NITR, JT, DT, TIN, NEFF, | |
254 | > UPD(JI(NFSN)+2,iset), UPD(JI(0)+2,iset), | |
255 | > UPD(JI(NFSN)+1,iset), UPD(JI(0)+1,iset)) | |
256 | IF (NEFF .EQ. 0) GOTO 88 | |
257 | DO 333 IFLV = 1, NEFF | |
258 | CALL CtLhNSEVL (LhCtNSRHSP, IKNL, NX, NITR, JT, DT, TIN, | |
259 | > NEFF, UPD(JI( IFLV)+2,iset), UPD(JI( IFLV)+1,iset)) | |
260 | IF (IFLV .LE. MXVAL) | |
261 | > CALL CtLhNSEVL (LhCtNSRHSM, IKNL, NX, NITR, JT, DT, TIN, | |
262 | > NEFF, UPD(JI(-IFLV)+2,iset), UPD(JI(-IFLV)+1,iset)) | |
263 | DO 55 IS = 0, NITR | |
264 | DO 56 IX = 0, NX | |
265 | TP = UPD (IS*(NX+1) + IX + 1 + JI( IFLV),iset) | |
266 | TS = UPD (IS*(NX+1) + IX + 1 + JI( NFSN),iset) / NEFF | |
267 | TP = TP + TS | |
268 | IF (IKNL .GT. 0) TP = MAX (TP, DZER) | |
269 | IF (IFLV .LE. MXVAL) THEN | |
270 | TM = UPD (IS*(NX+1) + IX + 1 + JI(-IFLV),iset) | |
271 | IF (IKNL .GT. 0) THEN | |
272 | TM = MAX (TM, DZER) | |
273 | TP = MAX (TP, TM) | |
274 | EndIf | |
275 | Else | |
276 | TM = 0. | |
277 | EndIf | |
278 | UPD (JI( IFLV) + IS*(NX+1) + IX + 1,iset) = (TP + TM)/2. | |
279 | UPD (JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = (TP - TM)/2. | |
280 | 56 CONTINUE | |
281 | 55 CONTINUE | |
282 | 333 CONTINUE | |
283 | DO 334 IFLV = NEFF + 1, NFMX | |
284 | DO 57 IS = 0, NITR | |
285 | DO 58 IX = 0, NX | |
286 | UPD(JI( IFLV) + IS*(NX+1) + IX + 1,iset) = 0 | |
287 | UPD(JI(-IFLV) + IS*(NX+1) + IX + 1,iset) = 0 | |
288 | 58 CONTINUE | |
289 | 57 CONTINUE | |
290 | 334 CONTINUE | |
291 | 88 CONTINUE | |
292 | IF (NFMX .EQ. NEFF) GOTO 21 | |
293 | DO 335 IFLV = -NFMX, NFMX+1 | |
294 | JI(IFLV) = JI(IFLV) + NITR * (NX+1) | |
295 | 335 CONTINUE | |
296 | CALL CtLhHQRK (NX, TT, NEFF+1, UPD(JI(0)+2,iset), | |
297 | > UPD(JI(NEFF+1)+2,iset)) | |
298 | DO 32 IZ = 1, NX | |
299 | QRKP (NEFF+1) = 2. * UPD(JI( NEFF+1) + IZ + 1,iset) | |
300 | UPD (JI(NFSN)+IZ+1,iset) = UPD (JI(NFSN)+IZ+1,iset) | |
301 | > + QRKP (NEFF+1) | |
302 | VS00 = UPD (JI(NFSN)+IZ+1,iset) / (NEFF+1) | |
303 | UPD(JI( NEFF+1) + IZ + 1,iset) = QRKP(NEFF+1) - VS00 | |
304 | DO 321 IFL = 1, NEFF | |
305 | A = UPD(JI( IFL)+IZ+1,iset) | |
306 | B = UPD(JI(-IFL)+IZ+1,iset) | |
307 | QRKP(IFL) = A + B | |
308 | UPD(JI( IFL)+IZ+1,iset) = QRKP(IFL) - VS00 | |
309 | IF (IFL .LE. MXVAL) UPD(JI(-IFL)+IZ+1,iset) = A - B | |
310 | 321 CONTINUE | |
311 | 32 CONTINUE | |
312 | 21 CONTINUE | |
313 | if(ientry .eq. 1) then | |
314 | nxsave = nx | |
315 | ntsave = nt | |
316 | jtsave = jt | |
317 | ngsave = ng | |
318 | xcrsave = xcr | |
319 | xminsave = xmin | |
320 | qinisave = qini | |
321 | qmaxsave = qmax | |
322 | endif | |
323 | if((nx .ne. nxsave) .or. | |
324 | & (nt .ne. ntsave) .or. | |
325 | & (jt .ne. jtsave) .or. | |
326 | & (ng .ne. ngsave) .or. | |
327 | & (xcr .ne. xcrsave) .or. | |
328 | & (xmin .ne. xminsave) .or. | |
329 | & (qini .ne. qinisave) .or. | |
330 | & (qmax .ne. qmaxsave)) then | |
331 | write(6,669) nx, nt, jt, ng, xcr, xmin, | |
332 | & qini, qmax, ientry | |
333 | 669 format(1x,'evolve.f: nx,nt,jt,ng=',4i4, | |
334 | & ' xcr,xmin=',2f9.6, | |
335 | & ' qini, qmax',f7.4,1x,e12.5,' ientry=',i6) | |
336 | nxsave = nx | |
337 | ntsave = nt | |
338 | jtsave = jt | |
339 | ngsave = ng | |
340 | qinisave = qini | |
341 | qmaxsave = qmax | |
342 | xcrsave = xcr | |
343 | xminsave = xmin | |
344 | endif | |
345 | Return | |
346 | End | |
347 | * | |
348 | FUNCTION CtLhFINTRP (FF, X0, DX, NX, XV, ERR, IR) | |
349 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
350 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
351 | PARAMETER (MX = 3) | |
352 | DIMENSION FF (0:NX), XX(MX) | |
353 | DATA SML, XX / 1.D-5, 0., 1.0, 2.0 / | |
354 | DATA IW1, IW3, IW5 / 3 * 0 / | |
355 | IR = 0 | |
356 | X = XV | |
357 | ERR = 0. | |
358 | ANX = NX | |
359 | CtLhFINTRP = 0. | |
360 | IF (NX .LT. 1) THEN | |
361 | CALL CtLhWARNI(IW1, 'Nx < 1, error in CtLhFINTRP.', | |
362 | > 'NX', NX, 1, 256, 1) | |
363 | IR = 1 | |
364 | RETURN | |
365 | ELSE | |
366 | MNX = MIN(NX+1, MX) | |
367 | ENDIF | |
368 | IF (DX .LE. 0) THEN | |
369 | CALL CtLhWARNR(IW3, 'DX < 0, error in CtLhFINTRP.', | |
370 | > 'DX', DX, D0, D1, 1) | |
371 | IR = 2 | |
372 | RETURN | |
373 | ENDIF | |
374 | XM = X0 + DX * NX | |
375 | IF (X .LT. X0-SML .OR. X .GT. XM+SML) THEN | |
376 | CALL CtLhWARNR(IW5, | |
377 | > 'X out of range in CtLhFINTRP, Extrapolation used.', | |
378 | > 'X',X,X0,XM,1) | |
379 | IR = 3 | |
380 | ENDIF | |
381 | TX = (X - X0) / DX | |
382 | IF (TX .LE. 1.) THEN | |
383 | IX = 0 | |
384 | ELSEIF (TX .GE. ANX-1.) THEN | |
385 | IX = NX - 2 | |
386 | ELSE | |
387 | IX = Int(TX) | |
388 | ENDIF | |
389 | DDX = TX - IX | |
390 | CALL CtLhRATINT (XX, FF(IX), MNX, DDX, TEM, ERR) | |
391 | CtLhFINTRP = TEM | |
392 | RETURN | |
393 | END | |
394 | FUNCTION CtLhGausInt(F,XL,XR,AERR,RERR,ERR,IRT) | |
395 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
396 | DIMENSION XLIMS(100), R(93), W(93) | |
397 | INTEGER PTR(4),NORD(4) | |
398 | external f | |
399 | DATA PTR,NORD/4,10,22,46, 6,12,24,48/ | |
400 | DATA R/.2386191860,.6612093865,.9324695142, | |
401 | 1 .1252334085,.3678314990,.5873179543,.7699026742,.9041172563, | |
402 | 1 .9815606342,.0640568929,.1911188675,.3150426797,.4337935076, | |
403 | 1 .5454214714,.6480936519,.7401241916,.8200019860,.8864155270, | |
404 | 1 .9382745520,.9747285560,.9951872200,.0323801710,.0970046992, | |
405 | 1 .1612223561,.2247637903,.2873624873,.3487558863,.4086864820, | |
406 | 1 .4669029048,.5231609747,.5772247261,.6288673968,.6778723796, | |
407 | 1 .7240341309,.7671590325,.8070662040,.8435882616,.8765720203, | |
408 | 1 .9058791367,.9313866907,.9529877032,.9705915925,.9841245837, | |
409 | 1 .9935301723,.9987710073,.0162767488,.0488129851,.0812974955, | |
410 | 1 .1136958501,.1459737146,.1780968824,.2100313105,.2417431561, | |
411 | 1 .2731988126,.3043649444,.3352085229,.3656968614,.3957976498, | |
412 | 1 .4254789884,.4547094222,.4834579739,.5116941772,.5393881083, | |
413 | 1 .5665104186,.5930323648,.6189258401,.6441634037,.6687183100, | |
414 | 1 .6925645366,.7156768123,.7380306437,.7596023411,.7803690438, | |
415 | 1 .8003087441,.8194003107,.8376235112,.8549590334,.8713885059, | |
416 | 1 .8868945174,.9014606353,.9150714231,.9277124567,.9393703398, | |
417 | 1 .9500327178,.9596882914,.9683268285,.9759391746,.9825172636, | |
418 | 1 .9880541263,.9925439003,.9959818430,.9983643759,.9996895039/ | |
419 | DATA W/.4679139346,.3607615730,.1713244924, | |
420 | 1 .2491470458,.2334925365,.2031674267,.1600783285,.1069393260, | |
421 | 1 .0471753364,.1279381953,.1258374563,.1216704729,.1155056681, | |
422 | 1 .1074442701,.0976186521,.0861901615,.0733464814,.0592985849, | |
423 | 1 .0442774388,.0285313886,.0123412298,.0647376968,.0644661644, | |
424 | 1 .0639242386,.0631141923,.0620394232,.0607044392,.0591148397, | |
425 | 1 .0572772921,.0551995037,.0528901894,.0503590356,.0476166585, | |
426 | 1 .0446745609,.0415450829,.0382413511,.0347772226,.0311672278, | |
427 | 1 .0274265097,.0235707608,.0196161605,.0155793157,.0114772346, | |
428 | 1 .0073275539,.0031533461,.0325506145,.0325161187,.0324471637, | |
429 | 1 .0323438226,.0322062048,.0320344562,.0318287589,.0315893308, | |
430 | 1 .0313164256,.0310103326,.0306713761,.0302999154,.0298963441, | |
431 | 1 .0294610900,.0289946142,.0284974111,.0279700076,.0274129627, | |
432 | 1 .0268268667,.0262123407,.0255700360,.0249006332,.0242048418, | |
433 | 1 .0234833991,.0227370697,.0219666444,.0211729399,.0203567972, | |
434 | 1 .0195190811,.0186606796,.0177825023,.0168854799,.0159705629, | |
435 | 1 .0150387210,.0140909418,.0131282296,.0121516047,.0111621020, | |
436 | 1 .0101607705,.0091486712,.0081268769,.0070964708,.0060585455, | |
437 | 1 .0050142027,.0039645543,.0029107318,.0018539608,.0007967921/ | |
438 | DATA TOLABS,TOLREL,NMAX/1.E-35,5.E-4,100/ | |
439 | TOLABS=AERR | |
440 | TOLREL=RERR | |
441 | ||
442 | CtLhGausInt=0. | |
443 | NLIMS=2 | |
444 | XLIMS(1)=XL | |
445 | XLIMS(2)=XR | |
446 | 10 AA=(XLIMS(NLIMS)-XLIMS(NLIMS-1))/2D0 | |
447 | BB=(XLIMS(NLIMS)+XLIMS(NLIMS-1))/2D0 | |
448 | TVAL=0. | |
449 | DO 15 I=1,3 | |
450 | 15 TVAL=TVAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I))) | |
451 | TVAL=TVAL*AA | |
452 | DO 25 J=1,4 | |
453 | VAL=0. | |
454 | DO 20 I=PTR(J),PTR(J)-1+NORD(J) | |
455 | 20 VAL=VAL+W(I)*(F(BB+AA*R(I))+F(BB-AA*R(I))) | |
456 | VAL=VAL*AA | |
457 | TOL=MAX(TOLABS,TOLREL*ABS(VAL)) | |
458 | IF (ABS(TVAL-VAL).LT.TOL) THEN | |
459 | CtLhGausInt=CtLhGausInt+VAL | |
460 | NLIMS=NLIMS-2 | |
461 | IF (NLIMS.NE.0) GO TO 10 | |
462 | RETURN | |
463 | END IF | |
464 | 25 TVAL=VAL | |
465 | IF (NMAX.EQ.2) THEN | |
466 | CtLhGausInt=VAL | |
467 | RETURN | |
468 | END IF | |
469 | IF (NLIMS.GT.(NMAX-2)) THEN | |
470 | write(*,50) CtLhGausInt,NMAX,BB-AA,BB+AA | |
471 | RETURN | |
472 | END IF | |
473 | XLIMS(NLIMS+1)=BB | |
474 | XLIMS(NLIMS+2)=BB+AA | |
475 | XLIMS(NLIMS)=BB | |
476 | NLIMS=NLIMS+2 | |
477 | GO TO 10 | |
478 | 50 FORMAT (' CtLhGausInt FAILS, CtLhGausInt,NMAX,XL,XR=', | |
479 | > G15.7,I5,2G15.7) | |
480 | END | |
481 | SUBROUTINE CtLhHINTEG (NX, F, H) | |
482 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
483 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
484 | PARAMETER (MXPN = MXF * 2 + 2) | |
485 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
486 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
487 | COMMON / LhCtHINTEC / GH(NDG, MXX) | |
488 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
489 | DIMENSION F(NX), H(NX), G(MXX) | |
490 | DZ = 1D0 / (NX-1) | |
491 | DO 20 I = 1, NX-2 | |
492 | NP = NX - I + 1 | |
493 | TEM = GH(1,I)*F(I) + GH(2,I)*F(I+1) + GH(3,I)*F(I+2) | |
494 | DO 30 KZ = 3, NP | |
495 | IY = I + KZ - 1 | |
496 | W = XA(I,1) / XA(IY,1) | |
497 | G(KZ) = DXTZ(IY)*(F(IY)-W*F(I))/(1.-W) | |
498 | 30 CONTINUE | |
499 | HTEM = CtLhSMPSNA (NP-2, DZ, G(3), ERR) | |
500 | TEM1 = F(I) * ELY(I) | |
501 | H(I) = TEM + HTEM + TEM1 | |
502 | 20 CONTINUE | |
503 | H(NX-1) = F(NX) - F(NX-1) + F(NX-1) * (ELY(NX-1) - XA(NX-1,0)) | |
504 | H(NX) = 0 | |
505 | RETURN | |
506 | END | |
507 | SUBROUTINE CtLhHQRK (NX, TT, NQRK, Y, F) | |
508 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
509 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
510 | DIMENSION Y(NX), F(NX) | |
511 | IF (NX .GT. 1) GOTO 11 | |
512 | 11 CONTINUE | |
513 | DO 230 IZ = 1, NX | |
514 | IF (NX .GT. 1) THEN | |
515 | F(IZ) = 0 | |
516 | GOTO 230 | |
517 | EndIf | |
518 | 230 CONTINUE | |
519 | RETURN | |
520 | END | |
521 | SUBROUTINE CtLhINTEGR (NX, M, F, G, IR) | |
522 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
523 | CHARACTER MSG*80 | |
524 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
525 | PARAMETER (MXPN = MXF * 2 + 2) | |
526 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
527 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
528 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
529 | COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) | |
530 | DIMENSION F(NX), G(NX) | |
531 | DATA IWRN1, IWRN2 / 0, 0 / | |
532 | IRR = 0 | |
533 | IF (NX .LT. 1 .OR. XA(NX-1,1) .EQ. 0D0) THEN | |
534 | MSG = 'NX out of range in CtLhINTEGR call' | |
535 | CALL CtLhWARNI (IWRN1, MSG, 'NX', NX, 0, MXX, 0) | |
536 | IRR = 1 | |
537 | EndIf | |
538 | IF (M .LT. M1 .OR. M .GT. M2) THEN | |
539 | MSG ='Exponent M out of range in CtLhINTEGR' | |
540 | CALL CtLhWARNI (IWRN2, MSG, 'M', M, M1, M2, 1) | |
541 | IRR = 2 | |
542 | EndIf | |
543 | G(NX) = 0D0 | |
544 | TEM = H(1, NX-1, -M) * F(NX-2) + H(2, NX-1, -M) * F(NX-1) | |
545 | > + H(3, NX-1, -M) * F(NX) | |
546 | IF (M .EQ. 0) THEN | |
547 | G(NX-1) = TEM | |
548 | Else | |
549 | G(NX-1) = TEM * XA(NX-1, M) | |
550 | EndIf | |
551 | DO 10 I = NX-2, 2, -1 | |
552 | TEM = TEM + H(1,I,-M)*F(I-1) + H(2,I,-M)*F(I) | |
553 | > + H(3,I,-M)*F(I+1) + H(4,I,-M)*F(I+2) | |
554 | IF (M .EQ. 0) THEN | |
555 | G(I) = TEM | |
556 | Else | |
557 | G(I) = TEM * XA(I, M) | |
558 | EndIf | |
559 | 10 CONTINUE | |
560 | TEM = TEM + H(2,1,-M)*F(1) + H(3,1,-M)*F(2) + H(4,1,-M)*F(3) | |
561 | IF (M .EQ. 0) THEN | |
562 | G(1) = TEM | |
563 | Else | |
564 | G(1) = TEM * XA(1, M) | |
565 | EndIf | |
566 | IR = IRR | |
567 | RETURN | |
568 | END | |
569 | SUBROUTINE CtLhKERNEL | |
570 | >(XX, FF1, FG1, GF1, GG1, PNSP, PNSM, FF2, FG2, GF2, GG2, NFL, IRT) | |
571 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
572 | PARAMETER (PI = 3.141592653589793d0, PI2 = PI ** 2) | |
573 | PARAMETER (D0 = 0.0, D1 = 1.0) | |
574 | DATA CF, CG, TR, IWRN / 1.33333333333333d0, 3.0d0, 0.5d0, 0 / | |
575 | IRT = 0 | |
576 | TRNF = TR * NFL | |
577 | X = XX | |
578 | IF (X .LE. 0. .OR. X .GE. 1.) THEN | |
579 | CALL CtLhWARNR(IWRN, 'X out of range in CtLhKERNEL', 'X', X, | |
580 | > D0, D1, 1) | |
581 | IRT = 1 | |
582 | RETURN | |
583 | EndIf | |
584 | XI = 1./ X | |
585 | X2 = X ** 2 | |
586 | XM1I = 1./ (1.- X) | |
587 | XP1I = 1./ (1.+ X) | |
588 | XLN = LOG (X) | |
589 | XLN2 = XLN ** 2 | |
590 | XLN1M = LOG (1.- X) | |
591 | SPEN2 = CtLhSPENC2 (X) | |
592 | FFP = (1.+ X2) * XM1I | |
593 | FGP = (2.- 2.* X + X2) / X | |
594 | GFP = 1. - 2.* X + 2.* X2 | |
595 | GGP = XM1I + XI - 2. + X - X2 | |
596 | FFM = (1.+ X2) * XP1I | |
597 | FGM = - (2.+ 2.* X + X2) / X | |
598 | GFM = 1. + 2.* X + 2.* X2 | |
599 | GGM = XP1I - XI - 2. - X - X2 | |
600 | FF1 = CF * FFP * (1.- X) | |
601 | FG1 = CF * FGP * X | |
602 | GF1 = 2.* TRNF * GFP | |
603 | GG1 = 2.* CG * GGP * X * (1.-X) | |
604 | PCF2 = -2.* FFP *XLN*XLN1M - (3.*XM1I + 2.*X)*XLN | |
605 | > - (1.+X)/2.*XLN2 - 5.*(1.-X) | |
606 | PCFG = FFP * (XLN2 + 11.*XLN/3.+ 67./9.- PI**2 / 3.) | |
607 | > + 2.*(1.+X) * XLN + 40.* (1.-X) / 3. | |
608 | PCFT = (FFP * (- XLN - 5./3.) - 2.*(1.-X)) * 2./ 3. | |
609 | PQQB = 2.* FFM * SPEN2 + 2.*(1.+X)*XLN + 4.*(1.-X) | |
610 | PQQB = (CF**2-CF*CG/2.) * PQQB | |
611 | PQQ2 = CF**2 * PCF2 + CF*CG * PCFG / 2. + CF*TRNF * PCFT | |
612 | PNSP = (PQQ2 + PQQB) * (1.-X) | |
613 | PNSM = (PQQ2 - PQQB) * (1.-X) | |
614 | FFCF2 = - 1. + X + (1.- 3.*X) * XLN / 2. - (1.+ X) * XLN2 / 2. | |
615 | > - FFP * (3.* XLN / 2. + 2.* XLN * XLN1M) | |
616 | > + FFM * 2.* SPEN2 | |
617 | FFCFG = 14./3.* (1.-X) | |
618 | > + FFP * (11./6.* XLN + XLN2 / 2. + 67./18. - PI2 / 6.) | |
619 | > - FFM * SPEN2 | |
620 | FFCFT = - 16./3. + 40./3.* X + (10.* X + 16./3.* X2 + 2.) * XLN | |
621 | > - 112./9.* X2 + 40./9./X - 2.* (1.+ X) * XLN2 | |
622 | > - FFP * (10./9. + 2./3. * XLN) | |
623 | FGCF2 = - 5./2.- 7./2.* X + (2.+ 7./2.* X) * XLN + (X/2.-1.)*XLN2 | |
624 | > - 2.* X * XLN1M | |
625 | > - FGP * (3.* XLN1M + XLN1M ** 2) | |
626 | FGCFG = 28./9. + 65./18.* X + 44./9. * X2 - (12.+ 5.*X + 8./3.*X2) | |
627 | > * XLN + (4.+ X) * XLN2 + 2.* X * XLN1M | |
628 | > + FGP * (-2.*XLN*XLN1M + XLN2/2. + 11./3.*XLN1M + XLN1M**2 | |
629 | > - PI2/6. + 0.5) | |
630 | > + FGM * SPEN2 | |
631 | FGCFT = -4./3.* X - FGP * (20./9.+ 4./3.*XLN1M) | |
632 | GFCFT = 4.- 9.*X + (-1.+ 4.*X)*XLN + (-1.+ 2.*X)*XLN2 + 4.*XLN1M | |
633 | > + GFP * (-4.*XLN*XLN1M + 4.*XLN + 2.*XLN2 - 4.*XLN1M | |
634 | > + 2.*XLN1M**2 - 2./3.* PI2 + 10.) | |
635 | GFCGT = 182./9.+ 14./9.*X + 40./9./X + (136./3.*X - 38./3.)*XLN | |
636 | > - 4.*XLN1M - (2.+ 8.*X)*XLN2 | |
637 | > + GFP * (-XLN2 + 44./3.*XLN - 2.*XLN1M**2 + 4.*XLN1M | |
638 | > + PI2/3. - 218./9.) | |
639 | > + GFM * 2. * SPEN2 | |
640 | GGCFT = -16.+ 8.*X + 20./3.*X2 + 4./3./X + (-6.-10.*X)*XLN | |
641 | > - 2.* (1.+ X) * XLN2 | |
642 | GGCGT = 2.- 2.*X + 26./9.*X2 - 26./9./X - 4./3.*(1.+X)*XLN | |
643 | > - GGP * 20./9. | |
644 | GGCG2 = 27./2.*(1.-X) + 67./9.*(X2-XI) + 4.*(1.+X)*XLN2 | |
645 | > + (-25.+ 11.*X - 44.*X2)/3.*XLN | |
646 | > + GGP * (67./9.- 4.*XLN*XLN1M + XLN2 - PI2/3.) | |
647 | > + GGM * 2.* SPEN2 | |
648 | FF2 = CF * TRNF * FFCFT + CF ** 2 * FFCF2 + CF * CG * FFCFG | |
649 | FG2 = CF * TRNF * FGCFT + CF ** 2 * FGCF2 + CF * CG * FGCFG | |
650 | GF2 = CF * TRNF * GFCFT + CG * TRNF * GFCGT | |
651 | GG2 = CF * TRNF * GGCFT + CG ** 2 * GGCG2 + CG * TRNF * GGCGT | |
652 | XLG = (LOG(1./(1.-X)) + 1.) | |
653 | XG2 = XLG ** 2 | |
654 | FF2 = FF2 * X * (1.- X) | |
655 | FG2 = FG2 * X / XG2 | |
656 | GF2 = GF2 * X / XG2 | |
657 | GG2 = GG2 * X * (1.- X) | |
658 | RETURN | |
659 | END | |
660 | SUBROUTINE CtLhLAMCWZ | |
661 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
662 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
663 | LOGICAL SET | |
664 | CALL CtLhSETL1 (NF, AL) | |
665 | END | |
666 | FUNCTION LhCtNAMQCD(NNAME) | |
667 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
668 | CHARACTER NNAME*(*), NAME*8 | |
669 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
670 | LOGICAL SET | |
671 | CHARACTER ONECH*(1) | |
672 | ONECH = '0' | |
673 | IASC0 = ICHAR(ONECH) | |
674 | NAME = NNAME | |
675 | LhCtNAMQCD=0 | |
676 | IF ( (NAME .EQ. 'ALAM') .OR. (NAME .EQ. 'LAMB') .OR. | |
677 | 1 (NAME .EQ. 'LAM') .OR. (NAME .EQ. 'LAMBDA') ) | |
678 | 2 LhCtNAMQCD=1 | |
679 | IF ( (NAME .EQ. 'NFL') .OR. (NAME(1:3) .EQ. '#FL') .OR. | |
680 | 1 (NAME .EQ. '# FL') ) | |
681 | 2 LhCtNAMQCD=2 | |
682 | DO 10 I=1, 9 | |
683 | IF (NAME .EQ. 'M'//CHAR(I+IASC0)) | |
684 | 1 LhCtNAMQCD=I+2 | |
685 | 10 CONTINUE | |
686 | DO 20 I= 0, NF | |
687 | IF (NAME .EQ. 'LAM'//CHAR(I+IASC0)) | |
688 | 1 LhCtNAMQCD=I+13 | |
689 | 20 CONTINUE | |
690 | IF (NAME(:3).EQ.'ORD' .OR. NAME(:3).EQ.'NRD') LhCtNAMQCD = 24 | |
691 | RETURN | |
692 | END | |
693 | FUNCTION LhCtNFL(AMU) | |
694 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
695 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
696 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
697 | LOGICAL SET | |
698 | IF (.NOT. SET) CALL CtLhLAMCWZ | |
699 | LhCtNFL = NF - NHQ | |
700 | IF ((LhCtNFL .EQ. NF) .OR. (AMU .LE. AMN)) GOTO 20 | |
701 | DO 10 I = NF - NHQ + 1, NF | |
702 | IF (AMU .GE. AMHAT(I)) THEN | |
703 | LhCtNFL = I | |
704 | ELSE | |
705 | GOTO 20 | |
706 | ENDIF | |
707 | 10 CONTINUE | |
708 | 20 RETURN | |
709 | END | |
710 | SUBROUTINE CtLhNSEVL (RHS, IKNL,NX,NT,JT,DT,TIN,NEFF,U0,UN) | |
711 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
712 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
713 | PARAMETER (MXPN = MXF * 2 + 2) | |
714 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
715 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
716 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
717 | DIMENSION U0(NX), UN(0:NX, 0:NT) | |
718 | DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX) | |
719 | external rhs | |
720 | DDT = DT / JT | |
721 | IF (NX .GT. MXX) THEN | |
722 | WRITE (*,*) 'Nx =', NX, ' greater than Max pts in CtLhNSEVL.' | |
723 | STOP 'Program stopped in CtLhNSEVL' | |
724 | EndIf | |
725 | c ++ remove unused quantities (jcp) | |
726 | c ++ TMD = TIN + DT * NT / 2. | |
727 | c ++ AMU = EXP(TMD) | |
728 | c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU) | |
729 | c ++ TLAM = TMD - TEM | |
730 | DO 9 IX = 1, NX | |
731 | UN(IX, 0) = U0(IX) | |
732 | 9 CONTINUE | |
733 | UN(0, 0) = 3D0*U0(1) - 3D0*U0(2) - U0(1) | |
734 | TT = TIN | |
735 | DO 10 IZ = 1, NX | |
736 | Y0(IZ) = U0(IZ) | |
737 | 10 CONTINUE | |
738 | DO 20 IS = 1, NT | |
739 | DO 202 JS = 1, JT | |
740 | IRND = (IS-1) * JT + JS | |
741 | IF (IRND .EQ. 1) THEN | |
742 | CALL RHS (TT, Neff, Y0, F0) | |
743 | DO 250 IZ = 1, NX | |
744 | Y0(IZ) = Y0(IZ) + DDT * F0(IZ) | |
745 | 250 CONTINUE | |
746 | TT = TT + DDT | |
747 | CALL RHS (TT, NEFF, Y0, F1) | |
748 | DO 251 IZ = 1, NX | |
749 | Y1(IZ) = U0(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0 | |
750 | 251 CONTINUE | |
751 | Else | |
752 | CALL RHS (TT, NEFF, Y1, F1) | |
753 | DO 252 IZ = 1, NX | |
754 | YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0 | |
755 | 252 CONTINUE | |
756 | TT = TT + DDT | |
757 | CALL RHS (TT, NEFF, YP, FP) | |
758 | DO 253 IZ = 1, NX | |
759 | Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0 | |
760 | F0(IZ) = F1(IZ) | |
761 | 253 CONTINUE | |
762 | EndIf | |
763 | 202 CONTINUE | |
764 | DO 260 IZ = 1, NX | |
765 | UN (IZ, IS) = Y1(IZ) | |
766 | 260 CONTINUE | |
767 | UN(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) | |
768 | 20 CONTINUE | |
769 | RETURN | |
770 | END | |
771 | SUBROUTINE LhCtNSRHSM (TT, NEFF, FI, FO) | |
772 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
773 | LOGICAL LSTX | |
774 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
775 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
776 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
777 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
778 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
779 | COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX), | |
780 | > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB | |
781 | COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) | |
782 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
783 | DIMENSION G1(MXX), FI(NX), FO(NX) | |
784 | DIMENSION W0(MXX), W1(MXX), WH(MXX) | |
785 | S = EXP(TT) | |
786 | Q = AL * EXP (S) | |
787 | CPL = CtLhALPI(Q) | |
788 | CPL2= CPL ** 2 / 2. * S | |
789 | CPL = CPL * S | |
790 | CALL CtLhINTEGR (NX, 0, FI, W0, IR1) | |
791 | CALL CtLhINTEGR (NX, 1, FI, W1, IR2) | |
792 | CALL CtLhHINTEG (NX, FI, WH) | |
793 | DO 230 IZ = 1, NX | |
794 | FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ)) | |
795 | FO(IZ) = CPL * FO(IZ) | |
796 | 230 CONTINUE | |
797 | IF (IKNL .EQ. 2) THEN | |
798 | DZ = 1./ (NX - 1) | |
799 | DO 21 IX = 1, NX-1 | |
800 | NP = NX - IX + 1 | |
801 | IS = NP | |
802 | DO 31 KZ = 2, NP | |
803 | IY = IX + KZ - 1 | |
804 | IT = NX - IY + 1 | |
805 | XY = ZZ (IS, IT) | |
806 | G1(KZ) = PNS (IS,IT) * (FI(IY) - XY * FI(IX)) | |
807 | 31 CONTINUE | |
808 | TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR) | |
809 | TMP2 = (TEM1 - FI(IX) * ANSM(IX)) * CPL2 | |
810 | FO(IX) = FO(IX) + TMP2 | |
811 | 21 CONTINUE | |
812 | EndIf | |
813 | RETURN | |
814 | END | |
815 | SUBROUTINE LhCtNSRHSP (TT, NEFF, FI, FO) | |
816 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
817 | LOGICAL LSTX | |
818 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
819 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
820 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
821 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
822 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
823 | COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX), | |
824 | > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB | |
825 | COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) | |
826 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
827 | DIMENSION G1(MXX), FI(NX), FO(NX) | |
828 | DIMENSION W0(MXX), W1(MXX), WH(MXX) | |
829 | S = EXP(TT) | |
830 | Q = AL * EXP (S) | |
831 | CPL = CtLhALPI(Q) | |
832 | CPL2= CPL ** 2 / 2. * S | |
833 | CPL = CPL * S | |
834 | CALL CtLhINTEGR (NX, 0, FI, W0, IR1) | |
835 | CALL CtLhINTEGR (NX, 1, FI, W1, IR2) | |
836 | CALL CtLhHINTEG (NX, FI, WH) | |
837 | DO 230 IZ = 1, NX | |
838 | FO(IZ) = 2.* FI(IZ) + 4./3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ)) | |
839 | FO(IZ) = CPL * FO(IZ) | |
840 | 230 CONTINUE | |
841 | IF (IKNL .EQ. 2) THEN | |
842 | DZ = 1./ (NX - 1) | |
843 | DO 21 IX = 1, NX-1 | |
844 | NP = NX - IX + 1 | |
845 | DO 31 KZ = 2, NP | |
846 | IY = IX + KZ - 1 | |
847 | XY = ZZ (NX-IX+1, NX-IY+1) | |
848 | G1(KZ) = PNS (IX,IY) * (FI(IY) - XY * FI(IX)) | |
849 | 31 CONTINUE | |
850 | TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR) | |
851 | TMP2 = (TEM1 + FI(IX) * (-ANSP(IX) + ZQQB)) * CPL2 | |
852 | FO(IX) = FO(IX) + TMP2 | |
853 | 21 CONTINUE | |
854 | EndIf | |
855 | RETURN | |
856 | END | |
857 | FUNCTION CtLhPARDIS (IPRTN, XX, QQ) | |
858 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
859 | include 'parmsetup.inc' | |
860 | Character Msg*80 | |
861 | LOGICAL LSTX | |
862 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
863 | PARAMETER (MXPN = MXF * 2 + 2) | |
864 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
865 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
866 | PARAMETER (Smll = 1D-9) | |
867 | parameter(nqvec = 4) | |
868 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
869 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
870 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
871 | COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG | |
872 | COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) | |
873 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
874 | COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt | |
875 | COMMON / LhCtCOMQMS / VALQMS(9) | |
876 | dimension fvec(4), fij(4) | |
877 | dimension xvpow(0:mxx) | |
878 | Data Iwrn1, Iwrn2, Iwrn3, OneP / 3*0, 1.00001 / | |
879 | data xpow / 0.3d0 / !**** choice of interpolation variable | |
880 | data nxsave / 0 / | |
881 | save xvpow, nxsave | |
882 | save xlast, qlast | |
883 | save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3 | |
884 | save const1 , const2, const3, const4, const5, const6 | |
885 | save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet | |
886 | c | |
887 | call getnset(iset) | |
888 | c | |
889 | if(nx .ne. nxsave) then | |
890 | xvpow(0) = 0D0 | |
891 | do i = 1, nx | |
892 | xvpow(i) = xv(i)**xpow | |
893 | enddo | |
894 | nxsave = nx | |
895 | endif | |
896 | ||
897 | X = XX | |
898 | Q = QQ | |
899 | ||
900 | if((x.lt.xmin).or.(x.gt.1.d0)) print 98,x | |
901 | 98 format(' WARNING: X=',e12.5,' OUT OF RANGE') | |
902 | if((q.lt.qini).or.(q.gt.qmax)) print 99,q | |
903 | 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE') | |
904 | ||
905 | c enforce threshold early to improve speed... | |
906 | ii = iabs(IPRTN) | |
907 | if(ii .ne. 0) then | |
908 | if(QQ .lt. VALQMS(ii) ) then | |
909 | ctlhpardis = 0.d0 | |
910 | return | |
911 | endif | |
912 | endif | |
913 | ||
914 | c force pardis = 0.0d0 at exactly =1.0d0 - added mrw 10/May/06 | |
915 | if(xx .eq. 1.0d0) then | |
916 | ctlhpardis = 0.0d0 | |
917 | return | |
918 | endif | |
919 | ||
920 | c skip the initialization in x if same as in the previous call. | |
921 | if(x .eq. xlast) goto 100 | |
922 | xlast = x | |
923 | ||
924 | JLx = -1 | |
925 | JU = Nx+1 | |
926 | 11 If (JU-JLx .GT. 1) Then | |
927 | JM = (JU+JLx) / 2 | |
928 | If (X .Ge. XV(JM)) Then | |
929 | JLx = JM | |
930 | Else | |
931 | JU = JM | |
932 | Endif | |
933 | Goto 11 | |
934 | Endif | |
935 | If (JLx .LE. -1) Then | |
936 | Print '(A,1pE12.4)','Severe error: x <= 0 in ParDis x=', x | |
937 | Stop | |
938 | ElseIf (JLx .Eq. 0) Then | |
939 | Jx = 0 | |
940 | Msg = '0 < X < Xmin in ParDis; extrapolation used!' | |
941 | CALL CtLhWARNR (IWRN1, Msg, 'X', X, Xmin, 1D0, 1) | |
942 | Elseif (JLx .LE. Nx-2) Then | |
943 | Jx = JLx - 1 | |
944 | Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then | |
945 | Jx = JLx - 2 | |
946 | Else | |
947 | Print '(A,1pE12.4)','Severe error: x > 1 in ParDis x=', x | |
948 | Stop | |
949 | Endif | |
950 | ss = x**xpow | |
951 | If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then | |
952 | svec1 = xvpow(jx) | |
953 | svec2 = xvpow(jx+1) | |
954 | svec3 = xvpow(jx+2) | |
955 | svec4 = xvpow(jx+3) | |
956 | s12 = svec1 - svec2 | |
957 | s13 = svec1 - svec3 | |
958 | s23 = svec2 - svec3 | |
959 | s24 = svec2 - svec4 | |
960 | s34 = svec3 - svec4 | |
961 | sy2 = ss - svec2 | |
962 | sy3 = ss - svec3 | |
963 | const1 = s13/s23 | |
964 | const2 = s12/s23 | |
965 | const3 = s34/s23 | |
966 | const4 = s24/s23 | |
967 | s1213 = s12 + s13 | |
968 | s2434 = s24 + s34 | |
969 | sdet = s12*s34 - s1213*s2434 | |
970 | tmp = sy2*sy3/sdet | |
971 | const5 = (s34*sy2-s2434*sy3)*tmp/s12 | |
972 | const6 = (s1213*sy2-s12*sy3)*tmp/s34 | |
973 | EndIf | |
974 | ||
975 | 100 continue | |
976 | ||
977 | c skip the initialization in q if same as in the previous call. | |
978 | if(q .eq. qlast) goto 110 | |
979 | qlast = q | |
980 | ||
981 | tt = log(log(Q/Al)) | |
982 | ||
983 | JLq = -1 | |
984 | JU = NT+1 | |
985 | 12 If (JU-JLq .GT. 1) Then | |
986 | JM = (JU+JLq) / 2 | |
987 | If (Q .GE. QV(JM)) Then | |
988 | JLq = JM | |
989 | Else | |
990 | JU = JM | |
991 | Endif | |
992 | Goto 12 | |
993 | Endif | |
994 | If (JLq .LE. 0) Then | |
995 | Jq = 0 | |
996 | If (JLq .LT. 0) Then | |
997 | Msg = 'Q < Q0 in ParDis; extrapolation used!' | |
998 | CALL CtLhWARNR (IWRN2, Msg, 'Q', Q, Qini, 1D0, 1) | |
999 | EndIf | |
1000 | Elseif (JLq .LE. Nt-2) Then | |
1001 | Jq = JLq - 1 | |
1002 | Else | |
1003 | Jq = Nt - 3 | |
1004 | If (JLq .GE. Nt) Then | |
1005 | Msg = 'Q > Qmax in ParDis; extrapolation used!' | |
1006 | CALL CtLhWARNR (IWRN3, Msg, 'Q', Q, Qmax, 1D0, 1) | |
1007 | Endif | |
1008 | Endif | |
1009 | ||
1010 | If (JLq.GE.1 .and. JLq.LE.Nt-2) Then | |
1011 | tvec1 = Tv(jq) | |
1012 | tvec2 = Tv(jq+1) | |
1013 | tvec3 = Tv(jq+2) | |
1014 | tvec4 = Tv(jq+3) | |
1015 | t12 = tvec1 - tvec2 | |
1016 | t13 = tvec1 - tvec3 | |
1017 | t23 = tvec2 - tvec3 | |
1018 | t24 = tvec2 - tvec4 | |
1019 | t34 = tvec3 - tvec4 | |
1020 | ty2 = tt - tvec2 | |
1021 | ty3 = tt - tvec3 | |
1022 | tmp1 = t12 + t13 | |
1023 | tmp2 = t24 + t34 | |
1024 | tdet = t12*t34 - tmp1*tmp2 | |
1025 | EndIf | |
1026 | ||
1027 | 110 continue | |
1028 | ||
1029 | jtmp = ((IPRTN + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1 | |
1030 | Do it = 1, nqvec | |
1031 | J1 = jtmp + it*(NX+1) | |
1032 | If (Jx .Eq. 0) Then | |
1033 | fij(1) = 0 | |
1034 | fij(2) = Upd(J1+1,iset) * Xa(1,2) | |
1035 | fij(3) = Upd(J1+2,iset) * Xa(2,2) | |
1036 | fij(4) = Upd(J1+3,iset) * Xa(3,2) | |
1037 | Call CtLhPolint4 (XVpow(0), Fij(1), 4, ss, Fx, Dfx) | |
1038 | ||
1039 | If (x .GT. 0D0) Fvec(it) = Fx / x**2 | |
1040 | ElseIf (JLx .Eq. Nx-1) Then | |
1041 | Call CtLhPolint4 (XVpow(Nx-3), Upd(J1,iset), 4, ss, Fx, Dfx) | |
1042 | Fvec(it) = Fx | |
1043 | Else | |
1044 | sf2 = Upd(J1+1,iset) | |
1045 | sf3 = Upd(J1+2,iset) | |
1046 | Fvec(it) = (const5*(Upd(J1,iset) | |
1047 | & - sf2*const1 + sf3*const2) | |
1048 | & + const6*(Upd(J1+3,iset) | |
1049 | & + sf2*const3 - sf3*const4) | |
1050 | & + sf2*sy3 - sf3*sy2) / s23 | |
1051 | Endif | |
1052 | enddo | |
1053 | If (JLq .LE. 0) Then | |
1054 | Call CtLhPolint4 (TV(0), Fvec(1), 4, tt, ff, Dfq) | |
1055 | ElseIf (JLq .GE. Nt-1) Then | |
1056 | Call CtLhPolint4 (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq) | |
1057 | Else | |
1058 | tf2 = fvec(2) | |
1059 | tf3 = fvec(3) | |
1060 | g1 = ( tf2*t13 - tf3*t12) / t23 | |
1061 | g4 = (-tf2*t34 + tf3*t24) / t23 | |
1062 | h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12 | |
1063 | & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34) | |
1064 | ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23 | |
1065 | EndIf | |
1066 | CtLhPARDIS = ff | |
1067 | Return | |
1068 | End | |
1069 | ||
1070 | SUBROUTINE CtLhPARPDF (IACT, NAME, VALUE, IRET) | |
1071 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1072 | CHARACTER NAME*(*), Uname*10 | |
1073 | LOGICAL START1 | |
1074 | DATA ILEVEL, LRET / 1, 1 / | |
1075 | JRET = IRET | |
1076 | CALL CtLhUPC (NAME, Ln, Uname) | |
1077 | IF (IACT .EQ. 0 .OR. IACT .EQ. 4) then | |
1078 | c > IVALUE = NINT (VALUE) !tentatively remove this since it seems not to be used | |
1079 | print *,'Fatal error: iact=',iact | |
1080 | stop | |
1081 | ENDIF | |
1082 | START1 = (IACT .NE. 1) .AND. (IACT .NE. 2) | |
1083 | c prepare to remove this stuff, since I think IACT=1 or 2 always | |
1084 | if(start1) then | |
1085 | print *,'Fatal error: start1=',start1 | |
1086 | stop | |
1087 | endif | |
1088 | IF (START1) ILEVEL = 1 | |
1089 | GOTO (1, 2), ILEVEL | |
1090 | 1 START1 = .TRUE. | |
1091 | ILEVEL = 0 | |
1092 | CALL CtLhParQcd (IACT, Uname(1:Ln), VALUE, JRET) | |
1093 | IF (JRET .EQ. 1) GOTO 11 | |
1094 | IF (JRET .EQ. 2) GOTO 12 | |
1095 | IF (JRET .EQ. 3) GOTO 13 | |
1096 | IF (JRET .GT. 4) GOTO 15 | |
1097 | ILEVEL = ILEVEL + 1 | |
1098 | 2 CALL CtLhEVLPAR (IACT, Uname(1:Ln), VALUE, JRET) | |
1099 | IF (JRET .EQ. 1) GOTO 11 | |
1100 | IF (JRET .EQ. 2) GOTO 12 | |
1101 | IF (JRET .EQ. 3) GOTO 13 | |
1102 | IF (JRET .GT. 4) GOTO 15 | |
1103 | ILEVEL = ILEVEL + 1 | |
1104 | IF (.NOT. START1) GOTO 1 | |
1105 | IF (JRET .EQ. 0) GOTO 10 | |
1106 | GOTO 14 | |
1107 | 10 CONTINUE | |
1108 | 11 CONTINUE | |
1109 | 12 CONTINUE | |
1110 | 13 CONTINUE | |
1111 | 14 CONTINUE | |
1112 | 15 CONTINUE | |
1113 | IF (JRET .NE. 4) LRET = JRET | |
1114 | IF (LRET.EQ.0 .OR. LRET.EQ.2 .OR. LRET.EQ.3) THEN | |
1115 | PRINT *, 'Error in CtLhPARPDF: IRET, IACT, NAME, VALUE =', | |
1116 | > LRET, IACT, NAME, VALUE | |
1117 | PRINT *, 'fatal error in CtLhparpdf' | |
1118 | stop | |
1119 | EndIf | |
1120 | IRET= JRET | |
1121 | RETURN | |
1122 | END | |
1123 | SUBROUTINE CtLhParQcd(IACT,NAME,VALUE,IRET) | |
1124 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1125 | INTEGER IACT,IRET | |
1126 | CHARACTER*(*) NAME | |
1127 | IRET=1 | |
1128 | IF (IACT.EQ.0) THEN | |
1129 | WRITE (NINT(VALUE), *) 'LAM(BDA), NFL, ORD(ER), Mi, ', | |
1130 | > '(i in 1 to 9), LAMi (i in 1 to NFL)' | |
1131 | IRET=4 | |
1132 | ELSEIF (IACT.EQ.1) THEN | |
1133 | CALL CtLhQCDSET (NAME,VALUE,IRET) | |
1134 | ELSEIF (IACT.EQ.2) THEN | |
1135 | CALL CtLhQCDGET (NAME,VALUE,IRET) | |
1136 | ELSE | |
1137 | IRET=3 | |
1138 | ENDIF | |
1139 | RETURN | |
1140 | END | |
1141 | FUNCTION CtLhPFF1 (XX) | |
1142 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1143 | LOGICAL LA, LB, LSTX | |
1144 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
1145 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
1146 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
1147 | PARAMETER (MX = 3) | |
1148 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
1149 | COMMON / LhCtKRNL00 / DZ, XL(MX), NNX | |
1150 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
1151 | COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX), | |
1152 | > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX), | |
1153 | > PNSP(0:MXX), PNSM(0:MXX) | |
1154 | SAVE | |
1155 | DATA LA, LB / 2 * .FALSE. / | |
1156 | LB = .TRUE. | |
1157 | ENTRY CtLhTFF1(ZZ) | |
1158 | LA = .TRUE. | |
1159 | IF (LA .AND. .NOT.LB) THEN | |
1160 | Z = ZZ | |
1161 | X = CtLhXFRMZ (Z) | |
1162 | Else | |
1163 | X = XX | |
1164 | EndIf | |
1165 | IF (X .GE. D1) THEN | |
1166 | CtLhPFF1 = 0 | |
1167 | RETURN | |
1168 | ElseIF (X .GE. XMIN) THEN | |
1169 | Z = CtLhZFRMX (X) | |
1170 | TEM = CtLhFINTRP (FF1, -DZ, DZ, NX, Z, ERR, IRT) | |
1171 | Else | |
1172 | CALL CtLhPOLIN1 (XL, FF1(1), MX, X, TEM, ERR) | |
1173 | EndIf | |
1174 | IF (LA) THEN | |
1175 | IF (LB) THEN | |
1176 | CtLhPFF1 = TEM / (1.-X) | |
1177 | LB =.FALSE. | |
1178 | Else | |
1179 | CtLhTFF1 = TEM / (1.-X) * CtLhDXDZ(Z) | |
1180 | EndIf | |
1181 | LA =.FALSE. | |
1182 | Else | |
1183 | IF (LB) THEN | |
1184 | c +++ something is wrong, since QFF1 and RFF1 are not used. | |
1185 | c +++ but this code appears to only be used for extrapolation | |
1186 | c +++ to small x, which is unreliable anyway, so ignore for now (jcp) | |
1187 | QFF1 = TEM | |
1188 | LB =.FALSE. | |
1189 | Else | |
1190 | RFF1 = TEM * X / (1.-X) | |
1191 | EndIf | |
1192 | EndIf | |
1193 | RETURN | |
1194 | ENTRY CtLhFNSP (XX) | |
1195 | X = XX | |
1196 | IF (X .GE. D1) THEN | |
1197 | CtLhFNSP = 0. | |
1198 | RETURN | |
1199 | ElseIF (X .GE. XMIN) THEN | |
1200 | Z = CtLhZFRMX (X) | |
1201 | TEM = CtLhFINTRP (PNSP, -DZ, DZ, NX, Z, ERR, IRT) | |
1202 | Else | |
1203 | CALL CtLhPOLIN1 (XL, PNSP(1), MX, X, TEM, ERR) | |
1204 | EndIf | |
1205 | CtLhFNSP = TEM / (1.- X) | |
1206 | RETURN | |
1207 | ENTRY CtLhFNSM (XX) | |
1208 | X = XX | |
1209 | IF (X .GE. D1) THEN | |
1210 | CtLhFNSM = 0. | |
1211 | RETURN | |
1212 | ElseIF (X .GE. XMIN) THEN | |
1213 | Z = CtLhZFRMX (X) | |
1214 | TEM = CtLhFINTRP (PNSM, -DZ, DZ, NX, Z, ERR, IRT) | |
1215 | Else | |
1216 | CALL CtLhPOLIN1 (XL, PNSM(1), MX, X, TEM, ERR) | |
1217 | EndIf | |
1218 | CtLhFNSM = TEM / (1.- X) | |
1219 | RETURN | |
1220 | ENTRY CtLhRGG1 (XX) | |
1221 | X = XX | |
1222 | IF (X .GE. D1) THEN | |
1223 | CtLhRGG1= 0 !error corrected? (jcp) | |
1224 | RETURN | |
1225 | ElseIF (X .GE. XMIN) THEN | |
1226 | Z = CtLhZFRMX (X) | |
1227 | TEM = CtLhFINTRP (GG1, -DZ, DZ, NX, Z, ERR, IRT) | |
1228 | Else | |
1229 | CALL CtLhPOLIN1 (XL, GG1(1), MX, X, TEM, ERR) | |
1230 | EndIf | |
1231 | IF (LA) THEN | |
1232 | CtLhRGG1 = TEM / X / (1.-X) !error corrected? (jcp) | |
1233 | LA =.FALSE. | |
1234 | Else | |
1235 | IF (LB) THEN | |
1236 | CtLhRGG1 = TEM / X !error corrected? (jcp) | |
1237 | LB =.FALSE. | |
1238 | Else | |
1239 | CtLhRGG1 = TEM / (1.-X) | |
1240 | EndIf | |
1241 | EndIf | |
1242 | RETURN | |
1243 | ENTRY CtLhRFF2 (XX) | |
1244 | X = XX | |
1245 | IF (X .GE. D1) THEN | |
1246 | CtLhRFF2 = 0 !error corrected? (jcp) | |
1247 | RETURN | |
1248 | ElseIF (X .GE. XMIN) THEN | |
1249 | Z = CtLhZFRMX (X) | |
1250 | TEM = CtLhFINTRP (FF2, -DZ, DZ, NX, Z, ERR, IRT) | |
1251 | Else | |
1252 | CALL CtLhPOLIN1 (XL, FF2(1), MX, X, TEM, ERR) | |
1253 | EndIf | |
1254 | IF (LA) THEN | |
1255 | CtLhRFF2 = TEM / X / (1.-X) !error corrected? (jcp) | |
1256 | LA =.FALSE. | |
1257 | Else | |
1258 | IF (LB) THEN | |
1259 | CtLhRFF2 = TEM / X !error corrected? (jcp) | |
1260 | LB =.FALSE. | |
1261 | Else | |
1262 | CtLhRFF2 = TEM / (1.-X) | |
1263 | EndIf | |
1264 | EndIf | |
1265 | RETURN | |
1266 | ENTRY CtLhRGG2 (XX) | |
1267 | X = XX | |
1268 | IF (X .GE. D1) THEN | |
1269 | CtLhRGG2 = 0 !error corrected? (jcp) | |
1270 | RETURN | |
1271 | ElseIF (X .GE. XMIN) THEN | |
1272 | Z = CtLhZFRMX (X) | |
1273 | TEM = CtLhFINTRP (GG2, -DZ, DZ, NX, Z, ERR, IRT) | |
1274 | Else | |
1275 | CALL CtLhPOLIN1 (XL, GG2(1), MX, X, TEM, ERR) | |
1276 | EndIf | |
1277 | IF (LA) THEN | |
1278 | CtLhRGG2 = TEM / X / (1.-X) !error corrected? (jcp) | |
1279 | LA =.FALSE. | |
1280 | Else | |
1281 | IF (LB) THEN | |
1282 | CtLhRGG2 = TEM / X !error corrected? (jcp) | |
1283 | LB =.FALSE. | |
1284 | Else | |
1285 | CtLhRGG2 = TEM / (1.-X) | |
1286 | EndIf | |
1287 | EndIf | |
1288 | RETURN | |
1289 | END | |
1290 | SUBROUTINE CtLhPOLIN1 (XA,YA,N,X,Y,DY) | |
1291 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1292 | PARAMETER (NMAX=10) | |
1293 | DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) | |
1294 | NS=1 | |
1295 | DIF=ABS(X-XA(1)) | |
1296 | DO 11 I=1,N | |
1297 | DIFT=ABS(X-XA(I)) | |
1298 | IF (DIFT.LT.DIF) THEN | |
1299 | NS=I | |
1300 | DIF=DIFT | |
1301 | ENDIF | |
1302 | C(I)=YA(I) | |
1303 | D(I)=YA(I) | |
1304 | 11 CONTINUE | |
1305 | Y=YA(NS) | |
1306 | NS=NS-1 | |
1307 | DO 13 M=1,N-1 | |
1308 | DO 12 I=1,N-M | |
1309 | HO=XA(I)-X | |
1310 | HP=XA(I+M)-X | |
1311 | W=C(I+1)-D(I) | |
1312 | DEN=HO-HP | |
1313 | DEN=W/DEN | |
1314 | D(I)=HP*DEN | |
1315 | C(I)=HO*DEN | |
1316 | 12 CONTINUE | |
1317 | IF (2*NS.LT.N-M)THEN | |
1318 | DY=C(NS+1) | |
1319 | ELSE | |
1320 | DY=D(NS) | |
1321 | NS=NS-1 | |
1322 | ENDIF | |
1323 | Y=Y+DY | |
1324 | 13 CONTINUE | |
1325 | RETURN | |
1326 | END | |
1327 | SUBROUTINE CtLhQARRAY (NINI) | |
1328 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1329 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
1330 | PARAMETER (MXPN = MXF * 2 + 2) | |
1331 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
1332 | COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG | |
1333 | COMMON / LhCtQARAY2 / TLN(MXF), DTN(MXF), NTL(MXF), NTN(MXF) | |
1334 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
1335 | NCNT = 0 | |
1336 | IF (NT .GE. mxq) NT = mxq - 1 | |
1337 | S = LOG(QINI/AL) | |
1338 | TINI = LOG(S) | |
1339 | S = LOG(QMAX/AL) | |
1340 | TMAX = LOG(S) | |
1341 | 1 DT0 = (TMAX - TINI) / float(NT) | |
1342 | NINI = LhCtNFL(QINI) | |
1343 | NFMX = LhCtNFL(QMAX) | |
1344 | Call CtLhParQcd (2, 'ORDER', Ord, Ir) | |
1345 | Call CtLhParQcd (2, 'ALAM', Al0, Ir) | |
1346 | Call CtLhParQcd (2, 'NFL', Afl0, Ir) | |
1347 | AFL = NfMx | |
1348 | Call CtLhParQcd (1, 'NFL', AFL, Ir) | |
1349 | Iordr = Nint (Ord) | |
1350 | Ifl0 = Nint (Afl0) | |
1351 | Call CtLhSetLam (Ifl0, Al0, Iordr) | |
1352 | NG = NFMX - NINI + 1 | |
1353 | QIN = QINI | |
1354 | QOUT = QINI | |
1355 | S = LOG(QIN/AL) | |
1356 | TIN = LOG(S) | |
1357 | TLN(1) = TIN | |
1358 | NTL(1) = 0 | |
1359 | QV(0) = QINI | |
1360 | TV(0) = Tin | |
1361 | DO 20 NEFF = NINI, NFMX | |
1362 | ICNT = NEFF - NINI + 1 | |
1363 | IF (NEFF .LT. NFMX) THEN | |
1364 | THRN = CtLhAMHATF (NEFF + 1) | |
1365 | QOUN = MIN (QMAX, THRN) | |
1366 | Else | |
1367 | QOUN = QMAX | |
1368 | EndIf | |
1369 | IF (QOUN-QOUT .LE. 0.0001) THEN | |
1370 | DT = 0 | |
1371 | NITR = 0 | |
1372 | Else | |
1373 | QOUT = QOUN | |
1374 | S = LOG(QOUT/AL) | |
1375 | TOUT = LOG(S) | |
1376 | TEM = TOUT - TIN | |
1377 | NITR = INT (TEM / DT0) + 1 | |
1378 | DT = TEM / NITR | |
1379 | EndIf | |
1380 | DTN (ICNT) = DT | |
1381 | NTN (ICNT) = NITR | |
1382 | TLN (ICNT) = TIN | |
1383 | NTL (ICNT+1) = NTL(ICNT) + NITR | |
1384 | IF (NITR .NE. 0) THEN | |
1385 | DO 205 I = 1, NITR | |
1386 | TV (NTL(ICNT)+I) = TIN + DT * I | |
1387 | S = EXP (TV(NTL(ICNT)+I)) | |
1388 | QV (NTL(ICNT)+I) = AL * EXP (S) | |
1389 | 205 CONTINUE | |
1390 | EndIf | |
1391 | QIN = QOUT | |
1392 | TIN = TOUT | |
1393 | 20 CONTINUE | |
1394 | NCNT = NCNT + 1 | |
1395 | NTP = NTL (NG + 1) | |
1396 | ND = NTP - NT | |
1397 | IF (NTP .GE. MXQ) THEN | |
1398 | NT = MXQ - ND - NCNT | |
1399 | GOTO 1 | |
1400 | EndIf | |
1401 | NT = NTP | |
1402 | RETURN | |
1403 | END | |
1404 | SUBROUTINE CtLhQCDGET(NAME,VALUE,IRET) | |
1405 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1406 | CHARACTER*(*) NAME | |
1407 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
1408 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
1409 | COMMON / LhCtCOMQMS / VALQMS(9) | |
1410 | LOGICAL SET | |
1411 | PARAMETER (PI=3.141592653589793d0, EULER=0.57721566) | |
1412 | ICODE = LhCtNAMQCD(NAME) | |
1413 | IRET = 1 | |
1414 | IF (ICODE .EQ. 1) THEN | |
1415 | VALUE = AL | |
1416 | ELSEIF (ICODE .EQ. 2) THEN | |
1417 | VALUE = NF | |
1418 | ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 12)) THEN | |
1419 | VALUE = VALQMS(ICODE - 2) | |
1420 | ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN | |
1421 | VALUE = ALAM(ICODE - 13) | |
1422 | ELSEIF (ICODE .EQ. 24) THEN | |
1423 | VALUE = NORDER | |
1424 | ELSE | |
1425 | IRET=0 | |
1426 | ENDIF | |
1427 | END | |
1428 | SUBROUTINE CtLhQCDSET (NAME,VALUE,IRET) | |
1429 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1430 | CHARACTER*(*) NAME | |
1431 | COMMON / LhCtCOMQMS / VALQMS(9) | |
1432 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
1433 | LOGICAL SET | |
1434 | PARAMETER (PI=3.141592653589793d0, EULER=0.57721566) | |
1435 | IVALUE = NINT(VALUE) | |
1436 | ICODE = LhCtNAMQCD(NAME) | |
1437 | IF (ICODE .EQ. 0) THEN | |
1438 | IRET=0 | |
1439 | c print *,'warning empty CtLhQCDSET call: NAME=', | |
1440 | c & NAME,' VALUE=',VALUE | |
1441 | ELSE | |
1442 | IRET = 1 | |
1443 | SET = .FALSE. | |
1444 | IF (ICODE .EQ. 1) THEN | |
1445 | IF (VALUE.LE.0) GOTO 12 | |
1446 | AL=VALUE | |
1447 | ELSEIF (ICODE .EQ. 2) THEN | |
1448 | IF ( (IVALUE .LT. 0) .OR. (IVALUE .GT. 9)) GOTO 12 | |
1449 | NF = IVALUE | |
1450 | ELSEIF ((ICODE .GE. 3) .AND. (ICODE .LE. 11)) THEN | |
1451 | IF (VALUE .LT. 0) GOTO 12 | |
1452 | Scle = Min (Value , VALQMS(ICODE - 2)) | |
1453 | AlfScle = CtLhALPI(Scle) * Pi | |
1454 | VALQMS(ICODE - 2) = VALUE | |
1455 | Call CtLhAlfSet (Scle, AlfScle) | |
1456 | ELSEIF ((ICODE .GE. 13) .AND. (ICODE .LE. 13+NF)) THEN | |
1457 | IF (VALUE .LE. 0) GOTO 12 | |
1458 | CALL CtLhSETL1 (ICODE-13, VALUE) | |
1459 | ELSEIF (ICODE .EQ. 24) THEN | |
1460 | IF ((IVALUE .LT. 1) .OR. (IVALUE .GT. 2)) GOTO 12 | |
1461 | NORDER = IVALUE | |
1462 | ENDIF | |
1463 | IF (.NOT. SET) CALL CtLhLAMCWZ | |
1464 | ENDIF | |
1465 | RETURN | |
1466 | 12 IRET=2 | |
1467 | RETURN | |
1468 | END | |
1469 | FUNCTION CtLhQZBRNT(FUNC, X1, X2, TOLIN, IRT) | |
1470 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1471 | PARAMETER (ITMAX = 1000, EPS = 3.E-12) | |
1472 | external func | |
1473 | TOL = ABS(TOLIN) | |
1474 | A=X1 | |
1475 | B=X2 | |
1476 | FA=FUNC(A) | |
1477 | FB=FUNC(B) | |
1478 | IF(FB*FA.GT.0.) THEN | |
1479 | WRITE (*, *) 'Root must be bracketed for CtLhQZBRNT.' | |
1480 | IRT = 1 | |
1481 | ENDIF | |
1482 | FC=FB | |
1483 | DO 11 ITER=1,ITMAX | |
1484 | IF(FB*FC.GT.0.) THEN | |
1485 | C=A | |
1486 | FC=FA | |
1487 | D=B-A | |
1488 | E=D | |
1489 | ENDIF | |
1490 | IF(ABS(FC).LT.ABS(FB)) THEN | |
1491 | A=B | |
1492 | B=C | |
1493 | C=A | |
1494 | FA=FB | |
1495 | FB=FC | |
1496 | FC=FA | |
1497 | ENDIF | |
1498 | TOL1=2.*EPS*ABS(B)+0.5*TOL | |
1499 | XM=.5*(C-B) | |
1500 | IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN | |
1501 | CtLhQZBRNT=B | |
1502 | RETURN | |
1503 | ENDIF | |
1504 | IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN | |
1505 | S=FB/FA | |
1506 | IF(A.EQ.C) THEN | |
1507 | P=2.*XM*S | |
1508 | Q=1.-S | |
1509 | ELSE | |
1510 | Q=FA/FC | |
1511 | R=FB/FC | |
1512 | P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.)) | |
1513 | Q=(Q-1.)*(R-1.)*(S-1.) | |
1514 | ENDIF | |
1515 | IF(P.GT.0.) Q=-Q | |
1516 | P=ABS(P) | |
1517 | IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN | |
1518 | E=D | |
1519 | D=P/Q | |
1520 | ELSE | |
1521 | D=XM | |
1522 | E=D | |
1523 | ENDIF | |
1524 | ELSE | |
1525 | D=XM | |
1526 | E=D | |
1527 | ENDIF | |
1528 | A=B | |
1529 | FA=FB | |
1530 | IF(ABS(D) .GT. TOL1) THEN | |
1531 | B=B+D | |
1532 | ELSE | |
1533 | B=B+SIGN(TOL1,XM) | |
1534 | ENDIF | |
1535 | FB=FUNC(B) | |
1536 | 11 CONTINUE | |
1537 | WRITE (*, *) 'CtLhQZBRNT exceeding maximum iterations.' | |
1538 | IRT = 2 | |
1539 | CtLhQZBRNT=B | |
1540 | RETURN | |
1541 | END | |
1542 | SUBROUTINE CtLhRATINT(XA,YA,N,X,Y,DY) | |
1543 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1544 | PARAMETER (NMAX=10,TINY=1.E-25) | |
1545 | DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) | |
1546 | NS=1 | |
1547 | HH=ABS(X-XA(1)) | |
1548 | DO 11 I=1,N | |
1549 | H=ABS(X-XA(I)) | |
1550 | IF (H.EQ.0.)THEN | |
1551 | Y=YA(I) | |
1552 | DY=0.0 | |
1553 | RETURN | |
1554 | ELSE IF (H.LT.HH) THEN | |
1555 | NS=I | |
1556 | HH=H | |
1557 | ENDIF | |
1558 | C(I)=YA(I) | |
1559 | D(I)=YA(I)+TINY | |
1560 | 11 CONTINUE | |
1561 | Y=YA(NS) | |
1562 | NS=NS-1 | |
1563 | DO 13 M=1,N-1 | |
1564 | DO 12 I=1,N-M | |
1565 | W=C(I+1)-D(I) | |
1566 | H=XA(I+M)-X | |
1567 | T=(XA(I)-X)*D(I)/H | |
1568 | DD=T-C(I+1) | |
1569 | DD=W/DD | |
1570 | D(I)=C(I+1)*DD | |
1571 | C(I)=T*DD | |
1572 | 12 CONTINUE | |
1573 | IF (2*NS.LT.N-M)THEN | |
1574 | DY=C(NS+1) | |
1575 | ELSE | |
1576 | DY=D(NS) | |
1577 | NS=NS-1 | |
1578 | ENDIF | |
1579 | Y=Y+DY | |
1580 | 13 CONTINUE | |
1581 | RETURN | |
1582 | END | |
1583 | FUNCTION CtLhRTALF (EFLLN) | |
1584 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1585 | include 'parmsetup.inc' | |
1586 | PARAMETER (PI = 3.141592653589793d0) | |
1587 | COMMON / CtLhRTALFC / ALFST, JORD, NEFF | |
1588 | EFMULM = EXP (EFLLN) | |
1589 | TEM1 = PI / ALFST | |
1590 | TEM2 = 1. / CtLhALPQCD (JORD, NEFF, EFMULM, I) | |
1591 | CtLhRTALF = TEM1 - TEM2 | |
1592 | END | |
1593 | Subroutine CtLhbldat1 | |
1594 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1595 | include 'parmsetup.inc' | |
1596 | LOGICAL LSTX | |
1597 | PARAMETER (MXX = 105, MXQ = 25, MxF = 6) | |
1598 | PARAMETER (MxPN = MxF * 2 + 2) | |
1599 | PARAMETER (MxQX= MXQ * MXX, MxPQX = MxQX * MxPN) | |
1600 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
1601 | COMMON / LhCtQARAY1 / QINI,QMAX, QV(0:MXQ),TV(0:MXQ), NT,JT,NG | |
1602 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
1603 | COMMON / LhCtPEVLDT / UPD(MXPQX,nmxset), KF, Nelmt | |
1604 | PARAMETER (NF0 = 4, Nshp = 8,NEX = Nshp+2) | |
1605 | XMIN = .999999D-4 | |
1606 | XCR = 1.5 | |
1607 | JT = 1 | |
1608 | Return | |
1609 | END | |
1610 | SUBROUTINE CtLhSETL1 (NEF, VLAM) | |
1611 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1612 | LOGICAL SET | |
1613 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
1614 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
1615 | COMMON / LhCtCOMQMS / VALQMS(9) | |
1616 | IF (NEF .LT. 0 .OR. NEF .GT. NF) THEN | |
1617 | WRITE(*,*)'NEF out of range in CtLhSETL1: NEF NF =',NEF,NF | |
1618 | STOP | |
1619 | ENDIF | |
1620 | AMHAT(0) = 0. | |
1621 | DO 5 N = 1, NF | |
1622 | AMHAT(N) = VALQMS(N) | |
1623 | 5 CONTINUE | |
1624 | ALAM(NEF) = VLAM | |
1625 | DO 10 N = NEF, 1, -1 | |
1626 | CALL CtLhTRNLAM(NORDER, N, -1, IR1) | |
1627 | 10 CONTINUE | |
1628 | DO 20 N = NEF, NF-1 | |
1629 | CALL CtLhTRNLAM(NORDER, N, 1, IR1) | |
1630 | 20 CONTINUE | |
1631 | DO 30, N = NF, 1, -1 | |
1632 | IF ((ALAM(N) .GE. 0.7 * AMHAT(N)) | |
1633 | > .OR. (ALAM(N-1) .GE. 0.7 * AMHAT(N)))THEN | |
1634 | NHQ = NF - N | |
1635 | GOTO 40 | |
1636 | ENDIF | |
1637 | 30 CONTINUE | |
1638 | NHQ = NF | |
1639 | 40 CONTINUE | |
1640 | DO 50, N = NF-NHQ, 1, -1 | |
1641 | AMHAT(N) = 0 | |
1642 | ALAM(N-1) = ALAM(N) | |
1643 | 50 CONTINUE | |
1644 | AMN = ALAM(NF) | |
1645 | DO 60, N = 0, NF-1 | |
1646 | IF (ALAM(N) .GT. AMN) AMN = ALAM(N) | |
1647 | 60 CONTINUE | |
1648 | AMN = AMN * 1.0001 | |
1649 | AL = ALAM(NF) | |
1650 | SET = .TRUE. | |
1651 | RETURN | |
1652 | END | |
1653 | SUBROUTINE CtLhSETLAM (NEF, WLAM, IRDR) | |
1654 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1655 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
1656 | LOGICAL SET | |
1657 | IF ((NEF .LT. 0) .OR. (NEF .GT. NF)) THEN | |
1658 | WRITE(*,*)'NEF out of range in CtLhSETLAM: NEF NF=',NEF,NF | |
1659 | STOP | |
1660 | ENDIF | |
1661 | VLAM = WLAM | |
1662 | IF (IRDR .NE. NORDER) then | |
1663 | PRINT *,'fatal error: wanted cnvl1' | |
1664 | stop | |
1665 | ENDIF | |
1666 | CALL CtLhSETL1 (NEF, VLAM) | |
1667 | END | |
1668 | Subroutine CtLhbldat2 | |
1669 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1670 | COMMON / LhCtCOMQMS / VALQMS(9) | |
1671 | COMMON / LhCtQCDPAR_LHA / AL, NF, NORDER, SET | |
1672 | LOGICAL SET | |
1673 | AL = .226d0 | |
1674 | NF = 5 | |
1675 | NORDER = 2 | |
1676 | SET = .FALSE. | |
1677 | VALQMS(1) = 0. | |
1678 | VALQMS(2) = 0. | |
1679 | VALQMS(3) = 0.2d0 | |
1680 | VALQMS(4) = 1.3d0 | |
1681 | VALQMS(5) = 4.5d0 | |
1682 | VALQMS(6) = 174.d0 | |
1683 | VALQMS(7) = 0. | |
1684 | VALQMS(8) = 0. | |
1685 | VALQMS(9) = 0. | |
1686 | Return | |
1687 | END | |
1688 | FUNCTION CtLhSMPNOL (NX, DX, FN, ERR) | |
1689 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1690 | DIMENSION FN(NX) | |
1691 | MS = MOD(NX, 2) | |
1692 | IF (NX .LE. 1 .OR. NX .GT. 1000) THEN | |
1693 | PRINT *, 'NX =', NX, ' OUT OF RANGE IN CtLhSMPNOL!' | |
1694 | STOP | |
1695 | ELSEIF (NX .EQ. 2) THEN | |
1696 | TEM = DX * FN(2) | |
1697 | ELSEIF (NX .EQ. 3) THEN | |
1698 | TEM = DX * FN(2) * 2. | |
1699 | ELSE | |
1700 | IF (MS .EQ. 0) THEN | |
1701 | TEM = DX * (23.* FN(2) - 16.* FN(3) + 5.* FN(4)) / 12. | |
1702 | TMP = DX * (3.* FN(2) - FN(3)) / 2. | |
1703 | ERR = ABS(TEM - TMP) | |
1704 | TEM = TEM + CtLhSMPSNA (NX-1, DX, FN(2), ER1) | |
1705 | ERR = ABS(ER1) + ERR | |
1706 | ELSE | |
1707 | TEM = DX * (8.* FN(2) - 4.* FN(3) + 8.* FN(4)) / 3. | |
1708 | TMP = DX * (3.* FN(2) + 2.* FN(3) + 3.* FN(4)) / 2. | |
1709 | ERR = ABS(TEM - TMP) | |
1710 | TEM = TEM + CtLhSMPSNA (NX-4, DX, FN(5), ER1) | |
1711 | ERR = ABS(ER1) + ERR | |
1712 | ENDIF | |
1713 | ENDIF | |
1714 | CtLhSMPNOL = TEM | |
1715 | RETURN | |
1716 | END | |
1717 | FUNCTION CtLhSMPSNA (NX, DX, F, ERR) | |
1718 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1719 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
1720 | PARAMETER (MAXX = 1000) | |
1721 | DIMENSION F(NX) | |
1722 | DATA IW1, IW2, TINY / 2*0, 1.E-35 / | |
1723 | IF (DX .LE. 0.) THEN | |
1724 | CALL CtLhWARNR(IW2,'DX cannot be < 0. in CtLhSMPSNA', 'DX', | |
1725 | > DX, D0, D1, 0) | |
1726 | CtLhSMPSNA = 0. | |
1727 | RETURN | |
1728 | ENDIF | |
1729 | IF (NX .LE. 0 .OR. NX .GT. MAXX) THEN | |
1730 | CALL CtLhWARNI(IW1, 'NX out of range in CtLhSMPSNA', 'NX', NX, | |
1731 | > 1, MAXX, 1) | |
1732 | SIMP = 0. | |
1733 | ELSEIF (NX .EQ. 1) THEN | |
1734 | SIMP = 0. | |
1735 | ELSEIF (NX .EQ. 2) THEN | |
1736 | SIMP = (F(1) + F(2)) / 2. | |
1737 | ERRD = (F(1) - F(2)) / 2. | |
1738 | ELSE | |
1739 | MS = MOD(NX, 2) | |
1740 | IF (MS .EQ. 0) THEN | |
1741 | ADD = (9.*F(NX) + 19.*F(NX-1) - 5.*F(NX-2) + F(NX-3)) / 24. | |
1742 | NZ = NX - 1 | |
1743 | ELSE | |
1744 | ADD = 0. | |
1745 | NZ = NX | |
1746 | ENDIF | |
1747 | IF (NZ .EQ. 3) THEN | |
1748 | SIMP = (F(1) + 4.* F(2) + F(3)) / 3. | |
1749 | TRPZ = (F(1) + 2.* F(2) + F(3)) / 2. | |
1750 | ELSE | |
1751 | SE = F(2) | |
1752 | SO = 0 | |
1753 | NM1 = NZ - 1 | |
1754 | DO 60 I = 4, NM1, 2 | |
1755 | IM1 = I - 1 | |
1756 | SE = SE + F(I) | |
1757 | SO = SO + F(IM1) | |
1758 | 60 CONTINUE | |
1759 | SIMP = (F(1) + 4.* SE + 2.* SO + F(NZ)) / 3. | |
1760 | TRPZ = (F(1) + 2.* (SE + SO) + F(NZ)) / 2. | |
1761 | ENDIF | |
1762 | ERRD = TRPZ - SIMP | |
1763 | SIMP = SIMP + ADD | |
1764 | ENDIF | |
1765 | CtLhSMPSNA = SIMP * DX | |
1766 | IF (ABS(SIMP) .GT. TINY) THEN | |
1767 | ERR = ERRD / SIMP | |
1768 | ELSE | |
1769 | ERR = 0. | |
1770 | ENDIF | |
1771 | RETURN | |
1772 | END | |
1773 | SUBROUTINE CtLhSNEVL(IKNL,NX,NT,JT,DT,TIN,NEFF,UI,GI,US,GS) | |
1774 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1775 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
1776 | PARAMETER (MXQX= MXQ * MXX) | |
1777 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
1778 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
1779 | DIMENSION UI(NX), US(0:NX, 0:NT) | |
1780 | DIMENSION GI(NX), GS(0:NX, 0:NT) | |
1781 | DIMENSION Y0(MXX), Y1(MXX), YP(MXX), F0(MXX), F1(MXX), FP(MXX) | |
1782 | DIMENSION Z0(MXX), Z1(MXX), ZP(MXX), G0(MXX), G1(MXX), GP(MXX) | |
1783 | DATA D0 / 0.0 / | |
1784 | JTT = 2 * JT | |
1785 | DDT = DT / JTT | |
1786 | IF (NX .GT. MXX) THEN | |
1787 | WRITE (*,*) 'Nx =', NX, ' too many pts in CtLhSNEVL' | |
1788 | STOP 'Program stopped in CtLhSNEVL' | |
1789 | EndIf | |
1790 | c ++ remove unused quantities (jcp) | |
1791 | c ++ TMD = TIN + DT * NT / 2. | |
1792 | c ++ AMU = EXP(TMD) | |
1793 | c ++ TEM = 6./ (33.- 2.* NEFF) / CtLhALPI(AMU) | |
1794 | c ++ TLAM = TMD - TEM | |
1795 | DO 9 IX = 1, NX | |
1796 | US (IX, 0) = UI(IX) | |
1797 | GS (IX, 0) = GI(IX) | |
1798 | 9 CONTINUE | |
1799 | US ( 0, 0) = (UI(1) - UI(2))* 3D0 + UI(3) | |
1800 | GS ( 0, 0) = (GI(1) - GI(2))* 3D0 + GI(3) | |
1801 | TT = TIN | |
1802 | DO 10 IZ = 1, NX | |
1803 | Y0(IZ) = UI(IZ) | |
1804 | Z0(IZ) = GI(IZ) | |
1805 | 10 CONTINUE | |
1806 | DO 20 IS = 1, NT | |
1807 | DO 202 JS = 1, JTT | |
1808 | IRND = (IS-1) * JTT + JS | |
1809 | IF (IRND .EQ. 1) THEN | |
1810 | CALL CtLhSNRHS (TT, NEFF, Y0,Z0, F0,G0) | |
1811 | DO 250 IZ = 1, NX | |
1812 | Y0(IZ) = Y0(IZ) + DDT * F0(IZ) | |
1813 | Z0(IZ) = Z0(IZ) + DDT * G0(IZ) | |
1814 | 250 CONTINUE | |
1815 | TT = TT + DDT | |
1816 | CALL CtLhSNRHS (TT, NEFF, Y0, Z0, F1, G1) | |
1817 | DO 251 IZ = 1, NX | |
1818 | Y1(IZ) = UI(IZ) + DDT * (F0(IZ) + F1(IZ)) / 2D0 | |
1819 | Z1(IZ) = GI(IZ) + DDT * (G0(IZ) + G1(IZ)) / 2D0 | |
1820 | 251 CONTINUE | |
1821 | Else | |
1822 | CALL CtLhSNRHS (TT, NEFF, Y1, Z1, F1, G1) | |
1823 | DO 252 IZ = 1, NX | |
1824 | YP(IZ) = Y1(IZ) + DDT * (3D0 * F1(IZ) - F0(IZ)) / 2D0 | |
1825 | ZP(IZ) = Z1(IZ) + DDT * (3D0 * G1(IZ) - G0(IZ)) / 2D0 | |
1826 | 252 CONTINUE | |
1827 | TT = TT + DDT | |
1828 | CALL CtLhSNRHS (TT, NEFF, YP, ZP, FP, GP) | |
1829 | DO 253 IZ = 1, NX | |
1830 | Y1(IZ) = Y1(IZ) + DDT * (FP(IZ) + F1(IZ)) / 2D0 | |
1831 | Z1(IZ) = Z1(IZ) + DDT * (GP(IZ) + G1(IZ)) / 2D0 | |
1832 | F0(IZ) = F1(IZ) | |
1833 | G0(IZ) = G1(IZ) | |
1834 | 253 CONTINUE | |
1835 | EndIf | |
1836 | 202 CONTINUE | |
1837 | DO 260 IX = 1, NX | |
1838 | IF (IKNL .GT. 0) THEN | |
1839 | US (IX, IS) = MAX(Y1(IX), D0) | |
1840 | GS (IX, IS) = MAX(Z1(IX), D0) | |
1841 | Else | |
1842 | US (IX, IS) = Y1(IX) | |
1843 | GS (IX, IS) = Z1(IX) | |
1844 | EndIf | |
1845 | 260 CONTINUE | |
1846 | US(0, IS) = 3D0*Y1(1) - 3D0*Y1(2) + Y1(3) | |
1847 | GS(0, IS) = 3D0*Z1(1) - 3D0*Z1(2) + Z1(3) | |
1848 | 20 CONTINUE | |
1849 | RETURN | |
1850 | END | |
1851 | SUBROUTINE CtLhSNRHS (TT, NEFF, FI, GI, FO, GO) | |
1852 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1853 | LOGICAL LSTX | |
1854 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
1855 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
1856 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
1857 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
1858 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
1859 | COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX), | |
1860 | > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB | |
1861 | COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) | |
1862 | COMMON / LhCtEVLPAC / AL, IKNL, IPD0, IHDN, NfMx | |
1863 | DIMENSION GI(NX), GO(NX), G1(MXX), G2(MXX), G3(MXX), G4(MXX) | |
1864 | DIMENSION FI(NX), FO(NX), W0(MXX), W1(MXX), WH(MXX), WM(MXX) | |
1865 | DIMENSION R0(MXX), R1(MXX), R2(MXX), RH(MXX), RM(MXX) | |
1866 | S = EXP(TT) | |
1867 | Q = AL * EXP (S) | |
1868 | CPL = CtLhALPI(Q) | |
1869 | CPL2= CPL ** 2 / 2. * S | |
1870 | CPL = CPL * S | |
1871 | CALL CtLhINTEGR (NX,-1, FI, WM, IR1) | |
1872 | CALL CtLhINTEGR (NX, 0, FI, W0, IR2) | |
1873 | CALL CtLhINTEGR (NX, 1, FI, W1, IR3) | |
1874 | CALL CtLhINTEGR (NX,-1, GI, RM, IR4) | |
1875 | CALL CtLhINTEGR (NX, 0, GI, R0, IR5) | |
1876 | CALL CtLhINTEGR (NX, 1, GI, R1, IR6) | |
1877 | CALL CtLhINTEGR (NX, 2, GI, R2, IR7) | |
1878 | CALL CtLhHINTEG (NX, FI, WH) | |
1879 | CALL CtLhHINTEG (NX, GI, RH) | |
1880 | IF (IKNL .GT. 0) THEN | |
1881 | DO 230 IZ = 1, NX | |
1882 | FO(IZ) = ( 2D0 * FI(IZ) | |
1883 | > + 4D0 / 3D0 * ( 2D0 * WH(IZ) - W0(IZ) - W1(IZ) )) | |
1884 | > + NEFF * ( R0(IZ) - 2D0 * R1(IZ) + 2D0 * R2(IZ) ) | |
1885 | FO(IZ) = FO(IZ) * CPL | |
1886 | GO(IZ) = 4D0 / 3D0 * ( 2D0 * WM(IZ) - 2D0 * W0(IZ) + W1(IZ) ) | |
1887 | > + (33D0 - 2D0 * NEFF) / 6D0 * GI(IZ) | |
1888 | > + 6D0 * (RH(IZ) + RM(IZ) - 2D0 * R0(IZ) + R1(IZ) - R2(IZ)) | |
1889 | GO(IZ) = GO(IZ) * CPL | |
1890 | 230 CONTINUE | |
1891 | Else | |
1892 | DO 240 IZ = 1, NX | |
1893 | FO(IZ) = NEFF * (-R0(IZ) + 2.* R1(IZ) ) | |
1894 | > + 2.* FI(IZ) + 4./ 3.* ( 2.* WH(IZ) - W0(IZ) - W1(IZ) ) | |
1895 | FO(IZ) = FO(IZ) * CPL | |
1896 | GO(IZ) = 4./ 3.* ( 2.* W0(IZ) - W1(IZ) ) | |
1897 | >+ (33.- 2.* NEFF) / 6.* GI(IZ) + 6.*(RH(IZ) + R0(IZ) - 2.* R1(IZ)) | |
1898 | GO(IZ) = GO(IZ) * CPL | |
1899 | 240 CONTINUE | |
1900 | EndIf | |
1901 | IF (IKNL .EQ. 2) THEN | |
1902 | DZ = 1./(NX - 1) | |
1903 | DO 21 I = 1, NX-1 | |
1904 | NP = NX - I + 1 | |
1905 | IS = NP | |
1906 | g2(1)=0d0 | |
1907 | g3(1)=0d0 | |
1908 | DO 31 KZ = 2, NP | |
1909 | IY = I + KZ - 1 | |
1910 | IT = NX - IY + 1 | |
1911 | XY = ZZ (IS, IT) | |
1912 | G1(KZ) = FFG(I, IY) * (FI(IY) - XY**2 *FI(I)) | |
1913 | G4(KZ) = GGF(I, IY) * (GI(IY) - XY**2 *GI(I)) | |
1914 | G2(KZ) = FFG(IS,IT) * (GI(IY) - xy*GI(I)) !FG | |
1915 | G3(KZ) = GGF(IS,IT) * (FI(IY) - XY*FI(I)) !GF (usual notations) | |
1916 | 31 CONTINUE | |
1917 | TEM1 = CtLhSMPNOL (NP, DZ, G1, ERR) | |
1918 | TEM2 = CtLhSMPSNA (NP, DZ, G2, ERR) | |
1919 | TEM3 = CtLhSMPSNA (NP, DZ, G3, ERR) | |
1920 | TEM4 = CtLhSMPNOL (NP, DZ, G4, ERR) | |
1921 | TEM1 = TEM1 - FI(I) * (AFF2(I) + ZGF2) | |
1922 | TEM4 = TEM4 - GI(I) * (AGG2(I) + ZFG2) | |
1923 | tem2 = tem2 + GI(I)*AFG2(I) | |
1924 | tem3= tem3 + FI(I)*AGF2(I) | |
1925 | TMF = TEM1 + TEM2 | |
1926 | TMG = TEM3 + TEM4 | |
1927 | FO(I) = FO(I) + TMF * CPL2 | |
1928 | GO(I) = GO(I) + TMG * CPL2 | |
1929 | 21 CONTINUE | |
1930 | EndIf | |
1931 | RETURN | |
1932 | END | |
1933 | FUNCTION CtLhSPENC2 (X) | |
1934 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1935 | EXTERNAL CtLhSPN2IN | |
1936 | COMMON / LhCtSPENCC / XX | |
1937 | DATA U1, AERR, RERR / 1.D0, 1.E-7, 5.E-3 / | |
1938 | XX = X | |
1939 | TEM = CtLhGausInt(CtLhSPN2IN, XX, U1, AERR, RERR, ERR, IRT) | |
1940 | CtLhSPENC2 = TEM + LOG (XX) ** 2 / 2. | |
1941 | RETURN | |
1942 | END | |
1943 | FUNCTION CtLhSPN2IN (ZZ) | |
1944 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1945 | COMMON / LhCtSPENCC / X | |
1946 | Z = ZZ | |
1947 | TEM = LOG (1.+ X - Z) / Z | |
1948 | CtLhSPN2IN = TEM | |
1949 | RETURN | |
1950 | END | |
1951 | SUBROUTINE CtLhSTUPKL (NFL) | |
1952 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
1953 | LOGICAL LSTX | |
1954 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
1955 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
1956 | PARAMETER (MX = 3) | |
1957 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
1958 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
1959 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
1960 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
1961 | COMMON / LhCtKRN1ST / FF1(0:MXX),FG1(0:MXX),GF1(0:MXX),GG1(0:MXX), | |
1962 | > FF2(0:MXX), FG2(0:MXX), GF2(0:MXX), GG2(0:MXX), | |
1963 | > PNSP(0:MXX), PNSM(0:MXX) | |
1964 | COMMON / LhCtKRN2ND / FFG(MXX, MXX), GGF(MXX, MXX), PNS(MXX, MXX) | |
1965 | COMMON / LhCtKRNL00 / DZ, XL(MX), NNX | |
1966 | COMMON / LhCtKRNL01 / AFF2(MXX),AFG2(MXX),AGF2(MXX),AGG2(MXX), | |
1967 | > ANSP (MXX), ANSM (MXX), ZFG2, ZGF2, ZQQB | |
1968 | EXTERNAL CtLhPFF1, CtLhRGG1, CtLhRFF2, CtLhRGG2 | |
1969 | EXTERNAL CtLhFNSP, CtLhFNSM | |
1970 | dimension aff1(mxx),agg1(mxx) | |
1971 | PARAMETER (PI = 3.141592653589793d0, PI2 = PI**2) | |
1972 | DATA CF, CG, TR / 1.33333333333333d0, 3.0, 0.5 / | |
1973 | data zeta3/1.20205690315959d0/ ! zeta(3.0) | |
1974 | SAVE | |
1975 | DATA AERR, RERR / 0.0, 0.02 / | |
1976 | NNX = NX | |
1977 | DZ = 1./ (NX - 1) | |
1978 | DO 5 I0 = 1, MX | |
1979 | XL(I0) = XV(I0) | |
1980 | 5 CONTINUE | |
1981 | DO 10 I = 1, NX-1 | |
1982 | XZ = XV(I) | |
1983 | CALL CtLhKERNEL (XZ, FF1(I), GF1(I), FG1(I), GG1(I), PNSP(I), | |
1984 | > PNSM(I), FF2(I), GF2(I), FG2(I), GG2(I), NFL, IRT) | |
1985 | 10 CONTINUE | |
1986 | FF1(0) = FF1(1) * 3. - FF1(2) * 3. + FF1(3) | |
1987 | FG1(0) = FG1(1) * 3. - FG1(2) * 3. + FG1(3) | |
1988 | GF1(0) = GF1(1) * 3. - GF1(2) * 3. + GF1(3) | |
1989 | GG1(0) = GG1(1) * 3. - GG1(2) * 3. + GG1(3) | |
1990 | PNSP(0) = PNSP(1) * 3. - PNSP(2) * 3. + PNSP(3) | |
1991 | PNSM(0) = PNSM(1) * 3. - PNSM(2) * 3. + PNSM(3) | |
1992 | FF2(0) = FF2(1) * 3. - FF2(2) * 3. + FF2(3) | |
1993 | FG2(0) = FG2(1) * 3. - FG2(2) * 3. + FG2(3) | |
1994 | GF2(0) = GF2(1) * 3. - GF2(2) * 3. + GF2(3) | |
1995 | GG2(0) = GG2(1) * 3. - GG2(2) * 3. + GG2(3) | |
1996 | FF1(NX) = FF1(NX-1) * 3. - FF1(NX-2) * 3. + FF1(NX-3) | |
1997 | FG1(NX) = FG1(NX-1) * 3. - FG1(NX-2) * 3. + FG1(NX-3) | |
1998 | GF1(NX) = GF1(NX-1) * 3. - GF1(NX-2) * 3. + GF1(NX-3) | |
1999 | GG1(NX) = GG1(NX-1) * 3. - GG1(NX-2) * 3. + GG1(NX-3) | |
2000 | PNSM(NX) = PNSM(NX-1) * 3. - PNSM(NX-2) * 3. + PNSM(NX-3) | |
2001 | PNSP(NX) = PNSP(NX-1) * 3. - PNSP(NX-2) * 3. + PNSP(NX-3) | |
2002 | FF2(NX) = FF2(NX-1) * 3. - FF2(NX-2) * 3. + FF2(NX-3) | |
2003 | FG2(NX) = FG2(NX-1) * 3. - FG2(NX-2) * 3. + FG2(NX-3) | |
2004 | GF2(NX) = GF2(NX-1) * 3. - GF2(NX-2) * 3. + GF2(NX-3) | |
2005 | GG2(NX) = GG2(NX-1) * 3. - GG2(NX-2) * 3. + GG2(NX-3) | |
2006 | RER = RERR * 4. | |
2007 | AFF1(1) = CtLhGausInt(CtLhPFF1,D0,XV(1),AERR,RERR,ER1,IRT) | |
2008 | DGG1 = NFL / 3. | |
2009 | TMPG = CtLhGausInt(CtLhRGG1,D0,XV(1),AERR,RERR,ER3,IRT) | |
2010 | AGG1(1) = TMPG + DGG1 | |
2011 | ANSM(1) = CtLhGausInt(CtLhFNSM,D0,XV(1),AERR,RER,ER2,IRT) | |
2012 | ANSP(1) = CtLhGausInt(CtLhFNSP,D0,XV(1),AERR,RER,ER2,IRT) | |
2013 | AER = AFF1(1) * RER | |
2014 | AFF2(1) = CtLhGausInt(CtLhRFF2, D0, XV(1), AER, RER, ER2, IRT) | |
2015 | AER = AGG1(1) * RER | |
2016 | AGG2(1) = CtLhGausInt(CtLhRGG2, D0, XV(1), AER, RER, ER4, IRT) | |
2017 | DO 20 I2 = 2, NX-1 | |
2018 | TEM =CtLhGausInt(CtLhPFF1,XV(I2-1),XV(I2),AERR,RERR,ER1,IRT) | |
2019 | AFF1(I2) = TEM + AFF1(I2-1) | |
2020 | AER = ABS(TEM * RER) | |
2021 | AFF2(I2)=CtLhGausInt(CtLhRFF2,XV(I2-1),XV(I2),AER,RER,ER2,IRT) | |
2022 | > +AFF2(I2-1) | |
2023 | TEM = CtLhGausInt(CtLhRGG1,XV(I2-1),XV(I2),AERR,RERR,ER3,IRT) | |
2024 | TMPG = TMPG + TEM | |
2025 | AGG1(I2) = TMPG + DGG1 | |
2026 | AER = ABS(TEM * RER) | |
2027 | AGG2(I2)=CtLhGausInt(CtLhRGG2,XV(I2-1),XV(I2),AER,RER,ER4,IRT) | |
2028 | > +AGG2(I2-1) | |
2029 | ANSP(I2)=CtLhGausInt(CtLhFNSP,XV(I2-1),XV(I2),AERR,RER,ER4,IRT) | |
2030 | > +ANSP(I2-1) | |
2031 | ANSM(I2)=CtLhGausInt(CtLhFNSM,XV(I2-1),XV(I2),AERR,RER,ER4,IRT) | |
2032 | > +ANSM(I2-1) | |
2033 | 20 CONTINUE | |
2034 | ANSP(NX)=CtLhGausInt(CtLhFNSP,XV(NX-1),D1,AERR,RER,ERR, | |
2035 | > IRT) + ANSP(NX-1) | |
2036 | ANSM(NX)=CtLhGausInt(CtLhFNSM,XV(NX-1),D1,AERR,RER,ERR, | |
2037 | > IRT) + ANSM(NX-1) | |
2038 | TRNF = TR * NFL | |
2039 | do i2=1,nx-1 !loop over x | |
2040 | x=xv(i2) | |
2041 | c XI = 1./ X !unused - jcp | |
2042 | X2 = X ** 2 | |
2043 | X3= x**3 | |
2044 | XLN = DLOG (X) | |
2045 | XLN2 = XLN ** 2 | |
2046 | XLN1M = DLOG (1.- X) | |
2047 | xLi2m=CtLhxLi(2,-x) | |
2048 | xLi2=CtLhxLi(2,x) | |
2049 | xLi3=CtLhxLi(3,x) | |
2050 | xLi31m=CtLhxLi(3,1d0-x) | |
2051 | xLi32=CtLhxLi(3,x2) | |
2052 | xln1m2=xln1m*xln1m | |
2053 | xln1p=dlog(1d0+x) | |
2054 | x1m=1d0-x | |
2055 | x1p=1d0+x | |
2056 | x3m=3d0-x | |
2057 | x3p=3d0+x | |
2058 | wgfcft= | |
2059 | > (9 + 4*Pi2 - 22*x + 13*x2 + 6*(3 - 4*x + x2)*xln1m + | |
2060 | > 40*xln - 24*xLi2)/9. | |
2061 | wgfcf2= | |
2062 | > (6*(2*(-9 + Pi2) + 3*x*(5 + x)) +4*(3 +2*Pi2+3*x*(-3 + 2*x))* | |
2063 | > xln1m + 6*x3m*x1m*xln1m2 - 6*(x*(8 + 3*x) + 4*xln1m2)* | |
2064 | > xln - 3*(-4 + x)*x*xln2)/12 - 2*(3 + 2*xln1m)*xLi2 - 4*xLi31m | |
2065 | wgfcfg= | |
2066 | > (3637-186*Pi2-x*(3198+72*Pi2+x*(231 + 208*x)))/108.- xln + | |
2067 | > (3*xln1m*(-33 - 4*Pi2 + (50 - 17*x)*x - 3*x3m*x1m*xln1m) + | |
2068 | > 2*(x*(198 + x*(27+8*x))+9*xln1m*(3 - 4*x + x2 + 2*xln1m))* | |
2069 | > xln - 9*x*(4 + x)*xln2)/18- x1p*x3p*xln*xln1p- | |
2070 | > (x1p*x3p - 4*xln)*xLi2m + (31d0/3d0 +4*xln1m- 4*xln)*xLi2 + | |
2071 | > 4*xLi31m + 12*xLi3 - 2*xLi32 - 10*zeta3 | |
2072 | wfgcft= | |
2073 | > (18 - 81*x + 6*Pi2*x + 123*x2 - 6*Pi2*x2 - 60*x3 + | |
2074 | > 4*Pi2*x3 - 6*(-2 + 3*x - 3*x2 + 2*x3)*xln1m2 -33*x*xln + | |
2075 | > 15*x2*xln - 24*x3*xln - 9*x*xln2 + 9*x2*xln2 - | |
2076 | > 12*x3*xln2 - 12*x1m*xln1m*(-1 + 2*x2 + 2*xln - x*xln + | |
2077 | > 2*x2*xln) - 24*xLi2)/9. | |
2078 | wfgcgt= | |
2079 | > (2*(-67 + 2*Pi2 + x*(64 + x*(-91 + 3*Pi2 + 94*x)) + | |
2080 | > x1m*(7+x*(-5+16*x))*xln1m -3*x1m*(2+ x*(-1+2*x))*xln1m2 - | |
2081 | > 20*xln - 3*x*xln*(13 + 16*x*x1p - 3*x1p*xln) + | |
2082 | > 6*x1p*(2+x+2*x2)*xln*xln1p+6*x1p*(2+x+2*x2)*xLi2m))/9. | |
2083 | AGF2(I2) = CF*TRNF*WGFCFT+CF**2* WGFCF2+CF*CG*WGFCFG | |
2084 | AFG2(I2) = CF*TRNF*WFGCFT +CG*TRNF*WFGCGT | |
2085 | enddo !i2 | |
2086 | AGF2(nx)=0d0 | |
2087 | AFG2(nx)=0d0 | |
2088 | ZGF2=-28./27.*Cf**2+94./27.*Cf*Cg -52./27.*Cf*TrNf | |
2089 | ZFG2= 37./27.*Cf*TrNf + 35./54.*Cg*TrNf | |
2090 | ZQQB=1.43862321154902*(Cf**2-0.5*Cf*Cg) | |
2091 | DO 21 IX = 1, NX-1 | |
2092 | X = XV(IX) | |
2093 | NP = NX - IX + 1 | |
2094 | IS = NP | |
2095 | XG2 = (LOG(1./(1.-X)) + 1.) ** 2 | |
2096 | FFG (IS, IS) = FG2(NX) * DXTZ(I) * XG2 | |
2097 | GGF (IS, IS) = GF2(NX) * DXTZ(I) * XG2 | |
2098 | PNS (IS, IS) =PNSM(NX) * DXTZ(I) | |
2099 | DO 31 KZ = 2, NP | |
2100 | IY = IX + KZ - 1 | |
2101 | IT = NX - IY + 1 | |
2102 | XY = X / XV(IY) | |
2103 | XM1 = 1.- XY | |
2104 | XG2 = (LOG(1./XM1) + 1.) ** 2 | |
2105 | Z = ZZ (IX, IY) | |
2106 | TZ = (Z + DZ) / DZ | |
2107 | IZ = Int(TZ) | |
2108 | IZ = MAX (IZ, 0) | |
2109 | IZ = MIN (IZ, NX-1) | |
2110 | DT = TZ - IZ | |
2111 | TEM = (FF2(IZ) * (1.- DT) + FF2(IZ+1) * DT) / XM1 / XY | |
2112 | FFG (IX, IY) = TEM * DXTZ(IY) | |
2113 | TEM = (FG2(IZ) * (1.- DT) + FG2(IZ+1) * DT) * XG2 / XY | |
2114 | FFG (IS, IT) = TEM * DXTZ(IY) | |
2115 | TEM = (GF2(IZ) * (1.- DT) + GF2(IZ+1) * DT) * XG2 / XY | |
2116 | GGF (IS, IT) = TEM * DXTZ(IY) | |
2117 | TEM = (GG2(IZ) * (1.- DT) + GG2(IZ+1) * DT) / XM1 / XY | |
2118 | GGF (IX, IY) = TEM * DXTZ(IY) | |
2119 | TEM = (PNSP(IZ) * (1.- DT) + PNSP(IZ+1) * DT) / XM1 | |
2120 | PNS (IX, IY) = TEM * DXTZ(IY) | |
2121 | TEM = (PNSM(IZ) * (1.- DT) + PNSM(IZ+1) * DT) / XM1 | |
2122 | PNS (IS, IT) = TEM * DXTZ(IY) | |
2123 | 31 CONTINUE | |
2124 | 21 CONTINUE | |
2125 | RETURN | |
2126 | END | |
2127 | SUBROUTINE CtLhTRNLAM (IRDR, NF, IACT, IRT) | |
2128 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2129 | COMMON / LhCtCWZPRM / ALAM(0:9), AMHAT(0:9), AMN, NHQ | |
2130 | COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1 | |
2131 | EXTERNAL CtLhZBRLAM | |
2132 | DATA ALM0, BLM0, RERR / 0.01, 10.0, 0.0001 / | |
2133 | DATA IR1, SML / 0, 1.E-5 / | |
2134 | IRT = 0 | |
2135 | N = NF | |
2136 | JRDR = IRDR | |
2137 | JACT = IACT | |
2138 | VLAM = ALAM(N) | |
2139 | IF (JACT .GT. 0) THEN | |
2140 | N1 = N + 1 | |
2141 | THMS = AMHAT(N1) | |
2142 | ALM = LOG (THMS/VLAM) | |
2143 | BLM = BLM0 | |
2144 | ELSE | |
2145 | N1 = N -1 | |
2146 | THMS = AMHAT(N) | |
2147 | ALM = ALM0 | |
2148 | THMS = MAX (THMS, SML) | |
2149 | BLM = LOG (THMS/VLAM) | |
2150 | ENDIF | |
2151 | IF (VLAM .GE. 0.7 * THMS) THEN | |
2152 | IF (JACT . EQ. 1) THEN | |
2153 | AMHAT(N1) = 0 | |
2154 | ELSE | |
2155 | AMHAT(N) = 0 | |
2156 | ENDIF | |
2157 | IRT = 4 | |
2158 | ALAM(N1) = VLAM | |
2159 | RETURN | |
2160 | ENDIF | |
2161 | IF (ALM .GE. BLM) THEN | |
2162 | WRITE (*, *) 'CtLhTRNLAM has ALM >= BLM: ', ALM, BLM | |
2163 | WRITE (*, *) 'I do not know how to continue' | |
2164 | STOP | |
2165 | ENDIF | |
2166 | VMULM = THMS/VLAM | |
2167 | ERR = RERR * LOG (VMULM) | |
2168 | WLLN = CtLhQZBRNT (CtLhZBRLAM, ALM, BLM, ERR, IR1) | |
2169 | ALAM(N1) = THMS / EXP (WLLN) | |
2170 | IF (IR1 .NE. 0) THEN | |
2171 | WRITE (*, *) 'CtLhQZBRNT failed in CtLhTRNLAM; ', | |
2172 | > 'NF, VLAM =', NF, VLAM | |
2173 | WRITE (*, *) 'I do not know how to continue' | |
2174 | STOP | |
2175 | ENDIF | |
2176 | RETURN | |
2177 | END | |
2178 | SUBROUTINE CtLhUPC (A, La, UpA) | |
2179 | CHARACTER A*(*), UpA*(*), C*(1) | |
2180 | INTEGER I, La, Ld | |
2181 | La = Len(A) | |
2182 | Lb = Len(UpA) | |
2183 | If (Lb .Lt. La) Stop 'UpCase conversion length mismatch!' | |
2184 | Ld = ICHAR('A')-ICHAR('a') | |
2185 | DO 1 I = 1, Lb | |
2186 | If (I .Le. La) Then | |
2187 | c = A(I:I) | |
2188 | IF ( LGE(C, 'a') .AND. LLE(C, 'z') ) THEN | |
2189 | UpA (I:I) = CHAR(Ichar(c) + ld) | |
2190 | Else | |
2191 | UpA (I:I) = C | |
2192 | ENDIF | |
2193 | Else | |
2194 | UpA (I:I) = ' ' | |
2195 | Endif | |
2196 | 1 CONTINUE | |
2197 | ||
2198 | RETURN | |
2199 | END | |
2200 | SUBROUTINE CtLhWARNI (IWRN, MSG, NMVAR, IVAB, | |
2201 | > IMIN, IMAX, IACT) | |
2202 | CHARACTER*(*) MSG, NMVAR | |
2203 | Save Iw | |
2204 | Data Nmax / 100 / | |
2205 | IW = IWRN | |
2206 | IV = IVAB | |
2207 | ||
2208 | IF (IW .EQ. 0) THEN | |
2209 | PRINT '(1X,A/1X, 2A,I10)', MSG, NMVAR,' = ', IV | |
2210 | IF (IACT .EQ. 1) THEN | |
2211 | PRINT '(A/2I10)', ' The limits are: ', IMIN, IMAX | |
2212 | ENDIF | |
2213 | ENDIF | |
2214 | If (Iw .LT. Nmax) Then | |
2215 | PRINT '(1X,A/1X, 2A,I10)', MSG, NMVAR,' = ', IV | |
2216 | Elseif (Iw .Eq. Nmax) Then | |
2217 | Print '(/A/)', 'CtLhWARNI Severe Warning: Too many errors' | |
2218 | Endif | |
2219 | IWRN = IW + 1 | |
2220 | RETURN | |
2221 | END | |
2222 | SUBROUTINE CtLhWARNR (IWRN, MSG, NMVAR, VARIAB, | |
2223 | > VMIN, VMAX, IACT) | |
2224 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2225 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
2226 | CHARACTER*(*) MSG, NMVAR | |
2227 | Save Iw | |
2228 | Data Nmax / 100 / | |
2229 | IW = IWRN | |
2230 | VR = VARIAB | |
2231 | IF (IW .EQ. 0) THEN | |
2232 | PRINT '(1X, A/1X,2A,1PD16.7)', MSG, NMVAR, ' = ', VR | |
2233 | IF (IACT .EQ. 1) THEN | |
2234 | PRINT '(A/2(1PE15.4))', ' The limits are: ', VMIN, VMAX | |
2235 | ENDIF | |
2236 | ENDIF | |
2237 | If (Iw .LT. Nmax) Then | |
2238 | PRINT '(I5, 2A/1X,2A,1PD16.7)', IW, ' ', MSG, | |
2239 | > NMVAR, ' = ', VR | |
2240 | Elseif (Iw .Eq. Nmax) Then | |
2241 | Print '(/A/)', 'CtLhWARNR Severe Warning: Too many errors' | |
2242 | Endif | |
2243 | IWRN = IW + 1 | |
2244 | RETURN | |
2245 | END | |
2246 | SUBROUTINE CtLhXARRAY | |
2247 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2248 | LOGICAL LSTX | |
2249 | PARAMETER (D0 = 0.0, D10=10.0) | |
2250 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
2251 | PARAMETER (MXPN = MXF * 2 + 2) | |
2252 | PARAMETER (MXQX= MXQ * MXX, MXPQX = MXQX * MXPN) | |
2253 | PARAMETER (M1=-3, M2=3, NDG=3, NDH=NDG+1, L1=M1-1, L2=M2+NDG-2) | |
2254 | Character Msg*80 | |
2255 | COMMON / LhCtVARIBX / XA(MXX, L1:L2), ELY(MXX), DXTZ(MXX) | |
2256 | COMMON / LhCtVARBAB / GB(NDG, NDH, MXX), H(NDH, MXX, M1:M2) | |
2257 | COMMON / LhCtHINTEC / GH(NDG, MXX) | |
2258 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
2259 | COMMON / LhCtXYARAY / ZZ(MXX, MXX), ZV(0:MXX) | |
2260 | DIMENSION G1(NDG,NDH), G2(NDG,NDH), A(NDG) | |
2261 | DATA F12, F22, F32 / 1D0, 1D0, 1D0 / | |
2262 | DATA (G1(I,NDH), G2(I,1), I=1,NDG) / 0.0,0.0,0.0,0.0,0.0,0.0 / | |
2263 | DATA PUNY / 1D-30 / | |
2264 | XV(0) = 0D0 | |
2265 | DZ = 1D0 / (NX-1) | |
2266 | DO 10 I = 1, NX - 1 | |
2267 | Z = DZ * (I-1) | |
2268 | ZV(I) = Z | |
2269 | X = CtLhXFRMZ (Z) | |
2270 | DXTZ(I) = CtLhDXDZ(Z) / X | |
2271 | XV (I) = X | |
2272 | XA(I, 1) = X | |
2273 | XA(I, 0) = LOG (X) | |
2274 | DO 20 L = L1, L2 | |
2275 | IF (L .NE. 0 .AND. L .NE. 1) XA(I, L) = X ** L | |
2276 | 20 CONTINUE | |
2277 | 10 CONTINUE | |
2278 | XV(1) = Xmin | |
2279 | XV(NX) = 1D0 | |
2280 | ZV(Nx) = 1D0 | |
2281 | DXTZ(NX) = CtLhDXDZ(1.D0) | |
2282 | DO 21 L = L1, L2 | |
2283 | XA (NX, L) = 1D0 | |
2284 | 21 CONTINUE | |
2285 | XA (NX, 0) = 0D0 | |
2286 | DO 11 I = 1, NX-1 | |
2287 | ELY(I) = LOG(1D0 - XV(I)) | |
2288 | 11 CONTINUE | |
2289 | ELY(NX) = 3D0* ELY(NX-1) - 3D0* ELY(NX-2) + ELY(NX-3) | |
2290 | DO 17 IX = 1, NX | |
2291 | ZZ (IX, IX) = 1. | |
2292 | DO 17 IY = IX+1, NX | |
2293 | XY = XV(IX) / XV(IY) | |
2294 | ZZ (IX, IY) = CtLhZFRMX (XY) | |
2295 | ZZ (NX-IX+1, NX-IY+1) = XY | |
2296 | 17 CONTINUE | |
2297 | DO 30 I = 1, NX-1 | |
2298 | IF (I .NE. NX-1) THEN | |
2299 | F11 = 1D0/XV(I) | |
2300 | F21 = 1D0/XV(I+1) | |
2301 | F31 = 1D0/XV(I+2) | |
2302 | F13 = XV(I) | |
2303 | F23 = XV(I+1) | |
2304 | F33 = XV(I+2) | |
2305 | DET = F11*F22*F33 + F21*F32*F13 + F31*F12*F23 | |
2306 | > - F31*F22*F13 - F21*F12*F33 - F11*F32*F23 | |
2307 | IF (ABS(DET) .LT. PUNY) THEN | |
2308 | Msg='Determinant close to zero; will be arbitrarily set to:' | |
2309 | CALL CtLhWARNR(IWRN, Msg, 'DET', PUNY, D0, D0, 0) | |
2310 | DET = PUNY | |
2311 | EndIf | |
2312 | G2(1,2) = (F22*F33 - F23*F32) / DET | |
2313 | G2(1,3) = (F32*F13 - F33*F12) / DET | |
2314 | G2(1,4) = (F12*F23 - F13*F22) / DET | |
2315 | G2(2,2) = (F23*F31 - F21*F33) / DET | |
2316 | G2(2,3) = (F33*F11 - F31*F13) / DET | |
2317 | G2(2,4) = (F13*F21 - F11*F23) / DET | |
2318 | G2(3,2) = (F21*F32 - F22*F31) / DET | |
2319 | G2(3,3) = (F31*F12 - F32*F11) / DET | |
2320 | G2(3,4) = (F11*F22 - F12*F21) / DET | |
2321 | B2 = LOG (XV(I+2)/XV(I)) | |
2322 | B3 = XV(I) * (B2 - 1.) + XV(I+2) | |
2323 | GH (1,I) = B2 * G2 (2,2) + B3 * G2 (3,2) | |
2324 | GH (2,I) = B2 * G2 (2,3) + B3 * G2 (3,3) | |
2325 | GH (3,I) = B2 * G2 (2,4) + B3 * G2 (3,4) | |
2326 | EndIf | |
2327 | DO 51 J = 1, NDH | |
2328 | DO 52 L = 1, NDG | |
2329 | IF (I .EQ. 1) THEN | |
2330 | GB(L,J,I) = G2(L,J) | |
2331 | ElseIF (I .EQ. NX-1) THEN | |
2332 | GB(L,J,I) = G1(L,J) | |
2333 | Else | |
2334 | GB(L,J,I) = (G1(L,J) + G2(L,J)) / 2D0 | |
2335 | EndIf | |
2336 | 52 CONTINUE | |
2337 | 51 CONTINUE | |
2338 | DO 35 MM = M1, M2 | |
2339 | DO 40 K = 1, NDG | |
2340 | KK = K + MM - 2 | |
2341 | IF (KK .EQ. 0) THEN | |
2342 | A(K) = XA(I+1, 0) - XA(I, 0) | |
2343 | Else | |
2344 | A(K) = (XA(I+1, KK) - XA(I, KK)) / DBLE(KK) | |
2345 | EndIf | |
2346 | 40 CONTINUE | |
2347 | DO 41 J = 1, NDH | |
2348 | TEM = 0 | |
2349 | DO 43 L = 1, NDG | |
2350 | TEM = TEM + A(L) * GB(L,J,I) | |
2351 | 43 CONTINUE | |
2352 | H(J,I,MM) = TEM | |
2353 | 41 CONTINUE | |
2354 | 35 CONTINUE | |
2355 | DO 42 J = 1, NDG | |
2356 | DO 44 L = 1, NDG | |
2357 | G1(L,J) = G2(L,J+1) | |
2358 | 44 CONTINUE | |
2359 | 42 CONTINUE | |
2360 | 30 CONTINUE | |
2361 | LSTX = .TRUE. | |
2362 | RETURN | |
2363 | END | |
2364 | FUNCTION CtLhXFRMZ (Z) | |
2365 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2366 | LOGICAL LSTX | |
2367 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
2368 | PARAMETER (MXX = 105) | |
2369 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
2370 | COMMON / LhCtINVERT / ZA | |
2371 | EXTERNAL CtLhZFXL | |
2372 | DATA TEM, RER / D1, 1E-3 / | |
2373 | DATA ZLOW, ZHIGH, IWRN2 / -10.0, 1.00002, 0 / | |
2374 | EPS = TEM * RER | |
2375 | ZA = Z | |
2376 | IF (Z .LE. ZHIGH .AND. Z .GT. ZLOW) THEN | |
2377 | XLA = LOG (XMIN) * 1.5 | |
2378 | XLB = 0.00001 | |
2379 | TEM = CtLhZBRNT (CtLhZFXL, XLA, XLB, EPS, IRT) | |
2380 | Else | |
2381 | CALL CtLhWARNR (IWRN2, 'Z out of range in CtLhXFRMZ, X set=0.', | |
2382 | > 'Z', Z, ZLOW, ZHIGH, 1) | |
2383 | TEM = 0 | |
2384 | EndIf | |
2385 | CtLhXFRMZ = EXP(TEM) | |
2386 | RETURN | |
2387 | END | |
2388 | FUNCTION CtLhxLi(n,x) | |
2389 | implicit NONE | |
2390 | integer NCUT, i,n,m3 | |
2391 | real*8 CtLhxLi,Out,x,pi2by6,zeta3,c1,c2 | |
2392 | real*8 r,xt,L,xln1m | |
2393 | parameter (m3=8) | |
2394 | dimension c1(2:m3),c2(2:m3) | |
2395 | data NCUT/27/ | |
2396 | data c1/0.75,-0.5833333333333333d0,0.454861111111111d0, | |
2397 | > -0.3680555555555555d0,0.3073611111111111d0, | |
2398 | > -0.2630555555555555d0,0.2294880243764172d0/ | |
2399 | data c2/-0.5d0,0.5d0,-0.4583333333333333d0,0.416666666666666d0, | |
2400 | > -0.3805555555555555d0,0.35d0,-0.3241071428571428d0/ | |
2401 | data zeta3,pi2by6 /1.20205690315959d0,1.64493406684823d0/ | |
2402 | L=0.0 | |
2403 | i=0 | |
2404 | r=1.0 | |
2405 | if (abs(x).gt.r) then | |
2406 | PRINT *,'Li: x out of range (-1,1) , x=',x | |
2407 | STOP | |
2408 | endif | |
2409 | if (n.lt.0) then | |
2410 | PRINT *,'Polylogarithm Li undefined for n=',n | |
2411 | STOP | |
2412 | elseif (n.eq.0) then | |
2413 | Out=x/(1d0-x) | |
2414 | elseif (n.eq.1) then | |
2415 | Out=-dlog(1-x) | |
2416 | elseif (n.eq.2) then | |
2417 | !Calculate dilogarithm | |
2418 | !separately for x<0.5 and x>0.5 | |
2419 | if (x.ge.(-0.5).and.x.le.0.5) then | |
2420 | do while(i.le.NCUT) | |
2421 | i=i+1 | |
2422 | r=r*x | |
2423 | L=L+r/i/i | |
2424 | enddo | |
2425 | Out=L | |
2426 | elseif (x.eq.0) then | |
2427 | Out=0d0 | |
2428 | elseif(x.gt.0.5) then !n.eq.2,x>0.5 | |
2429 | xt = 1.0-x | |
2430 | L = pi2by6 - dlog(x)*dlog(xt) | |
2431 | do while(i.le.NCUT) | |
2432 | i=i+1 | |
2433 | r=r*xt | |
2434 | L=L-r/i/i | |
2435 | enddo | |
2436 | Out=L | |
2437 | elseif (x.lt.(-0.5)) then | |
2438 | xt=-x/(1d0-x) | |
2439 | L=-0.5*dlog(1-x)**2 | |
2440 | do while (i.le.NCUT) | |
2441 | i=i+1 | |
2442 | r=r*xt | |
2443 | L=L-r/i/i | |
2444 | enddo | |
2445 | Out=L | |
2446 | endif | |
2447 | elseif (n.eq.3.and.x.ge.0.8) then !use the expansion of Li3 near x=1 | |
2448 | L=zeta3+pi2by6*dlog(x) | |
2449 | xt=(1d0-x) | |
2450 | xln1m=dlog(xt) | |
2451 | do i=2,m3 | |
2452 | L=L+(c1(i)+c2(i)*xln1m)*xt**i | |
2453 | enddo | |
2454 | Out=L | |
2455 | else !n>3 or x=3,x<0.8 | |
2456 | do while(i.le.NCUT) | |
2457 | i=i+1 | |
2458 | r=r*x | |
2459 | L=L+r/dble(i)**dble(n) | |
2460 | enddo | |
2461 | Out=L | |
2462 | endif | |
2463 | CtLhxLi=Out | |
2464 | End ! CtLhxLi | |
2465 | FUNCTION CtLhZBRLAM (WLLN) | |
2466 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2467 | COMMON / LhCtTRNCOM / VMULM, JRDR, N, N1 | |
2468 | WMULM = EXP (WLLN) | |
2469 | TEM1 = 1./ CtLhALPQCD(JRDR, N1, WMULM, I) | |
2470 | TEM2 = 1./ CtLhALPQCD(JRDR, N, VMULM, I) | |
2471 | CtLhZBRLAM = TEM1 - TEM2 | |
2472 | END | |
2473 | FUNCTION CtLhZBRNT(FUNC, X1, X2, TOL, IRT) | |
2474 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2475 | PARAMETER (ITMAX = 1000, EPS = 3.E-12) | |
2476 | external func | |
2477 | IRT = 0 | |
2478 | TOL = ABS(TOL) | |
2479 | A=X1 | |
2480 | B=X2 | |
2481 | FA=FUNC(A) | |
2482 | FB=FUNC(B) | |
2483 | IF(FB*FA.GT.0.) THEN | |
2484 | PRINT *, 'Root must be bracketed for CtLhZBRNT. Set = 0' | |
2485 | IRT = 1 | |
2486 | CtLhZBRNT=0. | |
2487 | RETURN | |
2488 | ENDIF | |
2489 | FC=FB | |
2490 | DO 11 ITER=1,ITMAX | |
2491 | IF(FB*FC.GT.0.) THEN | |
2492 | C=A | |
2493 | FC=FA | |
2494 | D=B-A | |
2495 | E=D | |
2496 | ENDIF | |
2497 | IF(ABS(FC).LT.ABS(FB)) THEN | |
2498 | A=B | |
2499 | B=C | |
2500 | C=A | |
2501 | FA=FB | |
2502 | FB=FC | |
2503 | FC=FA | |
2504 | ENDIF | |
2505 | TOL1=2.*EPS*ABS(B)+0.5*TOL | |
2506 | XM=.5*(C-B) | |
2507 | IF(ABS(XM).LE.TOL1 .OR. FB.EQ.0.)THEN | |
2508 | CtLhZBRNT=B | |
2509 | RETURN | |
2510 | ENDIF | |
2511 | IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN | |
2512 | S=FB/FA | |
2513 | IF(A.EQ.C) THEN | |
2514 | P=2.*XM*S | |
2515 | Q=1.-S | |
2516 | ELSE | |
2517 | Q=FA/FC | |
2518 | R=FB/FC | |
2519 | P=S*(2.*XM*Q*(Q-R)-(B-A)*(R-1.)) | |
2520 | Q=(Q-1.)*(R-1.)*(S-1.) | |
2521 | ENDIF | |
2522 | IF(P.GT.0.) Q=-Q | |
2523 | P=ABS(P) | |
2524 | IF(2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN | |
2525 | E=D | |
2526 | D=P/Q | |
2527 | ELSE | |
2528 | D=XM | |
2529 | E=D | |
2530 | ENDIF | |
2531 | ELSE | |
2532 | D=XM | |
2533 | E=D | |
2534 | ENDIF | |
2535 | A=B | |
2536 | FA=FB | |
2537 | IF(ABS(D) .GT. TOL1) THEN | |
2538 | B=B+D | |
2539 | ELSE | |
2540 | B=B+SIGN(TOL1,XM) | |
2541 | ENDIF | |
2542 | FB=FUNC(B) | |
2543 | 11 CONTINUE | |
2544 | PRINT *, 'CtLhZBRNT exceeding maximum iterations.' | |
2545 | IRT = 2 | |
2546 | CtLhZBRNT=B | |
2547 | RETURN | |
2548 | END | |
2549 | FUNCTION CtLhZFRMX (XX) | |
2550 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2551 | LOGICAL LSTX | |
2552 | PARAMETER (D0=0D0, D1=1D0, D2=2D0, D3=3D0, D4=4D0, D10=1D1) | |
2553 | PARAMETER (MXX = 105) | |
2554 | COMMON / LhCtXXARAY / XCR, XMIN, XV(0:MXX), LSTX, NX | |
2555 | DATA IWRN1, HUGE, TINY / 0, 1.E35, 1.E-35 / | |
2556 | F(X) = (XCR-XMIN) * LOG (X/XMIN) + LOG (XCR/XMIN) * (X-XMIN) | |
2557 | D(X) = (XCR-XMIN) / X + LOG (XCR/XMIN) | |
2558 | X = XX | |
2559 | IF (X .GE. XMIN) THEN | |
2560 | TEM = F(X) / F(D1) | |
2561 | ElseIF (X .GE. D0) THEN | |
2562 | X = MAX (X, TINY) | |
2563 | TEM = F(X) / F(D1) | |
2564 | Else | |
2565 | CALL CtLhWARNR(IWRN1, 'X out of range in CtLhZFRMX' | |
2566 | > , 'X', X, TINY, HUGE, 1) | |
2567 | TEM = 99. | |
2568 | STOP | |
2569 | EndIf | |
2570 | CtLhZFRMX = TEM | |
2571 | RETURN | |
2572 | ENTRY CtLhDZDX (XX) | |
2573 | X = XX | |
2574 | IF (X .GE. XMIN) THEN | |
2575 | TEM = D(X) / F(D1) | |
2576 | ElseIF (X .GE. D0) THEN | |
2577 | X = MAX (X, TINY) | |
2578 | TEM = D(X) / F(D1) | |
2579 | Else | |
2580 | CALL CtLhWARNR(IWRN1, 'X out of range in CtLhDZDX ' | |
2581 | > , 'X', X, TINY, HUGE, 1) | |
2582 | TEM = 99. | |
2583 | STOP | |
2584 | EndIf | |
2585 | CtLhDZDX = TEM | |
2586 | RETURN | |
2587 | END | |
2588 | FUNCTION CtLhZFXL (XL) | |
2589 | IMPLICIT DOUBLE PRECISION (A-H, O-Z) | |
2590 | COMMON / LhCtINVERT / ZA | |
2591 | X = EXP(XL) | |
2592 | TT = CtLhZFRMX (X) - ZA | |
2593 | CtLhZFXL = TT | |
2594 | RETURN | |
2595 | END | |
2596 |