/* mrsg.f -- translated by f2c (version 19961001).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Table of constant values */

static integer c__1 = 1;

/* Subroutine */ int mrseb_(x, scale, mode, upv, dnv, usea, dsea, str, chm, 
	bot, glu)
doublereal *x, *scale;
integer *mode;
doublereal *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Format strings */
    static char fmt_99[] = "(\002  WARNING:  Q^2 VALUE IS OUT OF RANGE   \
\002)";
    static char fmt_98[] = "(\002  WARNING:   X  VALUE IS OUT OF RANGE   \
\002)";

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int strc30_(), strc31_();
    static doublereal q2;

    /* Fortran I/O blocks */
    static cilist io___2 = { 0, 6, 0, fmt_99, 0 };
    static cilist io___3 = { 0, 6, 0, fmt_98, 0 };


/* ***************************************************************C */
/* 								C */
/*     This is a package for the new MRS(A',G) parton            C */
/*     distributions. The minimum Q^2  value is 5 GeV^2 and the  C */
/*     x range is, as before 10^-5 < x < 1. MSbar factorization  C */
/*     is used. The package reads 2 grids, which are in separate C */
/*     files (A'=for030.dat/ftn30, G=for031.dat/ftn31).          C */
/*     Note that x times the parton distribution is returned,    C */
/*     Q is the scale in GeV,                                    C */
/*     and Lambda(MSbar,nf=4) = 231/255 MeV for A'/G.            C */
/* 								C */
/* 	MODE=1 for MRS(A')                                      C */
/* 	MODE=2 for MRS(G)                                       C */
/* 								C */
/*         The reference is :                                    C */
/*         A.D. Martin, R.G. Roberts and W.J. Stirling,          C */
/*         Phys. Lett. B354 (1995) 155-162                       C */
/*                                                               C */
/*         Comments to : W.J.Stirling@durham.ac.uk               C */
/*                                                               C */
/*             >>>>>>>>  CROSS CHECK  <<<<<<<<                   C */
/*                                                               C */
/*         THE FIRST NUMBER IN THE 30 GRID IS 0.00341            C */
/*         THE FIRST NUMBER IN THE 31 GRID IS 0.00269            C */
/* 								C */
/* ***************************************************************C */
/* Computing 2nd power */
    d__1 = *scale;
    q2 = d__1 * d__1;
    if (q2 < 5. || q2 > 1310720.) {
	s_wsfe(&io___2);
	e_wsfe();
    }
    if (*x < 1e-5 || *x > 1.) {
	s_wsfe(&io___3);
	e_wsfe();
    }
    if (*mode == 1) {
	strc30_(x, scale, upv, dnv, usea, dsea, str, chm, bot, glu);
    }
    if (*mode == 2) {
	strc31_(x, scale, upv, dnv, usea, dsea, str, chm, bot, glu);
    }
    return 0;
} /* mrseb_ */


/* Subroutine */ int strc30_(x, scale, upv, dnv, usea, dsea, str, chm, bot, 
	glu)
doublereal *x, *scale, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static integer init = 0;
    static doublereal xx[47] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.02,.04,.06,.08,.1,.125,.15,
	    .175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,.5,.525,
	    .55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 5.;
    static doublereal qsqmax = 1310720.;
    static integer n0[8] = { 2,5,5,9,0,0,9,9 };

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    integer s_rsfe(), do_fio(), e_rsfe();
    double pow_di(), d_lg10(), log(), d_int(), pow_dd();

    /* Local variables */
    static doublereal a, b, f[7520]	/* was [8][47][20] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, rm, fac, qsq, xxx;

    /* Fortran I/O blocks */
    static cilist io___14 = { 0, 30, 0, fmt_50, 0 };


/*     THIS IS THE NEW  "Aprime" FIT -- Feb 1995 -- standard Q^2 range */
    xsave = *x;
    if (init != 0) {
	goto L10;
    }
    init = 1;
    for (n = 1; n <= 46; ++n) {
	for (m = 1; m <= 19; ++m) {
	    s_rsfe(&io___14);
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 384], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 383], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 382], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 381], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 380], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 378], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 379], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 377], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* 1=UV 2=DV 3=GLUE 4=UBAR 5=CBAR 7=BBAR 6=SBAR 8=DBAR */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 47 << 3) - 385] /= pow_di(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 20; ++j) {
	xx[j - 1] = d_lg10(&xx[j - 1]) + 1.1;
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 19; ++k) {
/* L30: */
		f[i__ + (j + k * 47 << 3) - 385] = log(f[i__ + (j + k * 47 << 
			3) - 385]) * f[i__ + (k * 47 + 21 << 3) - 385] / log(
			f[i__ + (k * 47 + 21 << 3) - 385]);
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 19; ++m) {
/* L40: */
	    f[i__ + (m * 47 + 47 << 3) - 385] = 0.;
	}
    }
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
/* Computing 2nd power */
    d__1 = *scale;
    qsq = d__1 * d__1;
    if (qsq < qsqmin) {
	qsq = qsqmin;
    }
    if (qsq > qsqmax) {
	qsq = qsqmax;
    }
    xxx = *x;
    if (*x < .1) {
	xxx = d_lg10(x) + 1.1;
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    rm = log(qsq / qsqmin) / log(2.);
    b = rm - d_int(&rm);
    m = (integer) rm + 1;
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 47 << 3) - 385] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 47 << 3) - 385] + a * (
		1. - b) * f[i__ + (n + 1 + m * 47 << 3) - 385] + a * b * f[
		i__ + (n + 1 + (m + 1) * 47 << 3) - 385];
	if (n >= 21) {
	    goto L65;
	}
	if (i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 47 + 21 << 3) - 385] + b * f[i__ + ((m 
		+ 1) * 47 + 21 << 3) - 385];
	d__1 = g[i__ - 1] / fac;
	g[i__ - 1] = pow_dd(&fac, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_di(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    return 0;
} /* strc30_ */


/* Subroutine */ int strc31_(x, scale, upv, dnv, usea, dsea, str, chm, bot, 
	glu)
doublereal *x, *scale, *upv, *dnv, *usea, *dsea, *str, *chm, *bot, *glu;
{
    /* Initialized data */

    static doublereal xx[47] = { 1e-5,2e-5,4e-5,6e-5,8e-5,1e-4,2e-4,4e-4,6e-4,
	    8e-4,.001,.002,.004,.006,.008,.01,.02,.04,.06,.08,.1,.125,.15,
	    .175,.2,.225,.25,.275,.3,.325,.35,.375,.4,.425,.45,.475,.5,.525,
	    .55,.575,.6,.65,.7,.75,.8,.9,1. };
    static doublereal xmin = 1e-5;
    static doublereal xmax = 1.;
    static doublereal qsqmin = 5.;
    static doublereal qsqmax = 1310720.;
    static integer n0[8] = { 2,5,5,9,0,0,9,9 };
    static integer init = 0;

    /* Format strings */
    static char fmt_50[] = "(8f10.5)";

    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    integer s_rsfe(), do_fio(), e_rsfe();
    double pow_di(), d_lg10(), log(), d_int(), pow_dd();

    /* Local variables */
    static doublereal a, b, f[7520]	/* was [8][47][20] */, g[8];
    static integer i__, j, k, m, n;
    static doublereal xsave, rm, fac, qsq, xxx;

    /* Fortran I/O blocks */
    static cilist io___36 = { 0, 31, 0, fmt_50, 0 };


/*     THIS IS THE NEW  "G" FIT -- Feb 1995 -- standard Q^2 range */
    xsave = *x;
    if (init != 0) {
	goto L10;
    }
    init = 1;
    for (n = 1; n <= 46; ++n) {
	for (m = 1; m <= 19; ++m) {
	    s_rsfe(&io___36);
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 384], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 383], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 382], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 381], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 380], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 378], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 379], (ftnlen)sizeof(
		    doublereal));
	    do_fio(&c__1, (char *)&f[(n + m * 47 << 3) - 377], (ftnlen)sizeof(
		    doublereal));
	    e_rsfe();
/* 1=UV 2=DV 3=GLUE 4=UBAR 5=CBAR 7=BBAR 6=SBAR 8=DBAR */
	    for (i__ = 1; i__ <= 8; ++i__) {
/* L25: */
		d__1 = 1. - xx[n - 1];
		f[i__ + (n + m * 47 << 3) - 385] /= pow_di(&d__1, &n0[i__ - 1]
			);
	    }
/* L20: */
	}
    }
    for (j = 1; j <= 20; ++j) {
	xx[j - 1] = d_lg10(&xx[j - 1]) + 1.1;
	for (i__ = 1; i__ <= 8; ++i__) {
	    if (i__ == 7) {
		goto L31;
	    }
	    for (k = 1; k <= 19; ++k) {
/* L30: */
		f[i__ + (j + k * 47 << 3) - 385] = log(f[i__ + (j + k * 47 << 
			3) - 385]) * f[i__ + (k * 47 + 21 << 3) - 385] / log(
			f[i__ + (k * 47 + 21 << 3) - 385]);
	    }
L31:
	    ;
	}
    }
    for (i__ = 1; i__ <= 8; ++i__) {
	for (m = 1; m <= 19; ++m) {
/* L40: */
	    f[i__ + (m * 47 + 47 << 3) - 385] = 0.;
	}
    }
L10:
    if (*x < xmin) {
	*x = xmin;
    }
    if (*x > xmax) {
	*x = xmax;
    }
/* Computing 2nd power */
    d__1 = *scale;
    qsq = d__1 * d__1;
    if (qsq < qsqmin) {
	qsq = qsqmin;
    }
    if (qsq > qsqmax) {
	qsq = qsqmax;
    }
    xxx = *x;
    if (*x < .1) {
	xxx = d_lg10(x) + 1.1;
    }
    n = 0;
L70:
    ++n;
    if (xxx > xx[n]) {
	goto L70;
    }
    a = (xxx - xx[n - 1]) / (xx[n] - xx[n - 1]);
    rm = log(qsq / qsqmin) / log(2.);
    b = rm - d_int(&rm);
    m = (integer) rm + 1;
    for (i__ = 1; i__ <= 8; ++i__) {
	g[i__ - 1] = (1. - a) * (1. - b) * f[i__ + (n + m * 47 << 3) - 385] + 
		(1. - a) * b * f[i__ + (n + (m + 1) * 47 << 3) - 385] + a * (
		1. - b) * f[i__ + (n + 1 + m * 47 << 3) - 385] + a * b * f[
		i__ + (n + 1 + (m + 1) * 47 << 3) - 385];
	if (n >= 21) {
	    goto L65;
	}
	if (i__ == 7) {
	    goto L65;
	}
	fac = (1. - b) * f[i__ + (m * 47 + 21 << 3) - 385] + b * f[i__ + ((m 
		+ 1) * 47 + 21 << 3) - 385];
	d__1 = g[i__ - 1] / fac;
	g[i__ - 1] = pow_dd(&fac, &d__1);
L65:
	d__1 = 1. - *x;
	g[i__ - 1] *= pow_di(&d__1, &n0[i__ - 1]);
/* L60: */
    }
    *upv = g[0];
    *dnv = g[1];
    *usea = g[3];
    *dsea = g[7];
    *str = g[5];
    *chm = g[4];
    *glu = g[2];
    *bot = g[6];
    *x = xsave;
    return 0;
} /* strc31_ */