/* Ctq5Pdf.f -- translated by f2c (version 20000121). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { integer nx, nt, nfmx; } ctqpar2_; #define ctqpar2_1 ctqpar2_ struct { doublereal alambda; integer nfl, iorder; } qcdtable_; #define qcdtable_1 qcdtable_ struct { doublereal al, xv[106], ql[26], upd[36750]; } ctqpar1_; #define ctqpar1_1 ctqpar1_ struct { doublereal qini, qmax, xmin; } xqrange_; #define xqrange_1 xqrange_ struct { doublereal amass[6]; } masstbl_; #define masstbl_1 masstbl_ /* Table of constant values */ static integer c__9 = 9; static integer c__1 = 1; static integer c__5 = 5; static integer c__3 = 3; /* ============================================================================ */ /* CTEQ Parton Distribution Functions: Version 5.0 */ /* March 1, 1999 */ /* Ref: "GLOBAL QCD ANALYSIS OF PARTON STRUCTURE OF THE NUCLEON: */ /* CTEQ5 PPARTON DISTRIBUTIONS" */ /* hep-ph/9903282 */ /* These PDF's use quadratic interpolation of attached tables. A parametrized */ /* version of the same PDF's without external tables is under construction. */ /* They will become available later. */ /* This package contains 7 sets of CTEQ5 PDF's. Details are: */ /* --------------------------------------------------------------------------- */ /* Iset PDF Description Alpha_s(Mz) Lam4 Lam5 Table_File */ /* --------------------------------------------------------------------------- */ /* 1 CTEQ5M Standard MSbar scheme 0.118 326 226 cteq5m.tbl */ /* 2 CTEQ5D Standard DIS scheme 0.118 326 226 cteq5d.tbl */ /* 3 CTEQ5L Leading Order 0.127 192 146 cteq5l.tbl */ /* 4 CTEQ5HJ Large-x gluon enhanced 0.118 326 226 cteq5hj.tbl */ /* 5 CTEQ5HQ Heavy Quark 0.118 326 226 cteq5hq.tbl */ /* 6 CTEQ5F3 Nf=3 FixedFlavorNumber 0.106 (Lam3=395) cteq5f3.tbl */ /* 7 CTEQ5F4 Nf=4 FixedFlavorNumber 0.112 309 XXX cteq5f4.tbl */ /* --------------------------------------------------------------------------- */ /* The available applied range is 10^-5 < x < 1 and 1.0 < Q < 10,000 (GeV). */ /* Lam5 (Lam4, Lam3) represents Lambda value (in MeV) for 5 (4,3) flavors. */ /* The matching alpha_s between 4 and 5 flavors takes place at Q=4.5 GeV, */ /* which is defined as the bottom quark mass, whenever it can be applied. */ /* The Table_Files are assumed to be in the working directory. */ /* Before using the PDF, it is necessary to do the initialization by */ /* Call SetCtq5(Iset) */ /* where Iset is the desired PDF specified in the above table. */ /* The function Ctq5Pdf (Iparton, X, Q) */ /* returns the parton distribution inside the proton for parton [Iparton] */ /* at [X] Bjorken_X and scale [Q] (GeV) in PDF set [Iset]. */ /* Iparton is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5) */ /* for (b, c, s, d, u, g, u_bar, ..., b_bar), */ /* whereas CTEQ5F3 has, by definition, only 3 flavors and gluon; */ /* CTEQ5F4 has only 4 flavors and gluon. */ /* For detailed information on the parameters used, e.q. quark masses, */ /* QCD Lambda, ... etc., see info lines at the beginning of the */ /* Table_Files. */ /* These programs, as provided, are in double precision. By removing the */ /* "Implicit Double Precision" lines, they can also be run in single */ /* precision. */ /* If you have detailed questions concerning these CTEQ5 distributions, */ /* or if you find problems/bugs using this package, direct inquires to */ /* Hung-Liang Lai(lai@phys.nthu.edu.tw) or Wu-Ki Tung(Tung@pa.msu.edu). */ /* =========================================================================== */ doublereal ctq5pdf_(iparton, x, q) integer *iparton; doublereal *x, *q; { /* Initialized data */ static logical warn = TRUE_; /* System generated locals */ doublereal ret_val; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(); /* Subroutine */ int s_stop(); /* Local variables */ extern doublereal partonx_(); /* Fortran I/O blocks */ static cilist io___2 = { 0, 6, 0, 0, 0 }; static cilist io___3 = { 0, 6, 0, 0, 0 }; static cilist io___4 = { 0, 6, 0, 0, 0 }; if (*x < 0. || *x > 1.) { s_wsle(&io___2); do_lio(&c__9, &c__1, "X out of range in Ctq5Pdf: ", (ftnlen)27); do_lio(&c__5, &c__1, (char *)&(*x), (ftnlen)sizeof(doublereal)); e_wsle(); s_stop("", (ftnlen)0); } if (*q < qcdtable_1.alambda) { s_wsle(&io___3); do_lio(&c__9, &c__1, "Q out of range in Ctq5Pdf: ", (ftnlen)27); do_lio(&c__5, &c__1, (char *)&(*q), (ftnlen)sizeof(doublereal)); e_wsle(); s_stop("", (ftnlen)0); } if (*iparton < -ctqpar2_1.nfmx || *iparton > ctqpar2_1.nfmx) { if (warn) { /* put a warning for calling extra flavor. */ warn = FALSE_; s_wsle(&io___4); do_lio(&c__9, &c__1, "Warning: Iparton out of range in Ctq5Pdf: ", (ftnlen)42); do_lio(&c__3, &c__1, (char *)&(*iparton), (ftnlen)sizeof(integer)) ; e_wsle(); } ret_val = 0.; return ret_val; } ret_val = partonx_(iparton, x, q); if (ret_val < 0.) { ret_val = 0.; } return ret_val; /* ******************** */ } /* ctq5pdf_ */ doublereal partonx_(iprtn, x, q) integer *iprtn; doublereal *x, *q; { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ doublereal ret_val; /* Builtin functions */ double log(); integer s_wsfe(), do_fio(), e_wsfe(); /* Local variables */ static doublereal ftmp; static integer j0, j1; static doublereal df[3]; static integer jl; static doublereal fq[3], qg; static integer jm, ip, jq, iq, ju, jx; extern /* Subroutine */ int polint_(); static doublereal ddf; static integer jfl; /* Fortran I/O blocks */ static cilist io___11 = { 0, 6, 0, "(A, 2(1pE12.4))", 0 }; static cilist io___13 = { 0, 6, 0, "(A, 2(1pE12.4))", 0 }; static cilist io___14 = { 0, 6, 0, "(A, 2(1pE12.4))", 0 }; /* Given the parton distribution function in the array Upd in */ /* COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of */ /* x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda). */ /* Work with Log (Q) */ qg = log(*q / ctqpar1_1.al); /* Find lower end of interval containing X */ jl = -1; ju = ctqpar2_1.nx + 1; L11: if (ju - jl > 1) { jm = (ju + jl) / 2; if (*x > ctqpar1_1.xv[jm]) { jl = jm; } else { ju = jm; } goto L11; } jx = jl; if (*x < xqrange_1.xmin && first) { first = FALSE_; s_wsfe(&io___11); do_fio(&c__1, " WARNING: X < Xmin, extrapolation used; X, Xmin =", ( ftnlen)49); do_fio(&c__1, (char *)&(*x), (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&xqrange_1.xmin, (ftnlen)sizeof(doublereal)); e_wsfe(); if (jx < 0) { jx = 0; } } else if (jx > ctqpar2_1.nx - 2) { jx = ctqpar2_1.nx - 2; } /* Find the interval where Q lies */ jl = -1; ju = ctqpar2_1.nt + 1; L12: if (ju - jl > 1) { jm = (ju + jl) / 2; if (qg > ctqpar1_1.ql[jm]) { jl = jm; } else { ju = jm; } goto L12; } jq = jl; if (jq < 0) { jq = 0; if (*q < xqrange_1.qini) { s_wsfe(&io___13); do_fio(&c__1, " WARNING: Q < Qini, extrapolation used; Q, Qini =", (ftnlen)49); do_fio(&c__1, (char *)&(*q), (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&xqrange_1.qini, (ftnlen)sizeof(doublereal)) ; e_wsfe(); } } else if (jq > ctqpar2_1.nt - 2) { jq = ctqpar2_1.nt - 2; if (*q > xqrange_1.qmax) { s_wsfe(&io___14); do_fio(&c__1, " WARNING: Q > Qmax, extrapolation used; Q, Qmax =", (ftnlen)49); do_fio(&c__1, (char *)&(*q), (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&xqrange_1.qmax, (ftnlen)sizeof(doublereal)) ; e_wsfe(); } } if (*iprtn >= 3) { ip = -(*iprtn); } else { ip = *iprtn; } /* Find the off-set in the linear array Upd */ jfl = ip + ctqpar2_1.nfmx; j0 = (jfl * (ctqpar2_1.nt + 1) + jq) * (ctqpar2_1.nx + 1) + jx; /* Now interpolate in x for M1 Q's */ for (iq = 1; iq <= 3; ++iq) { j1 = j0 + (ctqpar2_1.nx + 1) * (iq - 1) + 1; polint_(&ctqpar1_1.xv[jx], &ctqpar1_1.upd[j1 - 1], &c__3, x, &fq[iq - 1], &df[iq - 1]); /* L21: */ } /* Finish off by interpolating in Q */ polint_(&ctqpar1_1.ql[jq], fq, &c__3, &qg, &ftmp, &ddf); ret_val = ftmp; return ret_val; /* **************************** */ } /* partonx_ */ /* Subroutine */ int setctq5_(iset) integer *iset; { /* Initialized data */ static char flnm[12*7+1] = "cteq5m.tbl cteq5d.tbl cteq5l.tbl cteq5hj.\ tbl cteq5hq.tbl cteq5f3.tbl cteq5f4.tbl "; static char tablefile[40+1] = "test.tbl "; static integer isetold = -987; static integer isetmin = 1; static integer isettest = 911; /* System generated locals */ integer i__1; olist o__1; cllist cl__1; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(), f_open(); /* Subroutine */ int s_stop(), s_copy(); integer f_clos(), s_rsfe(), do_fio(), e_rsfe(); /* Local variables */ static integer iu; extern integer nextun_(); extern /* Subroutine */ int readtbl_(); /* Fortran I/O blocks */ static cilist io___30 = { 0, 6, 0, 0, 0 }; static cilist io___31 = { 0, 6, 0, 0, 0 }; static cilist io___32 = { 0, 6, 0, 0, 0 }; static cilist io___33 = { 0, 6, 0, 0, 0 }; static cilist io___34 = { 0, 6, 0, 0, 0 }; static cilist io___35 = { 0, 5, 0, "(A)", 0 }; /* If data file not initialized, do so. */ if (*iset != isetold) { iu = nextun_(); if (*iset == isettest) { s_wsle(&io___30); do_lio(&c__9, &c__1, "Opening ", (ftnlen)8); do_lio(&c__9, &c__1, tablefile, (ftnlen)40); e_wsle(); L21: o__1.oerr = 1; o__1.ounit = iu; o__1.ofnmlen = 40; o__1.ofnm = tablefile; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L101; } } else if (*iset < isetmin || *iset > 7) { s_wsle(&io___31); do_lio(&c__9, &c__1, "Invalid Iset number in SetCtq5 :", (ftnlen) 32); do_lio(&c__3, &c__1, (char *)&(*iset), (ftnlen)sizeof(integer)); e_wsle(); s_stop("", (ftnlen)0); } else { s_copy(tablefile, flnm + (*iset - 1) * 12, (ftnlen)40, (ftnlen)12) ; o__1.oerr = 1; o__1.ounit = iu; o__1.ofnmlen = 40; o__1.ofnm = tablefile; o__1.orl = 0; o__1.osta = "OLD"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; i__1 = f_open(&o__1); if (i__1 != 0) { goto L100; } } readtbl_(&iu); cl__1.cerr = 0; cl__1.cunit = iu; cl__1.csta = 0; f_clos(&cl__1); isetold = *iset; } return 0; L100: s_wsle(&io___32); do_lio(&c__9, &c__1, " Data file ", (ftnlen)11); do_lio(&c__9, &c__1, tablefile, (ftnlen)40); do_lio(&c__9, &c__1, " cannot be opened in SetCtq5!!", (ftnlen)30); e_wsle(); s_stop("", (ftnlen)0); L101: s_wsle(&io___33); do_lio(&c__9, &c__1, tablefile, (ftnlen)40); do_lio(&c__9, &c__1, " cannot be opened ", (ftnlen)18); e_wsle(); s_wsle(&io___34); do_lio(&c__9, &c__1, "Please input the .tbl file:", (ftnlen)27); e_wsle(); s_rsfe(&io___35); do_fio(&c__1, tablefile, (ftnlen)40); e_rsfe(); goto L21; /* ******************** */ } /* setctq5_ */ /* Subroutine */ int readtbl_(nu) integer *nu; { /* System generated locals */ integer i__1; /* Builtin functions */ integer s_rsfe(), do_fio(), e_rsfe(), s_rsle(), do_lio(), e_rsle(), i_dnnt(); double log(); /* Local variables */ static integer nblk; static char line[80]; static integer iret, npts, i__; static doublereal fl, dr; static integer iq; /* Fortran I/O blocks */ static cilist io___36 = { 0, 0, 0, "(A)", 0 }; static cilist io___38 = { 0, 0, 0, "(A)", 0 }; static cilist io___39 = { 0, 0, 0, 0, 0 }; static cilist io___43 = { 0, 0, 0, "(A)", 0 }; static cilist io___44 = { 0, 0, 0, 0, 0 }; static cilist io___45 = { 0, 0, 0, "(A)", 0 }; static cilist io___46 = { 0, 0, 0, 0, 0 }; static cilist io___47 = { 0, 0, 0, "(A)", 0 }; static cilist io___48 = { 0, 0, 0, 0, 0 }; static cilist io___52 = { 0, 0, 0, "(A)", 0 }; static cilist io___54 = { 1, 0, 1, 0, 0 }; io___36.ciunit = *nu; s_rsfe(&io___36); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___38.ciunit = *nu; s_rsfe(&io___38); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___39.ciunit = *nu; s_rsle(&io___39); do_lio(&c__5, &c__1, (char *)&dr, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&fl, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&ctqpar1_1.al, (ftnlen)sizeof(doublereal)); for (i__ = 1; i__ <= 6; ++i__) { do_lio(&c__5, &c__1, (char *)&masstbl_1.amass[i__ - 1], (ftnlen) sizeof(doublereal)); } e_rsle(); qcdtable_1.iorder = i_dnnt(&dr); qcdtable_1.nfl = i_dnnt(&fl); qcdtable_1.alambda = ctqpar1_1.al; io___43.ciunit = *nu; s_rsfe(&io___43); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___44.ciunit = *nu; s_rsle(&io___44); do_lio(&c__3, &c__1, (char *)&ctqpar2_1.nx, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&ctqpar2_1.nt, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&ctqpar2_1.nfmx, (ftnlen)sizeof(integer)); e_rsle(); io___45.ciunit = *nu; s_rsfe(&io___45); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___46.ciunit = *nu; s_rsle(&io___46); do_lio(&c__5, &c__1, (char *)&xqrange_1.qini, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&xqrange_1.qmax, (ftnlen)sizeof(doublereal)); i__1 = ctqpar2_1.nt; for (i__ = 0; i__ <= i__1; ++i__) { do_lio(&c__5, &c__1, (char *)&ctqpar1_1.ql[i__], (ftnlen)sizeof( doublereal)); } e_rsle(); io___47.ciunit = *nu; s_rsfe(&io___47); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___48.ciunit = *nu; s_rsle(&io___48); do_lio(&c__5, &c__1, (char *)&xqrange_1.xmin, (ftnlen)sizeof(doublereal)); i__1 = ctqpar2_1.nx; for (i__ = 0; i__ <= i__1; ++i__) { do_lio(&c__5, &c__1, (char *)&ctqpar1_1.xv[i__], (ftnlen)sizeof( doublereal)); } e_rsle(); i__1 = ctqpar2_1.nt; for (iq = 0; iq <= i__1; ++iq) { ctqpar1_1.ql[iq] = log(ctqpar1_1.ql[iq] / ctqpar1_1.al); /* L11: */ } /* Since quark = anti-quark for nfl>2 at this stage, */ /* we Read out only the non-redundent data points */ /* No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) */ nblk = (ctqpar2_1.nx + 1) * (ctqpar2_1.nt + 1); npts = nblk * (ctqpar2_1.nfmx + 3); io___52.ciunit = *nu; s_rsfe(&io___52); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___54.ciunit = *nu; iret = s_rsle(&io___54); if (iret != 0) { goto L100001; } i__1 = npts; for (i__ = 1; i__ <= i__1; ++i__) { iret = do_lio(&c__5, &c__1, (char *)&ctqpar1_1.upd[i__ - 1], (ftnlen) sizeof(doublereal)); if (iret != 0) { goto L100001; } } iret = e_rsle(); L100001: return 0; /* **************************** */ } /* readtbl_ */ integer nextun_() { /* System generated locals */ integer ret_val; inlist ioin__1; /* Builtin functions */ integer f_inqu(); /* Subroutine */ int s_stop(); /* Local variables */ static integer n; static logical ex; /* Returns an unallocated FORTRAN i/o unit. */ for (n = 10; n <= 300; ++n) { ioin__1.inerr = 0; ioin__1.inunit = n; ioin__1.infile = 0; ioin__1.inex = 0; ioin__1.inopen = &ex; ioin__1.innum = 0; ioin__1.innamed = 0; ioin__1.inname = 0; ioin__1.inacc = 0; ioin__1.inseq = 0; ioin__1.indir = 0; ioin__1.infmt = 0; ioin__1.inform = 0; ioin__1.inunf = 0; ioin__1.inrecl = 0; ioin__1.innrec = 0; ioin__1.inblank = 0; f_inqu(&ioin__1); if (! ex) { ret_val = n; return ret_val; } /* L10: */ } s_stop(" There is no available I/O unit. ", (ftnlen)33); /* ************************* */ return ret_val; } /* nextun_ */ /* Subroutine */ int polint_(xa, ya, n, x, y, dy) doublereal *xa, *ya; integer *n; doublereal *x, *y, *dy; { /* System generated locals */ integer i__1, i__2; doublereal d__1; /* Builtin functions */ /* Subroutine */ int s_paus(); /* Local variables */ static doublereal dift, c__[10], d__[10]; static integer i__, m; static doublereal w, ho, hp; static integer ns; static doublereal dif, den; /* Adapted from "Numerical Recipes" */ /* Parameter adjustments */ --ya; --xa; /* Function Body */ ns = 1; dif = (d__1 = *x - xa[1], abs(d__1)); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { dift = (d__1 = *x - xa[i__], abs(d__1)); if (dift < dif) { ns = i__; dif = dift; } c__[i__ - 1] = ya[i__]; d__[i__ - 1] = ya[i__]; /* L11: */ } *y = ya[ns]; --ns; i__1 = *n - 1; for (m = 1; m <= i__1; ++m) { i__2 = *n - m; for (i__ = 1; i__ <= i__2; ++i__) { ho = xa[i__] - *x; hp = xa[i__ + m] - *x; w = c__[i__] - d__[i__ - 1]; den = ho - hp; if (den == (float)0.) { s_paus("", (ftnlen)0); } den = w / den; d__[i__ - 1] = hp * den; c__[i__ - 1] = ho * den; /* L12: */ } if (ns << 1 < *n - m) { *dy = c__[ns]; } else { *dy = d__[ns - 1]; --ns; } *y += *dy; /* L13: */ } return 0; } /* polint_ */