]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | C----------------------------------------------------------------------- |
2 | #include "pdf/pilot.h" | |
3 | #if defined(CERNLIB_SINGLE) | |
4 | FUNCTION Ctq4PnX (IPRTN, X, Q) | |
5 | #endif | |
6 | #if defined(CERNLIB_DOUBLE) | |
7 | DOUBLE PRECISION FUNCTION Ctq4PnX (IPRTN, X, Q) | |
8 | #endif | |
9 | C | |
10 | C Given the parton distribution function in the array Upd in | |
11 | C COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of | |
12 | C x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda). | |
13 | C | |
14 | #include "pdf/impdp.inc" | |
15 | C | |
16 | PARAMETER (MXX = 105, MXQ = 25, MXF = 6) | |
17 | PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX) | |
18 | PARAMETER (M= 2, M1 = M + 1) | |
19 | C | |
20 | COMMON / W5051IC / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX) | |
21 | COMMON / W5051ID / Nx, Nt, NfMx | |
22 | COMMON / W5051IE / Qini, Qmax, Xmin | |
23 | C | |
24 | Dimension Fq(M1), Df(M1) | |
25 | C Work with Log (Q) | |
26 | QG = LOG (Q/AL) | |
27 | ||
28 | C Find lower end of interval containing X | |
29 | JL = -1 | |
30 | JU = Nx+1 | |
31 | 11 If (JU-JL .GT. 1) Then | |
32 | JM = (JU+JL) / 2 | |
33 | If (X .GT. XV(JM)) Then | |
34 | JL = JM | |
35 | Else | |
36 | JU = JM | |
37 | Endif | |
38 | Goto 11 | |
39 | Endif | |
40 | ||
41 | Jx = JL - (M-1)/2 | |
42 | If (X .lt. Xmin) Then | |
43 | C Print '(A, 2(1pE12.4))', | |
44 | C > ' WARNING: X < Xmin, extrapolation used; X, Xmin =', X, Xmin | |
45 | If (Jx .LT. 0) Jx = 0 | |
46 | Elseif (Jx .GT. Nx-M) Then | |
47 | Jx = Nx - M | |
48 | Endif | |
49 | C Find the interval where Q lies | |
50 | JL = -1 | |
51 | JU = NT+1 | |
52 | 12 If (JU-JL .GT. 1) Then | |
53 | JM = (JU+JL) / 2 | |
54 | If (QG .GT. QL(JM)) Then | |
55 | JL = JM | |
56 | Else | |
57 | JU = JM | |
58 | Endif | |
59 | Goto 12 | |
60 | Endif | |
61 | ||
62 | Jq = JL - (M-1)/2 | |
63 | If (Jq .LT. 0) Then | |
64 | Jq = 0 | |
65 | C If (Q .lt. Qini) Print '(A, 2(1pE12.4))', | |
66 | C > ' WARNING: Q < Qini, extrapolation used; Q, Qini =', Q, Qini | |
67 | Elseif (Jq .GT. Nt-M) Then | |
68 | Jq = Nt - M | |
69 | C If (Q .gt. Qmax) Print '(A, 2(1pE12.4))', | |
70 | C > ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax | |
71 | Endif | |
72 | ||
73 | If (Iprtn .GE. 3) Then | |
74 | Ip = - Iprtn | |
75 | Else | |
76 | Ip = Iprtn | |
77 | EndIf | |
78 | C Find the off-set in the linear array Upd | |
79 | JFL = Ip + NfMx | |
80 | J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx | |
81 | C | |
82 | C Now interpolate in x for M1 Q's | |
83 | Do 21 Iq = 1, M1 | |
84 | J1 = J0 + (Nx+1)*(Iq-1) + 1 | |
85 | Call Dpolin (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq)) | |
86 | 21 Continue | |
87 | C Finish off by interpolating in Q | |
88 | Call Dpolin (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf) | |
89 | ||
90 | Ctq4PnX = Ftmp | |
91 | C | |
92 | RETURN | |
93 | C **************************** | |
94 | END |