Functions/Subroutines
function	ctq6pdf (Iparton, X, Q)

subroutine	setctq6 (Iset)

subroutine	readtbl (Nu)

function	nextun ()

subroutine	polint (XA, YA, N, X, Y, DY)

function	partonx6 (IPRTN, XX, QQ)

Function/Subroutine Documentation

◆ ctq6pdf()

function ctq6pdf	(	Iparton,
		X,
		Q
	)

Definition at line 93 of file Cteq61Pdf.f.

       Implicit Double Precision (a-h,o-z)
       Logical Warn
       Common
      > / ctqpar2 / nx, nt, nfmx
      > / qcdtable /  alambda, nfl, iorder
  
       Data warn /.true./
       save warn
  
       If (x .lt. 0d0 .or. x .gt. 1d0) Then
         print *, 'X out of range in Ctq6Pdf: ', x
         stop
       Endif
       If (q .lt. alambda) Then
         print *, 'Q out of range in Ctq6Pdf: ', q
         stop
       Endif
       If ((iparton .lt. -nfmx .or. iparton .gt. nfmx)) Then
          If (warn) Then
 C        put a warning for calling extra flavor.
              print *, 'Warning: Iparton out of range in Ctq6Pdf! '
              print *, 'Iparton, MxFlvN0: ', iparton, nfmx
          Endif
          ctq6pdf = 0d0
          Return
       Endif
  
       ctq6pdf = partonx6(iparton, x, q)
       if(ctq6pdf.lt.0.d0)  ctq6pdf = 0.d0
       if(isnan(ctq6pdf)) then
          if (warn) then
             print *, 'Warning: Ctq6Pdf is NaN!'
             print *, 'Iparton, MxFlvN0, X, Q: ', iparton, nfmx, x, q
          endif
          ctq6pdf = 0.d0
       endif
  
       if (warn) warn = .false.
  
       Return
  
 C                             ********************

◆ nextun()

function nextun

Definition at line 224 of file Cteq61Pdf.f.

 C                                 Returns an unallocated FORTRAN i/o unit.
       Logical EX
 C
       Do 10 n = 10, 300
          INQUIRE (unit=n, opened=ex)
          If (.NOT. ex) then
             nextun = n
             Return
          Endif
  10   Continue
       stop ' There is no available I/O unit. '
 C               *************************

◆ partonx6()

function partonx6	(	IPRTN,
		XX,
		QQ
	)

Definition at line 281 of file Cteq61Pdf.f.

  
 c  Given the parton distribution function in the array U in
 c  COMMON / PEVLDT / , this routine interpolates to find
 c  the parton distribution at an arbitray point in x and q.
 c
       Implicit Double Precision (a-h,o-z)
  
       parameter(mxx = 96, mxq = 20, mxf = 5)
       parameter(mxqx= mxq * mxx,   mxpqx = mxqx * (mxf+3))
  
       Common
      > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
      > / ctqpar2 / nx, nt, nfmx
      > / xqrange / qini, qmax, xmin
  
       dimension fvec(4), fij(4)
       dimension xvpow(0:mxx)
       Data onep / 1.00001 /
       Data xpow / 0.3d0 /       !**** choice of interpolation variable
       Data nqvec / 4 /
       Data ientry / 0 /
       Save ientry,xvpow
  
 c store the powers used for interpolation on first call...
       if(ientry .eq. 0) then
          ientry = 1
  
          xvpow(0) = 0d0
          do i = 1, nx
             xvpow(i) = xv(i)**xpow
          enddo
       endif
  
       x = xx
       q = qq
       tt = log(log(q/al))
  
 c      -------------    find lower end of interval containing x, i.e.,
 c                       get jx such that xv(jx) .le. x .le. xv(jx+1)...
       jlx = -1
       ju = nx+1
  11   If (ju-jlx .GT. 1) Then
          jm = (ju+jlx) / 2
          If (x .Ge. xv(jm)) Then
             jlx = jm
          Else
             ju = jm
          Endif
          Goto 11
       Endif
 C                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
 C                           |---|---|---|...|---|-x-|---|...|---|---|
 C                     x     0  Xmin               x                 1
 C
       If     (jlx .LE. -1) Then
         print '(A,1pE12.4)', 'Severe error: x <= 0 in PartonX6! x = ', x
         stop
       ElseIf (jlx .Eq. 0) Then
          jx = 0
       Elseif (jlx .LE. nx-2) Then
  
 C                For interrior points, keep x in the middle, as shown above
          jx = jlx - 1
       Elseif (jlx.Eq.nx-1 .or. x.LT.onep) Then
  
 C                  We tolerate a slight over-shoot of one (OneP=1.00001),
 C              perhaps due to roundoff or whatever, but not more than that.
 C                                      Keep at least 4 points >= Jx
          jx = jlx - 2
       Else
         print '(A,1pE12.4)', 'Severe error: x > 1 in PartonX6! x = ', x
         stop
       Endif
 C          ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
  
 C                       This is the variable to be interpolated in
       ss = x**xpow
  
       If (jlx.Ge.2 .and. jlx.Le.nx-2) Then
  
 c     initiation work for "interior bins": store the lattice points in s...
       svec1 = xvpow(jx)
       svec2 = xvpow(jx+1)
       svec3 = xvpow(jx+2)
       svec4 = xvpow(jx+3)
  
       s12 = svec1 - svec2
       s13 = svec1 - svec3
       s23 = svec2 - svec3
       s24 = svec2 - svec4
       s34 = svec3 - svec4
  
       sy2 = ss - svec2
       sy3 = ss - svec3
  
 c constants needed for interpolating in s at fixed t lattice points...
       const1 = s13/s23
       const2 = s12/s23
       const3 = s34/s23
       const4 = s24/s23
       s1213 = s12 + s13
       s2434 = s24 + s34
       sdet = s12*s34 - s1213*s2434
       tmp = sy2*sy3/sdet
       const5 = (s34*sy2-s2434*sy3)*tmp/s12
       const6 = (s1213*sy2-s12*sy3)*tmp/s34
  
       EndIf
  
 c         --------------Now find lower end of interval containing Q, i.e.,
 c                          get jq such that qv(jq) .le. q .le. qv(jq+1)...
       jlq = -1
       ju = nt+1
  12   If (ju-jlq .GT. 1) Then
          jm = (ju+jlq) / 2
          If (tt .GE. tv(jm)) Then
             jlq = jm
          Else
             ju = jm
          Endif
          Goto 12
        Endif
  
       If     (jlq .LE. 0) Then
          jq = 0
       Elseif (jlq .LE. nt-2) Then
 C                                  keep q in the middle, as shown above
          jq = jlq - 1
       Else
 C                         JLq .GE. Nt-1 case:  Keep at least 4 points >= Jq.
         jq = nt - 3
  
       Endif
 C                                   This is the interpolation variable in Q
  
       If (jlq.GE.1 .and. jlq.LE.nt-2) Then
 c                                        store the lattice points in t...
       tvec1 = tv(jq)
       tvec2 = tv(jq+1)
       tvec3 = tv(jq+2)
       tvec4 = tv(jq+3)
  
       t12 = tvec1 - tvec2
       t13 = tvec1 - tvec3
       t23 = tvec2 - tvec3
       t24 = tvec2 - tvec4
       t34 = tvec3 - tvec4
  
       ty2 = tt - tvec2
       ty3 = tt - tvec3
  
       tmp1 = t12 + t13
       tmp2 = t24 + t34
  
       tdet = t12*t34 - tmp1*tmp2
  
       EndIf
  
  
 c get the pdf function values at the lattice points...
  
       If (iprtn .GE. 3) Then
          ip = - iprtn
       Else
          ip = iprtn
       EndIf
       jtmp = ((ip + nfmx)*(nt+1)+(jq-1))*(nx+1)+jx+1
  
       Do it = 1, nqvec
  
          j1  = jtmp + it*(nx+1)
  
        If (jx .Eq. 0) Then
 C                          For the first 4 x points, interpolate x^2*f(x,Q)
 C                           This applies to the two lowest bins JLx = 0, 1
 C            We can not put the JLx.eq.1 bin into the "interrior" section
 C                           (as we do for q), since Upd(J1) is undefined.
          fij(1) = 0
          fij(2) = upd(j1+1) * xv(1)**2
          fij(3) = upd(j1+2) * xv(2)**2
          fij(4) = upd(j1+3) * xv(3)**2
 C
 C                 Use Polint which allows x to be anywhere w.r.t. the grid
  
          Call polint (xvpow(0), fij(1), 4, ss, fx, dfx)
  
          If (x .GT. 0d0)  fvec(it) =  fx / x**2
 C                                              Pdf is undefined for x.eq.0
        ElseIf  (jlx .Eq. nx-1) Then
 C                                                This is the highest x bin:
  
         Call polint (xvpow(nx-3), upd(j1), 4, ss, fx, dfx)
  
         fvec(it) = fx
  
        Else
 C                       for all interior points, use Jon's in-line function
 C                              This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
          sf2 = upd(j1+1)
          sf3 = upd(j1+2)
  
          g1 =  sf2*const1 - sf3*const2
          g4 = -sf2*const3 + sf3*const4
  
          fvec(it) = (const5*(upd(j1)-g1)
      &               + const6*(upd(j1+3)-g4)
      &               + sf2*sy3 - sf3*sy2) / s23
  
        Endif
  
       enddo
 C                                   We now have the four values Fvec(1:4)
 c     interpolate in t...
  
       If (jlq .LE. 0) Then
 C                         1st Q-bin, as well as extrapolation to lower Q
         Call polint (tv(0), fvec(1), 4, tt, ff, dfq)
  
       ElseIf (jlq .GE. nt-1) Then
 C                         Last Q-bin, as well as extrapolation to higher Q
         Call polint (tv(nt-3), fvec(1), 4, tt, ff, dfq)
       Else
 C                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
 C       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
 C                         the full range QV(0:Nt)  (in contrast to XV)
         tf2 = fvec(2)
         tf3 = fvec(3)
  
         g1 = ( tf2*t13 - tf3*t12) / t23
         g4 = (-tf2*t34 + tf3*t24) / t23
  
         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
  
         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
       EndIf
  
       partonx6 = ff
  
       Return
 C                                       ********************

◆ polint()

subroutine polint	(	dimension(n)	XA,
		dimension(n)	YA,
			N,
			X,
			Y,
			DY
	)

Definition at line 240 of file Cteq61Pdf.f.

  
       IMPLICIT DOUBLE PRECISION (a-h, o-z)
 C                                        Adapted from "Numerical Recipes"
       parameter(nmax=10)
       dimension xa(n),ya(n),c(nmax),d(nmax)
       ns=1
       dif=abs(x-xa(1))
       DO 11 i=1,n
         dift=abs(x-xa(i))
         IF (dift.LT.dif) THEN
           ns=i
           dif=dift
         ENDIF
         c(i)=ya(i)
         d(i)=ya(i)
 11    CONTINUE
       y=ya(ns)
       ns=ns-1
       DO 13 m=1,n-1
         DO 12 i=1,n-m
           ho=xa(i)-x
           hp=xa(i+m)-x
           w=c(i+1)-d(i)
           den=ho-hp
 !           IF(DEN.EQ.0.)PAUSE
           den=w/den
           d(i)=hp*den
           c(i)=ho*den
 12      CONTINUE
         IF (2*ns.LT.n-m)THEN
           dy=c(ns+1)
         ELSE
           dy=d(ns)
           ns=ns-1
         ENDIF
         y=y+dy
 13    CONTINUE
       RETURN

◆ readtbl()

subroutine readtbl ( Nu )

Definition at line 179 of file Cteq61Pdf.f.

       Implicit Double Precision (a-h,o-z)
       Character Line*80
       parameter(mxx = 96, mxq = 20, mxf = 5)
       parameter(mxpqx = (mxf + 3) * mxq * mxx)
       Common
      > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
      > / ctqpar2 / nx, nt, nfmx
      > / xqrange / qini, qmax, xmin
      > / qcdtable /  alambda, nfl, iorder
      > / masstbl / amass(6)
  
       Read  (nu, '(A)') line
       Read  (nu, '(A)') line
       Read  (nu, *) dr, fl, al, (amass(i),i=1,6)
       iorder = nint(dr)
       nfl = nint(fl)
       alambda = al
  
       Read  (nu, '(A)') line
       Read  (nu, *) nx,  nt, nfmx
  
       Read  (nu, '(A)') line
       Read  (nu, *) qini, qmax, (tv(i), i =0, nt)
  
       Read  (nu, '(A)') line
       Read  (nu, *) xmin, (xv(i), i =0, nx)
  
       Do 11 iq = 0, nt
          tv(iq) = log(log(tv(iq) /al))
    11 Continue
 C
 C                  Since quark = anti-quark for nfl>2 at this stage,
 C                  we Read  out only the non-redundent data points
 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
  
       nblk = (nx+1) * (nt+1)
       npts =  nblk  * (nfmx+3)
       Read  (nu, '(A)') line
       Read  (nu, *, iostat=iret) (upd(i), i=1,npts)
  
       Return
 C                        ****************************

◆ setctq6()

subroutine setctq6 ( Iset )

Definition at line 138 of file Cteq61Pdf.f.

       Implicit Double Precision (a-h,o-z)
       parameter(isetmax0=5)
       Character Flnm(Isetmax0)*6, nn*3, Tablefile*40
       Data (flnm(i), i=1,isetmax0)
      > / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l','ctq61.'/
       Data isetold, isetmin0, isetmin1, isetmax1 /-987,1,100,140/
       Data isetmin2,isetmax2 /200,240/
       save
  
 C             If data file not initialized, do so.
       If(iset.ne.isetold) then
          iu= nextun()
          If (iset.ge.isetmin0 .and. iset.le.3) Then
             tablefile=flnm(iset)//'.tbl'
          Elseif (iset.eq.4) Then
             tablefile=flnm(iset)//'1.tbl'
          Elseif (iset.ge.isetmin1 .and. iset.le.isetmax1) Then
             write(nn,'(I3)') iset
             tablefile=flnm(1)//nn//'.tbl'
          Elseif (iset.ge.isetmin2 .and. iset.le.isetmax2) Then
             write(nn,'(I3)') iset
             tablefile=flnm(5)//nn(2:3)//'.tbl'
          Else
             print *, 'Invalid Iset number in SetCtq6 :', iset
             stop
          Endif
          Open(iu, file='pdfs/'//tablefile
      .            ,status='OLD', err=100)
  21      Call readtbl (iu)
          Close (iu)
          isetold=iset
       Endif
       Return
  
  100  print *, ' Data file ', tablefile, ' cannot be opened '
      >//'in SetCtq6!!'
       stop
 C                             ********************