/* Ctq4Fn.f -- translated by f2c (version 19980913). 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 4.3 */ /* June 21, 1996 */ /* Modified on October 17, 1996 */ /* Modified on January 7, 1997 */ /* Modified on January 15, 1997 */ /* Ref[1]: "IMPROVED PARTON DISTRIBUTIONS FROM GLOBAL ANALYSIS OF RECENT DEEP */ /* INELASTIC SCATTERING AND INCLUSIVE JET DATA" */ /* By: H.L. Lai, J. Huston, S. Kuhlmann, F. Olness, J. Owens, D. Soper */ /* W.K. Tung, H. Weerts */ /* MSUHEP-60426, CTEQ-604, e-Print Archive: hep-ph/9606399, */ /* to appear in Phys. Rev. D55 (1997) */ /* Ref[2]: "CHARM PRODUCTION AND PARTON DISTRIBUTIONS" */ /* By: H.L. Lai and W.K. Tung */ /* MSU-HEP-61222, CTEQ-622, e-Print Archive: hep-ph/9701256 */ /* This package contains 12 sets of CTEQ4 PDF's. Details are: */ /* --------------------------------------------------------------------------- */ /* Iset PDF Description Alpha_s(Mz) Lam4 Lam5 Table_File */ /* --------------------------------------------------------------------------- */ /* 1 CTEQ4M Standard MSbar scheme 0.116 298 202 cteq4m.tbl */ /* 2 CTEQ4D Standard DIS scheme 0.116 298 202 cteq4d.tbl */ /* 3 CTEQ4L Leading Order 0.132 236 181 cteq4l.tbl */ /* 4 CTEQ4A1 Alpha_s series 0.110 215 140 cteq4a1.tbl */ /* 5 CTEQ4A2 Alpha_s series 0.113 254 169 cteq4a2.tbl */ /* 6 CTEQ4A3 ( same as CTEQ4M ) */ /* 7 CTEQ4A4 Alpha_s series 0.119 346 239 cteq4a4.tbl */ /* 8 CTEQ4A5 Alpha_s series 0.122 401 282 cteq4a5.tbl */ /* 9 CTEQ4HJ High Jet 0.116 303 206 cteq4hj.tbl */ /* 10 CTEQ4LQ Low Q0 0.114 261 174 cteq4lq.tbl */ /* --------------------------------------------------------------------------- */ /* 11 CTEQ4HQ Heavy Quark 0.116 298 202 cteq4hq.tbl */ /* 12 CTEQ4F3 Nf=3 FFN 0.106 (Lam3=385) cteq4f3.tbl */ /* 13 CTEQ4F4 Nf=4 FFN 0.110 282 XXX cteq4f4.tbl */ /* --------------------------------------------------------------------------- */ /* The available applied range is 10^-5 < x < 1 and 1.6 < Q < 10,000 (GeV) */ /* except CTEQ4LQ for which Q starts at a lower value of 0.7 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=5.0 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 SetCtq4(Iset) */ /* where Iset is the desired PDF specified in the above table. */ /* The function Ctq4Pdf (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 CTEQ4F3 has, by definition, only 3 flavors and gluon; */ /* CTEQ4F4 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 CTEQ4 distributions, */ /* or if you find problems/bugs using this package, direct inquires to */ /* Hung-Liang Lai(Lai_H@pa.msu.edu) or Wu-Ki Tung(Tung@pa.msu.edu). */ /* =========================================================================== */ doublereal ctq4pdf_(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 partonx4_(); /* 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 Ctq4Pdf: ", (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 Ctq4Pdf: ", (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 Ctq4Pdf: ", (ftnlen)42); do_lio(&c__3, &c__1, (char *)&(*iparton), (ftnlen)sizeof(integer)) ; e_wsle(); } ret_val = 0.; return ret_val; } ret_val = partonx4_(iparton, x, q); if (ret_val < 0.) { ret_val = 0.; } return ret_val; /* ******************** */ } /* ctq4pdf_ */ doublereal partonx4_(iprtn, x, q) integer *iprtn; doublereal *x, *q; { /* 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 polint4_(); static doublereal ddf; static integer jfl; /* Fortran I/O blocks */ static cilist io___10 = { 0, 6, 0, "(A, 2(1pE12.4))", 0 }; static cilist io___12 = { 0, 6, 0, "(A, 2(1pE12.4))", 0 }; static cilist io___13 = { 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) { s_wsfe(&io___10); 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___12); 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___13); 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; polint4_(&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 */ polint4_(&ctqpar1_1.ql[jq], fq, &c__3, &qg, &ftmp, &ddf); ret_val = ftmp; return ret_val; /* **************************** */ } /* partonx4_ */ /* Subroutine */ int setctq4_(iset) integer *iset; { /* Initialized data */ static char flnm[11*13+1] = "cteq4m.tbl cteq4d.tbl cteq4l.tbl cteq4a1.tb\ lcteq4a2.tblcteq4m.tbl cteq4a4.tblcteq4a5.tblcteq4hj.tblcteq4lq.tblcteq4hq.t\ blcteq4f3.tblcteq4f4.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 nextun4_(); extern /* Subroutine */ int readtbl4_(); /* Fortran I/O blocks */ static cilist io___29 = { 0, 6, 0, 0, 0 }; 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, 5, 0, "(A)", 0 }; /* If data file not initialized, do so. */ if (*iset != isetold) { iu = nextun4_(); if (*iset == isettest) { s_wsle(&io___29); 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 > 13) { s_wsle(&io___30); do_lio(&c__9, &c__1, "Invalid Iset number in SetCtq4 :", (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) * 11, (ftnlen)40, (ftnlen)11) ; 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; } } readtbl4_(&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___31); 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 SetCtq4!!", (ftnlen)30); e_wsle(); s_stop("", (ftnlen)0); L101: s_wsle(&io___32); 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___33); do_lio(&c__9, &c__1, "Please input the .tbl file:", (ftnlen)27); e_wsle(); s_rsfe(&io___34); do_fio(&c__1, tablefile, (ftnlen)40); e_rsfe(); goto L21; /* ******************** */ } /* setctq4_ */ /* Subroutine */ int readtbl4_(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___35 = { 0, 0, 0, "(A)", 0 }; static cilist io___37 = { 0, 0, 0, "(A)", 0 }; static cilist io___38 = { 0, 0, 0, 0, 0 }; static cilist io___42 = { 0, 0, 0, "(A)", 0 }; static cilist io___43 = { 0, 0, 0, 0, 0 }; static cilist io___44 = { 0, 0, 0, "(A)", 0 }; static cilist io___45 = { 0, 0, 0, 0, 0 }; static cilist io___46 = { 0, 0, 0, "(A)", 0 }; static cilist io___47 = { 0, 0, 0, 0, 0 }; static cilist io___51 = { 0, 0, 0, "(A)", 0 }; static cilist io___53 = { 1, 0, 1, 0, 0 }; io___35.ciunit = *nu; s_rsfe(&io___35); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___37.ciunit = *nu; s_rsfe(&io___37); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___38.ciunit = *nu; s_rsle(&io___38); 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___42.ciunit = *nu; s_rsfe(&io___42); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___43.ciunit = *nu; s_rsle(&io___43); 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___44.ciunit = *nu; s_rsfe(&io___44); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___45.ciunit = *nu; s_rsle(&io___45); 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___46.ciunit = *nu; s_rsfe(&io___46); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___47.ciunit = *nu; s_rsle(&io___47); 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___51.ciunit = *nu; s_rsfe(&io___51); do_fio(&c__1, line, (ftnlen)80); e_rsfe(); io___53.ciunit = *nu; iret = s_rsle(&io___53); 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; /* **************************** */ } /* readtbl4_ */ integer nextun4_() { /* 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; } /* nextun4_ */ /* Subroutine */ int polint4_(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; } /* polint4_ */