#include <math.h>
#include <stdio.h>

main(argc,argv)
     int argc;
     char *argv[];
{
  double  scale, upv, dnv, usea, dsea, str, chm, bot, glu;
  double  x, dbar, ubar, gluon, u, d, s, c, b, t;
  double  xmin,xmax, dx, u_int, d_int;
  int mrseb(), nx, n;

  
  /* Pick a range of *x* values to plot */
  xmin =   0.0001;
  xmax =   0.9999;
  nx   =  200;

  if (argc != 2) {
    printf("Usage: %s <scale (in GeV)>\n",argv[0]);
    exit(1);
  }
  sscanf(argv[1],"%lf",&scale);
  
  /* While looping on *x* we check if 
     
     Integral[ U(x)/x, {x,0,1} ] = 2
     Integral[ D(x)/x, {x,0,1} ] = 1
     
     as appropriate for a proton 
     */
  
  u_int = d_int = 0.0;
  
  /* Loop on *x* */
  dx = (xmax-xmin)/nx;
  x  = xmin - dx;
  for (n=0; n<nx; n++) {
    x += dx;
    mrseb(x,scale,&upv,&dnv,&usea,&dsea,&str,&chm,&bot,&glu);
    dbar  =  dsea;
    ubar  =  usea;
    gluon =  glu;
    u     =  upv+usea;
    d     =  dnv+dsea;
    s     =  str;
    c     =  chm;
    b     =  bot;
    printf("%lf %lf\n",x,upv);
    
    /* Do simple Riemann sum for the integrals */
    u_int += (upv/x)*dx;
    d_int += (dnv/x)*dx;
  }
  /*  printf("Up integral: %lf    Down integral: %lf\n",u_int,d_int);*/
}


/* C translation of MRS routine, which came with the comment ...

  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  c                                                                    c
  c  This program returns the fitted parton distributions for MRS(A')  c
  c  ... it is an analytic approximation to the "grid version" and     c
  c  gives much smoother distributions. It is, however, much slower.   c
  c  Full details of the fit are contained in the preprint "Pinning    c
  c  Down the Gluon in the Proton", by A.D. Martin, R.G. Roberts and   c
  c  W.J. Stirling, Phys. Lett. B354 (1995) 155-162,  to which         c
  c  reference should be made if appropriate.                          c  
  c                                                                    c
  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

*/

/* Table of fit parameters for parton distributions */
static double z[8][18] = {
  2.385035, -1.750356,  0.370491,  0.581567,  0.033300,  0.012245,  /* U Valence */
  3.925280, -0.226187, -0.062425, -0.862021,  1.178930,  0.038682,
  4.894751, -1.377950,  0.280011,  0.061201,  0.749858,  0.092008,

  0.270620, -0.195488,  0.040248,  0.357825,  0.156226, -0.020230,  /* D Valence */
  4.339434, -0.525551, -0.016963,  6.583414, -0.366604,  1.082098,
  2.174049,  0.263928, -1.336412,  0.101026,  1.140400, -0.026956,
  
  0.852604, -0.683552,  0.214748, -0.105954, -0.211982,  0.010073,  /* Glue */
  4.305152,  0.992946,  0.019575, -1.042052,  0.027041, -0.255378,
  4.441180, -1.800216,  0.129922,  0.639220,  0.833309, -0.341173,

  0.328730, -0.173874,  0.006988, -0.228390,  0.025648, -0.083929,  /* U+D Sea */
  9.748830, -0.727386,  0.397301, -1.384145, -0.538920,  0.432677,
  10.440065, -4.280664, 1.752318, -0.192712,  0.860909, -0.191209,

  0.012332, -0.002796,  0.000000, -0.205990,  0.123876,  0.070738,  /* U-D Sea */
  10.590572, 0.702136, -0.391293,  0.000000, 16.654520,  8.695223,  
  190.873402,-92.774861,-17.858415, -1.078797,0.650973,  0.080318,

  0.042355, -0.025763,  0.006208, -0.117133, -0.090933, -0.031820,  /* Strange */
  8.486509, -1.178987,  0.706619, -1.478869,  0.416645,  0.202853,
  12.931922, -4.845551,  0.221202,  0.563591, 0.954011, -0.329131,

  0.001518,  0.018562, -0.005142, -0.067847, -0.212202,  0.017218,  /* Charm */ 
  5.612756,  3.786456, -2.084374, -0.140671, -3.278622,  1.815593,
  6.910856,  0.785627, -3.351047,  1.384535, -0.069658, -0.078327,

  -0.004515, 0.013875, -0.004472, -0.146553, -0.034326, -0.043137,  /* Bottom */
  2.593083,  3.836694, -0.763188,  1.095127, -2.988716,  1.042826,
  -2.628259,  7.263671, -2.674676, 1.660073,  0.206697, -0.081092};

/* These data come a fit using Lamba_MS = 231 MeV */
#define LAMBDA 0.231

mrseb(x,scale,upv,dnv,usea,dsea,str,chm,bot,glu)
     double x,scale;               /* Inputs */
     double *upv,*dnv,*usea,*dsea; /* Outputs */
     double *str,*chm,*bot,*glu;
{
  extern double z[8][18];
  double q2, q0, arg,s;
  double a,b,c,d,e,h,dusum,dudif;
  double f[8];
  int j;

  q2 = scale*scale;
  q0 = 2.0;
  arg = log(scale/LAMBDA)/log(q0/LAMBDA);
  s = log(arg);

  for (j=0; j<8; j++) {
    a = z[j][ 0] + z[j][ 1]*s + z[j][ 2]*s*s;
    b = z[j][ 3] + z[j][ 4]*s + z[j][ 5]*s*s;
    c = z[j][ 6] + z[j][ 7]*s + z[j][ 8]*s*s;
    d = z[j][ 9] + z[j][10]*s + z[j][11]*s*s;
    e = z[j][12] + z[j][13]*s + z[j][14]*s*s;
    h = z[j][15] + z[j][16]*s + z[j][17]*s*s;
    f[j] = a*pow(x,b) * pow((1.-x),c) * (1.+d*sqrt(x)+e*x)* pow(log(1./x),h);
  }
  *upv=f[0];
  *dnv=f[1];
  *glu=f[2];
  dusum=f[3];
  dudif=f[4];
  *str=f[5];
  *chm=f[6];
  *bot=f[7];
  *dsea=0.5*(dusum+dudif);
  *usea=0.5*(dusum-dudif);
  if (q2 < 2.7) *chm=0.;  /* Can't make charm and bottom if Q**2 < M**2 */
  if (q2 < 30.) *bot=0.;
  if (*chm < 0.) *chm=0.; /* Evidently the fits can go negative */
  if (*bot < 0.) *bot=0.;
}