partonx5.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file partonx5.F

      FUNCTION partonx5 (IPRTN, X, Q)
C
C   Given the parton distribution function in the array Upd in
C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
C
      IMPLICIT DOUBLE PRECISION (a-h, o-z)
C
      parameter(mxx = 105, mxq = 25, mxf = 6)
      parameter(mxpqx = (mxf *2 +2) * mxq * mxx)
      parameter(m= 2, m1 = m + 1)
C
      Logical first
      Common 
     > / ctqpar1 / al, xv(0:mxx), ql(0:mxq), upd(mxpqx)
     > / ctqpar2 / nx, nt, nfmx
     > / xqrange / qini, qmax, xmin
C
      dimension fq(m1), df(m1)

      Data first /.true./
      save first
C                                                 Work with Log (Q)
      qg  = log(q/al)

C                           Find lower end of interval containing X
      jl = -1
      ju = nx+1
 11   If (ju-jl .GT. 1) Then
         jm = (ju+jl) / 2
         If (x .GT. xv(jm)) Then
            jl = jm
         Else
            ju = jm
         Endif
         goto 11
      Endif

      jx = jl - (m-1)/2
      If (x .lt. xmin .and. first ) Then
         first = .false.
         print '(A, 2(1pE12.4))', 
     >     ' WARNING: X << Xmin, extrapolation used; X, Xmin =', x, xmin
         If (jx .LT. 0) jx = 0
      Elseif (jx .GT. nx-m) Then
         jx = nx - m
      Endif
C                                    Find the interval where Q lies
      jl = -1
      ju = nt+1
 12   If (ju-jl .GT. 1) Then
         jm = (ju+jl) / 2
         If (qg .GT. ql(jm)) Then
            jl = jm
         Else
            ju = jm
         Endif
         goto 12
      Endif

      jq = jl - (m-1)/2
      If (jq .LT. 0) Then
         jq = 0
         If (q .lt. qini)  print '(A, 2(1pE12.4))', 
     >     ' WARNING: Q << Qini, extrapolation used; Q, Qini =', q, qini
      Elseif (jq .GT. nt-m) Then
         jq = nt - m
         If (q .gt. qmax)  print '(A, 2(1pE12.4))', 
     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', q, qmax
      Endif

      If (iprtn .GE. 3) Then
         ip = - iprtn
      Else
         ip = iprtn
      EndIf
C                             Find the off-set in the linear array Upd
      jfl = ip + nfmx
      j0  = (jfl * (nt+1) + jq) * (nx+1) + jx
C
C                                           Now interpolate in x for M1 Q's
      Do 21 iq = 1, m1
         j1 = j0 + (nx+1)*(iq-1) + 1
         Call polint5(xv(jx), upd(j1), m1, x, fq(iq), df(iq))
 21   Continue
C                                          Finish off by interpolating in Q
      Call polint5(ql(jq), fq(1), m1, qg, ftmp, ddf)

      partonx5 = ftmp
C
      RETURN
C                        ****************************
      END