/* grv.f -- translated by f2c (version 19961001). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { integer iinip; } intinip_; #define intinip_1 intinip_ struct { integer iinif; } intinif_; #define intinif_1 intinif_ /* Table of constant values */ static integer c__1 = 1; static doublereal c_b32 = .5; static doublereal c_b33 = -.2; static integer c__2 = 2; static doublereal c_b71 = .2; /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* * */ /* G R V - P R O T O N - P A R A M E T R I Z A T I O N S * */ /* * */ /* 1998 UPDATE * */ /* * */ /* For a detailed explanation see * */ /* M. Glueck, E. Reya, A. Vogt : * */ /* hep-ph/9806404 = DO-TH 98/07 = WUE-ITP-98-019 * */ /* (To appear in Eur. Phys. J. C) * */ /* * */ /* This package contains subroutines returning the light-parton * */ /* distributions in NLO (for the MSbar and DIS schemes) and LO; * */ /* the respective light-parton, charm, and bottom contributions * */ /* to F2(electromagnetic); and the scale dependence of alpha_s. * */ /* * */ /* The parton densities and F2 values are calculated from inter- * */ /* polation grids covering the regions * */ /* Q^2/GeV^2 between 0.8 and 1.E6 ( 1.E4 for F2 ) * */ /* x between 1.E-9 and 1. * */ /* Any call outside these regions stops the program execution. * */ /* * */ /* At Q^2 = MZ^2, alpha_s reads 0.114 (0.125) in NLO (LO); the * */ /* heavy quark thresholds, QH^2 = mh^2, in the beta function are * */ /* mc = 1.4 GeV, mb = 4.5 GeV, mt = 175 GeV. * */ /* Note that the NLO alpha_s running is different from GRV(94). * */ /* * */ /* Questions, comments etc to: avogt@physik.uni-wuerzburg.de * */ /* * */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* Subroutine */ int grv98pa_(iset, x, q2, uv, dv, us, ds, ss, gl, iinit) integer *iset, *iinit; doublereal *x, *q2, *uv, *dv, *us, *ds, *ss, *gl; { /* Initialized data */ static doublereal qs[27] = { .8,1.,1.3,1.8,2.7,4.,6.4,10.,16.,25.,40.,64., 100.,180.,320.,570.,1e3,1800.,3200.,5700.,1e4,2.2e4,4.6e4,1e5, 2.2e5,4.6e5,1e6 }; static doublereal xb[68] = { 1e-9,1.8e-9,3.2e-9,5.7e-9,1e-8,1.8e-8,3.2e-8, 5.7e-8,1e-7,1.8e-7,3.2e-7,5.7e-7,1e-6,1.4e-6,2e-6,3e-6,4.5e-6, 6.7e-6,1e-5,1.4e-5,2e-5,3e-5,4.5e-5,6.7e-5,1e-4,1.4e-4,2e-4,3e-4, 4.5e-4,6.7e-4,.001,.0014,.002,.003,.0045,.0067,.01,.014,.02,.03, .045,.06,.08,.1,.125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375, .4,.45,.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1. }; /* Format strings */ static char fmt_91[] = "(2x,\002PARTON INTERPOLATION: X OUT OF RANGE\002)" ; static char fmt_92[] = "(2x,\002PARTON INTERPOLATION: Q2 OUT OF RANGE\ \002)"; static char fmt_93[] = "(2x,\002NO OR INVALID PARTON SET CHOICE\002)"; static char fmt_89[] = "(a80)"; static char fmt_90[] = "(6(1pe10.3))"; /* System generated locals */ doublereal d__1, d__2; olist o__1; cllist cl__1; /* Builtin functions */ integer s_wsfe(), e_wsfe(); /* Subroutine */ int s_stop(); integer f_open(), s_rsfe(), do_fio(), e_rsfe(), f_clos(); double pow_dd(), log(); /* Local variables */ static doublereal xdef[1836] /* was [68][27] */; static char line[80]; static doublereal arrf[95]; extern doublereal fint_(); static doublereal xudf[1836] /* was [68][27] */, xdvf[1836] /* was [68][27] */, xuvf[1836] /* was [68][27] */; static integer m, n; static doublereal x1, de; static integer na[2]; static doublereal ud; static integer iq, ix; static doublereal xt[2], xs, xv, parton[10854] /* was [6][27][67] */, xb1, xgf[1836] /* was [68][27] */, xsf[1836] /* was [68][ 27] */, xb0s, xb0v; /* Fortran I/O blocks */ static cilist io___3 = { 0, 6, 0, fmt_91, 0 }; static cilist io___4 = { 0, 6, 0, fmt_92, 0 }; static cilist io___5 = { 0, 6, 0, fmt_93, 0 }; static cilist io___6 = { 0, 11, 0, fmt_89, 0 }; static cilist io___10 = { 0, 11, 0, fmt_90, 0 }; /* ******************************************************************** */ /* * */ /* THE PARTON ROUTINE. * */ /* __ * */ /* INPUT: ISET = 1 (LO), 2 (NLO, MS), or 3 (NLO, DIS) * */ /* X = Bjorken-x (between 1.E-9 and 1.) * */ /* Q2 = scale in GeV**2 (between 0.8 and 1.E6) * */ /* * */ /* OUTPUT: UV = u - u(bar), DV = d - d(bar), US = u(bar), * */ /* DS = d(bar), SS = s = s(bar), GL = gluon. * */ /* Always x times the distribution is returned. * */ /* * */ /* COMMON: The main program or the calling routine has to have * */ /* a common block COMMON / INTINIP / IINIP , and the * */ /* integer variable IINIP has always to be zero when * */ /* GRV98PA is called for the first time or when ISET * */ /* has been changed. * */ /* * */ /* GRIDS: 1. grv98lo.grid, 2. grv98nlm.grid, 3. grv98nld.grid, * */ /* (1+1809 lines with 6 columns, 4 significant figures) * */ /* * */ /* ******************************************************i************* */ /* ...BJORKEN-X AND Q**2 VALUES OF THE GRID : */ /* printf("\n hello %lf %lf %d",*x,*q2, *iinit);*/ /* ...CHECK OF X AND Q2 VALUES : */ if (*x < 9.9e-10 || *x > 1.) { s_wsfe(&io___3); e_wsfe(); s_stop("", 0L); } if (*q2 < (float).799 || *q2 > (float)1.01e6) { s_wsfe(&io___4); e_wsfe(); s_stop("", 0L); } if (*iinit != 0) { goto L16; } /* ...INITIALIZATION, IF REQUIRED : */ /* SELECTION AND READING OF THE GRID : */ /* (COMMENT: FIRST NUMBER IN THE FIRST LINE OF THE GRID) */ if (*iset == 1) { o__1.oerr = 0; o__1.ounit = 11; o__1.ofnmlen = 12; o__1.ofnm = "grv98lo.grid"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* 7.332E-05 */ } else if (*iset == 2) { o__1.oerr = 0; o__1.ounit = 11; o__1.ofnmlen = 13; o__1.ofnm = "grv98nlm.grid"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* 1.015E-04 */ } else if (*iset == 3) { o__1.oerr = 0; o__1.ounit = 11; o__1.ofnmlen = 13; o__1.ofnm = "grv98nld.grid"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* 1.238E-04 */ } else { s_wsfe(&io___5); e_wsfe(); s_stop("", 0L); } /* iinit = 1;*/ s_rsfe(&io___6); do_fio(&c__1, line, 80L); e_rsfe(); for (m = 1; m <= 67; ++m) { for (n = 1; n <= 27; ++n) { s_rsfe(&io___10); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 168], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 167], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 166], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 165], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 164], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&parton[(n + m * 27) * 6 - 163], (ftnlen) sizeof(doublereal)); e_rsfe(); /* L15: */ } } cl__1.cerr = 0; cl__1.cunit = 11; cl__1.csta = 0; f_clos(&cl__1); /* ....ARRAYS FOR THE INTERPOLATION SUBROUTINE : */ for (iq = 1; iq <= 27; ++iq) { for (ix = 1; ix <= 67; ++ix) { xb0v = pow_dd(&xb[ix - 1], &c_b32); xb0s = pow_dd(&xb[ix - 1], &c_b33); xb1 = (float)1. - xb[ix - 1]; /* Computing 3rd power */ d__1 = xb1, d__2 = d__1; xuvf[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 168] / ( d__2 * (d__1 * d__1) * xb0v); /* Computing 4th power */ d__1 = xb1, d__1 *= d__1; xdvf[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 167] / ( d__1 * d__1 * xb0v); /* Computing 7th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; xdef[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 166] / ( d__2 * (d__1 * d__1) * xb0v); /* Computing 7th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; xudf[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 165] / ( d__2 * (d__1 * d__1) * xb0s); /* Computing 7th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; xsf[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 164] / (d__2 * (d__1 * d__1) * xb0s); /* Computing 5th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1; xgf[ix + iq * 68 - 69] = parton[(iq + ix * 27) * 6 - 163] / (d__2 * (d__1 * d__1) * xb0s); /* L20: */ } xuvf[iq * 68 - 1] = (float)0.; xdvf[iq * 68 - 1] = (float)0.; xdef[iq * 68 - 1] = (float)0.; xudf[iq * 68 - 1] = (float)0.; xsf[iq * 68 - 1] = (float)0.; xgf[iq * 68 - 1] = (float)0.; /* L10: */ } na[0] = 68; na[1] = 27; for (ix = 1; ix <= 68; ++ix) { arrf[ix - 1] = log(xb[ix - 1]); /* L30: */ } for (iq = 1; iq <= 27; ++iq) { arrf[iq + 67] = log(qs[iq - 1]); /* L40: */ } /* ...CONTINUATION, IF INITIALIZATION WAS DONE PREVIOUSLY. */ L16: /* ...INTERPOLATION : */ xt[0] = log(*x); xt[1] = log(*q2); x1 = (float)1. - *x; xv = pow_dd(x, &c_b32); xs = pow_dd(x, &c_b33); /* Computing 3rd power */ d__1 = x1, d__2 = d__1; *uv = fint_(&c__2, xt, na, arrf, xuvf) * (d__2 * (d__1 * d__1)) * xv; /* Computing 4th power */ d__1 = x1, d__1 *= d__1; *dv = fint_(&c__2, xt, na, arrf, xdvf) * (d__1 * d__1) * xv; /* Computing 7th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; de = fint_(&c__2, xt, na, arrf, xdef) * (d__2 * (d__1 * d__1)) * xv; /* Computing 7th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; ud = fint_(&c__2, xt, na, arrf, xudf) * (d__2 * (d__1 * d__1)) * xs; *us = (ud - de) * (float).5; *ds = (ud + de) * (float).5; /* Computing 7th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; *ss = fint_(&c__2, xt, na, arrf, xsf) * (d__2 * (d__1 * d__1)) * xs; /* Computing 5th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1; *gl = fint_(&c__2, xt, na, arrf, xgf) * (d__2 * (d__1 * d__1)) * xs; /* L60: */ return 0; } /* grv98pa_ */ /* Subroutine */ int grv98f2_(iset, x, q2, f2l, f2c, f2b, f2) integer *iset; doublereal *x, *q2, *f2l, *f2c, *f2b, *f2; { /* Initialized data */ static doublereal qs[21] = { .8,1.,1.3,1.8,2.7,4.,6.4,10.,16.,25.,40.,64., 100.,180.,320.,570.,1e3,1800.,3200.,5700.,1e4 }; static doublereal xb[68] = { 1e-9,1.8e-9,3.2e-9,5.7e-9,1e-8,1.8e-8,3.2e-8, 5.7e-8,1e-7,1.8e-7,3.2e-7,5.7e-7,1e-6,1.4e-6,2e-6,3e-6,4.5e-6, 6.7e-6,1e-5,1.4e-5,2e-5,3e-5,4.5e-5,6.7e-5,1e-4,1.4e-4,2e-4,3e-4, 4.5e-4,6.7e-4,.001,.0014,.002,.003,.0045,.0067,.01,.014,.02,.03, .045,.06,.08,.1,.125,.15,.175,.2,.225,.25,.275,.3,.325,.35,.375, .4,.45,.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1. }; /* Format strings */ static char fmt_91[] = "(2x,\002STR.FCT. INTERPOLATION: X OUT OF RANG\ E\002)"; static char fmt_92[] = "(2x,\002STR.FCT. INTERPOLATION: Q2 OUT OF RANG\ E\002)"; static char fmt_93[] = "(2x,\002NO OR INVALID STR.FCT. SET CHOICE\002)"; static char fmt_89[] = "(a80)"; static char fmt_90[] = "(3(1pe10.3))"; /* System generated locals */ doublereal d__1, d__2; olist o__1; cllist cl__1; /* Builtin functions */ integer s_wsfe(), e_wsfe(); /* Subroutine */ int s_stop(); integer f_open(), s_rsfe(), do_fio(), e_rsfe(), f_clos(); double pow_dd(), log(); /* Local variables */ static char line[80]; static doublereal arrf[89]; extern doublereal fint_(); static doublereal xf2bf[1428] /* was [68][21] */, xf2cf[1428] /* was [68][21] */, xf2lf[1428] /* was [68][21] */; static integer m, n; static doublereal x1; static integer na[2], iq, ix; static doublereal xs, xt[2], strfct[4221] /* was [3][21][67] */, xb1, xbs; /* Fortran I/O blocks */ static cilist io___33 = { 0, 6, 0, fmt_91, 0 }; static cilist io___34 = { 0, 6, 0, fmt_92, 0 }; static cilist io___35 = { 0, 6, 0, fmt_93, 0 }; static cilist io___36 = { 0, 11, 0, fmt_89, 0 }; static cilist io___40 = { 0, 11, 0, fmt_90, 0 }; /* ******************************************************************** */ /* * */ /* THE F2 ROUTINE. * */ /* * */ /* INPUT : ISET = 1 (LO), 2 (NLO). * */ /* X = Bjorken-x (between 1.E-9 and 1.) * */ /* Q2 = scale in GeV**2 (between 0.8 and 1.E4) * */ /* * */ /* OUTPUT: F2L = F2(light), F2C = F2(charm), F2B = F2(bottom,) * */ /* F2 = sum, all given for electromagnetic proton DIS. * */ /* * */ /* COMMON: The main program or the calling routine has to have * */ /* a common block COMMON / INTINIF / IINIF , and the * */ /* integer variable IINIF has always to be zero when * */ /* GRV98F2 is called for the first time or when ISET * */ /* has been changed. * */ /* * */ /* GRIDS: 1. grv98lof.grid, 2. grv98nlf.grid. * */ /* (1+1407 lines with 3 columns, 4 significant figures) * */ /* * */ /* ******************************************************i************* */ /* ...BJORKEN-X AND Q**2 VALUES OF THE GRID : */ /* ...CHECK OF X AND Q2 VALUES : */ if (*x < 9.9e-10 || *x > 1.) { s_wsfe(&io___33); e_wsfe(); s_stop("", 0L); } if (*q2 < (float).799 || *q2 > (float)10100.) { s_wsfe(&io___34); e_wsfe(); s_stop("", 0L); } if (intinif_1.iinif != 0) { goto L16; } /* ...INITIALIZATION, IF REQUIRED : */ /* SELECTION AND READING OF THE GRID : */ /* (COMMENT: FIRST NUMBER IN THE FIRST LINE OF THE GRID) */ if (*iset == 1) { o__1.oerr = 0; o__1.ounit = 11; o__1.ofnmlen = 13; o__1.ofnm = "grv98lof.grid"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* 7.907E-01 */ } else if (*iset == 2) { o__1.oerr = 0; o__1.ounit = 11; o__1.ofnmlen = 13; o__1.ofnm = "grv98nlf.grid"; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = 0; o__1.ofm = 0; o__1.oblnk = 0; f_open(&o__1); /* 9.368E-01 */ } else { s_wsfe(&io___35); e_wsfe(); s_stop("", 0L); } intinif_1.iinif = 1; s_rsfe(&io___36); do_fio(&c__1, line, 80L); e_rsfe(); for (m = 1; m <= 67; ++m) { for (n = 1; n <= 21; ++n) { s_rsfe(&io___40); do_fio(&c__1, (char *)&strfct[(n + m * 21) * 3 - 66], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&strfct[(n + m * 21) * 3 - 65], (ftnlen) sizeof(doublereal)); do_fio(&c__1, (char *)&strfct[(n + m * 21) * 3 - 64], (ftnlen) sizeof(doublereal)); e_rsfe(); /* L15: */ } } cl__1.cerr = 0; cl__1.cunit = 11; cl__1.csta = 0; f_clos(&cl__1); /* ....ARRAYS FOR THE INTERPOLATION SUBROUTINE : */ for (iq = 1; iq <= 21; ++iq) { for (ix = 1; ix <= 67; ++ix) { xbs = pow_dd(&xb[ix - 1], &c_b71); xb1 = (float)1. - xb[ix - 1]; /* Computing 3rd power */ d__1 = xb1, d__2 = d__1; xf2lf[ix + iq * 68 - 69] = strfct[(iq + ix * 21) * 3 - 66] / ( d__2 * (d__1 * d__1)) * xbs; /* Computing 7th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; xf2cf[ix + iq * 68 - 69] = strfct[(iq + ix * 21) * 3 - 65] / ( d__2 * (d__1 * d__1)) * xbs; /* Computing 7th power */ d__1 = xb1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; xf2bf[ix + iq * 68 - 69] = strfct[(iq + ix * 21) * 3 - 64] / ( d__2 * (d__1 * d__1)) * xbs; /* L20: */ } xf2lf[iq * 68 - 1] = (float)0.; xf2cf[iq * 68 - 1] = (float)0.; xf2bf[iq * 68 - 1] = (float)0.; /* L10: */ } na[0] = 68; na[1] = 21; for (ix = 1; ix <= 68; ++ix) { arrf[ix - 1] = log(xb[ix - 1]); /* L30: */ } for (iq = 1; iq <= 21; ++iq) { arrf[iq + 67] = log(qs[iq - 1]); /* L40: */ } /* ...CONTINUATION, IF INITIALIZATION WAS DONE PREVIOUSLY. */ L16: /* ...INTERPOLATION : */ xt[0] = log(*x); xt[1] = log(*q2); x1 = (float)1. - *x; xs = pow_dd(x, &c_b33); /* Computing 3rd power */ d__1 = x1, d__2 = d__1; *f2l = fint_(&c__2, xt, na, arrf, xf2lf) * (d__2 * (d__1 * d__1)) * xs; /* Computing 7th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; *f2c = fint_(&c__2, xt, na, arrf, xf2cf) * (d__2 * (d__1 * d__1)) * xs; /* Computing 7th power */ d__1 = x1, d__2 = d__1, d__1 *= d__1, d__2 *= d__1; *f2b = fint_(&c__2, xt, na, arrf, xf2bf) * (d__2 * (d__1 * d__1)) * xs; *f2 = *f2l + *f2c + *f2b; /* L60: */ return 0; } /* grv98f2_ */ doublereal fint_(narg, arg, nent, ent, table) integer *narg; doublereal *arg; integer *nent; doublereal *ent, *table; { /* System generated locals */ integer i__1, i__2; doublereal ret_val; /* Local variables */ static integer iadr, ient[5]; static doublereal d__[5]; static integer i__, j, k, m, ifadr, ncomb[5], ja, jb, kd, il, jr; static doublereal fac; /* ******************************************************************** */ /* * */ /* THE INTERPOLATION ROUTINE (CERN LIBRARY ROUTINE E104) * */ /* * */ /* ******************************************************************** */ /* Parameter adjustments */ --table; --ent; --nent; --arg; /* Function Body */ kd = 1; m = 1; ja = 1; i__1 = *narg; for (i__ = 1; i__ <= i__1; ++i__) { ncomb[i__ - 1] = 1; jb = ja - 1 + nent[i__]; i__2 = jb; for (j = ja; j <= i__2; ++j) { if (arg[i__] <= ent[j]) { goto L3; } /* L2: */ } j = jb; L3: if (j != ja) { goto L4; } ++j; L4: jr = j - 1; d__[i__ - 1] = (ent[j] - arg[i__]) / (ent[j] - ent[jr]); ient[i__ - 1] = j - ja; kd += ient[i__ - 1] * m; m *= nent[i__]; /* L5: */ ja = jb + 1; } ret_val = (float)0.; L10: fac = (float)1.; iadr = kd; ifadr = 1; i__1 = *narg; for (i__ = 1; i__ <= i__1; ++i__) { if (ncomb[i__ - 1] == 0) { goto L12; } fac *= (float)1. - d__[i__ - 1]; goto L15; L12: fac *= d__[i__ - 1]; iadr -= ifadr; L15: ifadr *= nent[i__]; } ret_val += fac * table[iadr]; il = *narg; L40: if (ncomb[il - 1] == 0) { goto L80; } ncomb[il - 1] = 0; if (il == *narg) { goto L10; } ++il; i__1 = *narg; for (k = il; k <= i__1; ++k) { /* L50: */ ncomb[k - 1] = 1; } goto L10; L80: --il; if (il != 0) { goto L40; } return ret_val; } /* fint_ */ doublereal alphas_(q2, naord) doublereal *q2; integer *naord; { /* Initialized data */ static doublereal q2thr[3] = { 1.96,20.25,30625. }; static doublereal lambdal[4] = { .2041,.175,.132,.0665 }; static doublereal lambdan[4] = { .2994,.246,.1677,.0678 }; /* Format strings */ static char fmt_91[] = "(\002INVALID CHOICE FOR ORDER IN ALPHA_S\002)"; /* System generated locals */ doublereal ret_val; /* Builtin functions */ double log(); integer s_wsfe(), e_wsfe(); /* Subroutine */ int s_stop(); /* Local variables */ static integer i__, k; static doublereal y, b0, b1, y1, b10; static integer nf; static doublereal xl, lq2, alp, xlm, xlp, lam2; /* Fortran I/O blocks */ static cilist io___79 = { 0, 6, 0, fmt_91, 0 }; /* ******************************************************************** */ /* * */ /* THE ALPHA_S ROUTINE. * */ /* * */ /* INPUT : Q2 = scale in GeV**2 (not too low, of course); * */ /* NAORD = 1 (LO), 2 (NLO). * */ /* * */ /* OUTPUT: alphas_s/(4 pi) for use with the GRV(98) partons. * */ /* * */ /* ******************************************************i************* */ /* ...HEAVY QUARK THRESHOLDS AND LAMBDA VALUES : */ /* ...DETERMINATION OF THE APPROPRIATE NUMBER OF FLAVOURS : */ nf = 3; for (k = 1; k <= 3; ++k) { if (*q2 > q2thr[k - 1]) { ++nf; } else { goto L20; } /* L10: */ } /* ...LO ALPHA_S AND BETA FUNCTION FOR NLO CALCULATION : */ L20: b0 = (float)11. - nf * (float).66666666666666663; b1 = (float)102. - nf * (float)12.666666666666666; b10 = b1 / (b0 * b0); if (*naord == 1) { lam2 = lambdal[nf - 3] * lambdal[nf - 3]; alp = (float)1. / (b0 * log(*q2 / lam2)); goto L1; } else if (*naord == 2) { lam2 = lambdan[nf - 3] * lambdan[nf - 3]; b1 = (float)102. - nf * (float)12.666666666666666; b10 = b1 / (b0 * b0); } else { s_wsfe(&io___79); e_wsfe(); s_stop("", 0L); } /* ...START VALUE FOR NLO ITERATION : */ lq2 = log(*q2 / lam2); alp = (float)1. / (b0 * lq2) * ((float)1. - b10 * log(lq2) / lq2); /* ...EXACT NLO VALUE, FOUND VIA NEWTON PROCEDURE : */ for (i__ = 1; i__ <= 6; ++i__) { xl = log((float)1. / (b0 * alp) + b10); xlp = log((float)1. / (b0 * alp * (float)1.01) + b10); xlm = log((float)1. / (b0 * alp * (float).99) + b10); y = lq2 - (float)1. / (b0 * alp) + b10 * xl; y1 = ((float)-1. / (b0 * alp * (float)1.01) + b10 * xlp + (float)1. / (b0 * alp * (float).99) - b10 * xlp) / (alp * .02); alp -= y / y1; /* L2: */ } /* ...OUTPUT : */ L1: ret_val = alp; return ret_val; } /* alphas_ */