PDF/spdf/ctq4pnx.F

   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