***********************************************************************
*                                                                     *
*    Program for generating Electromagnetic F2 Structure Functions    *
*    using parton distributions in the DIS scheme according to MRS99  *
*    The F2's generated here correspond to the F2's generated by MRST *
*    partons in the MS bar scheme.         Therefore included is a    *
*    consistent treatment of charm and bottom structure functions     *
*    Not included are the effects due to NLO corrections to photon-   *
*    gluon fusion.   Charm mass = 1.43 GeV   Bottom mass = 4.3 GeV    *
*    A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne         *
*    Durham preprint DTP/99/64      hep-ph/9907231                    *
*    Rutherford Appleton Laboratory preprint RAL-TR-1999-047          *
*    Setting MODE=1 gives MRST      L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=2 gives MRST(g^^) L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=3 gives MRST(gvv) L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=4 gives MRST(avv) L(4)=229 MeV a_s(M_Z)=0.1125      *
*    Setting MODE=5 gives MRST(a^^) L(4)=383 MeV a_s(M_Z)=0.1225      *
*    Setting MODE=6 gives MRST(q^)  L(4)=303 MeV a_s(M_Z)=0.1178      *
*    Setting MODE=7 gives MRST(qv)  L(4)=293 MeV a_s(M_Z)=0.1171      *
*    Setting MODE=8 gives MRST(s^)  L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=9 gives MRST(sv)  L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=10 gives MRST(cv) L(4)=300 MeV a_s(M_Z)=0.1175      *
*    Setting MODE=11 gives MRST(c^) L(4)=300 MeV a_s(M_Z)=0.1175      *
*                                                                     *
*    The program should be run only with iord set to 1                *
*                                                                     *
***********************************************************************

      implicit real*8(a-h,o-z)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord

      call wate96
      iord=1
      flavor=3
      qsdt=8.18
      qsct=74.
   99 read(5,*) mode,x,q2

      if(mode.eq.1) alambda=0.300 
      if(mode.eq.2) alambda=0.300 
      if(mode.eq.3) alambda=0.300 
      if(mode.eq.4) alambda=0.229 
      if(mode.eq.5) alambda=0.383
      if(mode.eq.6) alambda=0.303 
      if(mode.eq.7) alambda=0.293 
      if(mode.eq.8) alambda=0.300 
      if(mode.eq.9) alambda=0.300 
      if(mode.eq.10) alambda=0.300
      if(mode.eq.11) alambda=0.300 

      call disfun(x,q2,mode,f2p,f2n,f2c,f2b)

      print 101,x,q2,f2p,f2n,f2c,f2b

  101 format(1x,f10.6,1x,f12.5,7(1x,f10.4))
      go to 99
      end

      subroutine disfun(x,q2,mode,f2p,f2n,f2c,f2b)
      implicit real*8(a-h,o-z)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord
      data pi,pi2/3.14159,9.8696/
      XLAM=alambda
      xlam2=xlam*xlam
      t=dlog(q2/xlam2)
      al=alpha(t)/(4.*pi)
      scale=dsqrt(q2)
      cf=4./3.
      ca=3.

      epsc4=qsdt/q2
      epsc=epsc4/4.
      epsb4=qsct/q2
      epsb=epsb4/4.
      diffc=(1.-epsc)
      if(diffc.lt.0.) diffc=0.
      smochm=diffc**0.5
      diffb=(1.-epsb)
      if(diffb.lt.0.) diffb=0.
      smobot=diffb**0.5

      call mrs99dis(x,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      fp=(4.*upv+dnv+8.*usea+2.*dsea+2.*str)/9.
      fn=(4.*dnv+upv+2.*usea+8.*dsea+2.*str)/9.
      fc=8.*chm/9.
      fb=2.*bot/9.

      print 97,glu,chm
   97 format(1x,'glu=',f10.6,1x,'chm=',f10.6)

      ffp=fp
      ffn=fn
      ffc=0.
      ffb=0.

      facc=4./9.
      facb=1./9.

      AL1=dLOG(1.-X)

      DO 23 I=1,NTERMS
      Y=0.5*(1.-X)*XI(I)+0.5*(1.+X)
      XY=X/Y
      AL1=dLOG(1.-Y)
      call mrs99dis(xy,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      gluxy=glu
      CG2=(-1.+8.*Y*(1.-Y)+(1.-2.*Y+2.*Y*Y)*dLOG(1./Y-1.))
      if(epsc.lt.1.) then
      ffc=ffc - smochm*0.5*(1.-X)*WI(I)*AL*2.*facc*CG2*gluxy
      endif
      if(epsb.lt.1.) then
      ffb=ffb - smobot*0.5*(1.-X)*WI(I)*AL*2.*facb*CG2*gluxy
      endif
   23 CONTINUE
   21 CONTINUE

      print 98,ffc
   98 format(1x,'98 ',f10.6)

      xcmax=1./(1.+epsc4)
      if(xcmax.le.x) go to 321
      DO 323 I=1,NTERMS
      Y=0.5*(xcmax-X)*XI(I)+0.5*(xcmax+X)
      XY=X/Y
      call mrs99dis(xy,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsc.gt.1.) chm=0.
      gluxy=glu
      fcxy=8.*chm/9.
      del=0.01
      xyu=xy*(1.+0.5*del)
      if(xyu.gt.1.d0) then
      dfcxy=0.d0
      go to 1324
      endif
      call mrs99dis(xyu,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsc.gt.1.) chm=0.
      fcxyu=8.*chm/9.
      xyl=xy*(1.-0.5*del)
      call mrs99dis(xyl,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsc.gt.1.) chm=0.
      fcxyl=8.*chm/9.
      dfcxy=(fcxyu-fcxyl)/del
 1324 c0c=cheavy(2,y,epsc)
      if(epsc.gt.1.d0) c0c=0.d0
      cg21c=2.*facc*cheavy(1,y,epsc)
      cg22c=2.*facc*c0c*dlog(1./epsc)

      ffc=ffc+0.5*(xcmax-x)*wi(i)*c0c*(-dfcxy+3.*fcxy)
      ffc=ffc+0.5*(xcmax-x)*wi(i)*al*(cg21c-cg22c)*gluxy

      fcextra=coeff2(y,epsc)
      ffc=ffc+0.5*(xcmax-x)*wi(i)*fcxy*fcextra

  323 CONTINUE
  321 CONTINUE

      xbmax=1./(1.+epsb4)
      if(xbmax.le.x) go to 421
      DO 423 I=1,NTERMS
      Y=0.5*(xbmax-X)*XI(I)+0.5*(xbmax+X)
      XY=X/Y
      call mrs99dis(xy,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsb.gt.1.) bot=0.
      gluxy=glu
      fbxy=2.*bot/9.
      del=0.01
      xyu=xy*(1.+0.5*del)
      if(xyu.gt.1.d0) then
      dfbxy=0.d0
      go to 1424
      endif
      call mrs99dis(xyu,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsb.gt.1.) bot=0.
      fbxyu=2.*bot/9.
      xyl=xy*(1.-0.5*del)
      call mrs99dis(xyl,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      if(epsb.gt.1.) bot=0.
      fbxyl=2.*bot/9.
      dfbxy=(fbxyu-fbxyl)/del

 1424 c0b=cheavy(2,y,epsb)
      if(epsb.gt.1.d0) c0b=0.d0
      cg21b=2.*facb*cheavy(1,y,epsb)
      cg22b=2.*facb*c0b*dlog(1./epsb)

      ffb=ffb+0.5*(xbmax-x)*wi(i)*c0b*(-dfbxy+3.*fbxy)
      ffb=ffb+0.5*(xbmax-x)*wi(i)*al*(cg21b-cg22b)*gluxy

      fbextra=coeff2(y,epsb)
      ffb=ffb+0.5*(xbmax-x)*wi(i)*fbxy*fbextra

  423 CONTINUE
  421 CONTINUE

      if(ffc.lt.0.) ffc=0.
      if(ffb.lt.0.) ffb=0.
      f2p=ffp+ffc+ffb
      f2n=ffn+ffc+ffb
      f2c=ffc
      f2b=ffb

   27 RETURN
      END

      real*8 function coeff2(y,eps)
      implicit real*8(a-h,o-z)
      dimension z1(17),z2(17)
  
      data z1/-0.183E+01,0.400E+01,0.159E+02,-0.357E+02,
     .-0.186E+00,-0.988E+02,0.712E+02,0.631E+02,-0.136E+01,
     .0.175E+03,-0.158E+03,-0.433E+02,0.375E+01,-0.913E+02,
     .0.842E+02,0.107E+02,0.629E+00/
      data z2/-0.204E+01,-0.127E+01,0.117E+01,-0.208E+00,
     .-0.228E+02,0.261E+03,-0.340E+03,0.153E+03,-0.591E+02,
     .-0.199E+04,-0.113E+04,0.299E+04,-0.408E+04,-0.446E+04,
     . 0.307E+05,-0.282E+05,0.144E+02/

      epsl=eps
      eps4=4.d0*eps
      x0=1.d0/(1.d0+eps4)
      yx0=y/x0
      yx01=1.d0-yx0
      epsl2=epsl*epsl
      epsl3=epsl2*epsl
      if(y.lt.0.05) go to 100    
      a0=z1(1)+z1(2)*epsl+z1(3)*epsl2+z1(4)*epsl3
      a1=z1(5)+z1(6)*epsl+z1(7)*epsl2+z1(8)*epsl3
      a2=z1(9)+z1(10)*epsl+z1(11)*epsl2+z1(12)*epsl3
      a3=z1(13)+z1(14)*epsl+z1(15)*epsl2+z1(16)*epsl3
      fac=a0+a1*yx0+a2*yx0*yx0+a3*yx0*yx0*yx0
      coeff2=fac*yx01**z1(17)
      return
  100 continue
      a0=z2(1)+z2(2)*epsl+z2(3)*epsl2+z2(4)*epsl3
      a1=z2(5)+z2(6)*epsl+z2(7)*epsl2+z2(8)*epsl3
      a2=z2(9)+z2(10)*epsl+z2(11)*epsl2+z2(12)*epsl3
      a3=z2(13)+z2(14)*epsl+z2(15)*epsl2+z2(16)*epsl3
      fac=a0+a1*yx0+a2*yx0*yx0+a3*yx0*yx0*yx0
      coeff2=fac*yx01**z2(17)
      return
      end

      real*8 function coeff3(y,eps)
      implicit real*8(a-h,o-z)
      dimension z1(17),z2(17)

      data z1/-0.857E+00,0.217E+02,-0.199E+02,0.162E+02,
     .0.204E+00,-0.142E+03,0.101E+03,-0.575E+02,0.832E+00,
     .0.175E+03,-0.139E+03,0.604E+02,0.344E-01,-0.648E+02,
     .0.538E+02,-0.160E+02,0.000E+00/
      data z2/-0.334E-03,0.677E-02,-0.108E-02,0.825E-03,
     .-0.425E+04,0.883E+05,-0.562E+05,0.601E+05,0.562E+07,
     .-0.122E+09,0.912E+08,-0.928E+08,-0.236E+10,0.499E+11,
     .-0.449E+11,0.412E+11,0.000E+00/

      epsl=eps
      eps4=4.d0*eps
      x0=1.d0/(1.d0+eps)
      yx0=y/x0
      yx01=1.d0-yx0
      epsl2=epsl*epsl
      epsl3=epsl2*epsl
      if(y.lt.0.05) go to 100    
      a0=z1(1)+z1(2)*epsl+z1(3)*epsl2+z1(4)*epsl3
      a1=z1(5)+z1(6)*epsl+z1(7)*epsl2+z1(8)*epsl3
      a2=z1(9)+z1(10)*epsl+z1(11)*epsl2+z1(12)*epsl3
      a3=z1(13)+z1(14)*epsl+z1(15)*epsl2+z1(16)*epsl3
      fac=a0+a1*yx0+a2*yx0*yx0+a3*yx0*yx0*yx0
      coeff3=fac*yx01**z1(17)
      return
  100 continue
      a0=z2(1)+z2(2)*epsl+z2(3)*epsl2+z2(4)*epsl3
      a1=z2(5)+z2(6)*epsl+z2(7)*epsl2+z2(8)*epsl3
      a2=z2(9)+z2(10)*epsl+z2(11)*epsl2+z2(12)*epsl3
      a3=z2(13)+z2(14)*epsl+z2(15)*epsl2+z2(16)*epsl3
      fac=a0+a1*yx0+a2*yx0*yx0+a3*yx0*yx0*yx0
      coeff3=fac*yx01**z2(17)
      return
      end

      DOUBLE PRECISION FUNCTION	ALPHA(T)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord
      DATA PI/3.14159/
      DATA TOL/.0005/
      ITH=0
      TT=T
      qsdtt=qsdt/4.
      qsctt=qsct/4.
      AL=ALAMBDA
      AL2=AL*AL
      FLAV=4.
      QS=AL2*dEXP(T)

      if(qs.lt.0.5d0) then   !!  running stops below 0.5
          qs=0.5d0
          t=dlog(qs/al2)
          tt=t
      endif

      IF(QS.gt.QSCTT) GO	TO 12  
      IF(QS.lt.QSDTT) GO	TO 312  
   11 CONTINUE
      B0=11-2.*FLAV/3. 
      IF(IORD)1,1,2
c     IF(IORD)2,2,2	!TAKE CARE !!
    1 CONTINUE
      ALPHA=4.*PI/B0/T
      RETURN
    2 CONTINUE
      X1=4.*PI/B0
      B1=102.-38.*FLAV/3.
      X2=B1/B0**2
      AS=X1/T*(1.-X2*dLOG(T)/T)
    5 CONTINUE
      F=-T+X1/AS-X2*dLOG(X1/AS+X2)
      FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
      AS2=AS-F/FP
      DEL=ABS(F/FP/AS)
      IF(DEL-TOL)3,3,4
    3 CONTINUE
      ALPHA=AS2
      IF(ITH.EQ.0) RETURN
      GO TO (13,14,15) ITH
    4 CONTINUE
      AS=AS2
      GO TO 5
   12 ITH=1
      T=dLOG(QSCTT/AL2)
      GO TO 11
   13 ALFQC4=ALPHA
      FLAV=5.   
      ITH=2
      GO TO 11
   14 ALFQC5=ALPHA
      ITH=3
      T=TT
      GO TO 11
   15 ALFQS5=ALPHA
      ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
      ALPHA=1./ALFINV
      RETURN

  311 CONTINUE
      B0=11-2.*FLAV/3. 
      IF(IORD)31,31,32
c     IF(IORD)32,32,32	!TAKE CARE !!
   31 CONTINUE
      ALPHA=4.*PI/B0/T
      RETURN
   32 CONTINUE
      X1=4.*PI/B0
      B1=102.-38.*FLAV/3.
      X2=B1/B0**2
      AS=X1/T*(1.-X2*dLOG(T)/T)
   35 CONTINUE
      F=-T+X1/AS-X2*dLOG(X1/AS+X2)
      FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
      AS2=AS-F/FP
      DEL=ABS(F/FP/AS)
      IF(DEL-TOL)33,33,34
   33 CONTINUE
      ALPHA=AS2
      IF(ITH.EQ.0) RETURN
      GO TO (313,314,315) ITH
   34 CONTINUE
      AS=AS2
      GO TO 35
  312 ITH=1
      T=dLOG(QSDTT/AL2)
      GO TO 311
  313 ALFQC4=ALPHA
      FLAV=3.   
      ITH=2
      GO TO 311
  314 ALFQC3=ALPHA
      ITH=3
      T=TT
      GO TO 311
  315 ALFQS3=ALPHA
      ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
      ALPHA=1./ALFINV
      RETURN
      END


      SUBROUTINE WATE96
C*******************************************************************
C*****							       *****
C***** THE X(I)	AND W(I) ARE THE DIRECT	OUTPUT FROM A PROGRAM  *****
C***** USING NAG ROUTINE D01BCF	TO CALCULATE THE	       *****
C***** GAUSS-LEGENDRE WEIGHTS FOR 96 POINT INTEGRATION.	       *****
C***** THEY AGREE TO TYPICALLY 14 DECIMAL PLACES WITH THE      *****
C***** TABLE IN	ABRAMOWITZ & STEGUN, PAGE 919.		       *****
C*****							       *****
C***** ---->   PETER HARRIMAN, APRIL 3RD 1990.		       *****
C*****							       *****
C*******************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION	X(48),W(48)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      NTERMS=96

      X( 1)=   0.01627674484960183561
      X( 2)=   0.04881298513604856015
      X( 3)=   0.08129749546442434360
      X( 4)=   0.11369585011066471632
      X( 5)=   0.14597371465489567682
      X( 6)=   0.17809688236761733390
      X( 7)=   0.21003131046056591064
      X( 8)=   0.24174315616383866556
      X( 9)=   0.27319881259104774468
      X(10)=   0.30436494435449495954
      X(11)=   0.33520852289262397655
      X(12)=   0.36569686147231213885
      X(13)=   0.39579764982890709712
      X(14)=   0.42547898840729897474
      X(15)=   0.45470942216774136446
      X(16)=   0.48345797392059470382
      X(17)=   0.51169417715466604391
      X(18)=   0.53938810832435567233
      X(19)=   0.56651041856139533470
      X(20)=   0.59303236477757022282
      X(21)=   0.61892584012546672523
      X(22)=   0.64416340378496526886
      X(23)=   0.66871831004391424358
      X(24)=   0.69256453664216964528
      X(25)=   0.71567681234896561582
      X(26)=   0.73803064374439816819
      X(27)=   0.75960234117664555964
      X(28)=   0.78036904386743123629
      X(29)=   0.80030874413913884180
      X(30)=   0.81940031073792957139
      X(31)=   0.83762351122818502758
      X(32)=   0.85495903343459936363
      X(33)=   0.87138850590929436968
      X(34)=   0.88689451740241818933
      X(35)=   0.90146063531585023110
      X(36)=   0.91507142312089592706
      X(37)=   0.92771245672230655266
      X(38)=   0.93937033975275308073
      X(39)=   0.95003271778443564022
      X(40)=   0.95968829144874048809
      X(41)=   0.96832682846326217918
      X(42)=   0.97593917458513455843
      X(43)=   0.98251726356301274934
      X(44)=   0.98805412632962202890
      X(45)=   0.99254390032376081654
      X(46)=   0.99598184298720747465
      X(47)=   0.99836437586317963722
      X(48)=   0.99968950388322870559
      W( 1)=   0.03255061449236316962
      W( 2)=   0.03251611871386883307
      W( 3)=   0.03244716371406427668
      W( 4)=   0.03234382256857594104
      W( 5)=   0.03220620479403026124
      W( 6)=   0.03203445623199267876
      W( 7)=   0.03182875889441101874
      W( 8)=   0.03158933077072719007
      W( 9)=   0.03131642559686137819
      W(10)=   0.03101033258631386231
      W(11)=   0.03067137612366917839
      W(12)=   0.03029991542082762553
      W(13)=   0.02989634413632842385
      W(14)=   0.02946108995816795100
      W(15)=   0.02899461415055528410
      W(16)=   0.02849741106508543861
      W(17)=   0.02797000761684838950
      W(18)=   0.02741296272602931385
      W(19)=   0.02682686672559184485
      W(20)=   0.02621234073567250055
      W(21)=   0.02557003600534944960
      W(22)=   0.02490063322248370695
      W(23)=   0.02420484179236479915
      W(24)=   0.02348339908592633665
      W(25)=   0.02273706965832950717
      W(26)=   0.02196664443874448477
      W(27)=   0.02117293989219144572
      W(28)=   0.02035679715433347898
      W(29)=   0.01951908114014518992
      W(30)=   0.01866067962741165898
      W(31)=   0.01778250231604547316
      W(32)=   0.01688547986424539715
      W(33)=   0.01597056290256253144
      W(34)=   0.01503872102699521608
      W(35)=   0.01409094177231515264
      W(36)=   0.01312822956696188190
      W(37)=   0.01215160467108866759
      W(38)=   0.01116210209983888144
      W(39)=   0.01016077053500880978
      W(40)=   0.00914867123078384552
      W(41)=   0.00812687692569928101
      W(42)=   0.00709647079115442616
      W(43)=   0.00605854550423662775
      W(44)=   0.00501420274292825661
      W(45)=   0.00396455433844564804
      W(46)=   0.00291073181793626202
      W(47)=   0.00185396078894924657
      W(48)=   0.00079679206555731759
      DO 1 I=1,48
      XI(I)=-X(49-I)
      WI(I)=W(49-I)
      XI(I+48)=X(I)
      WI(I+48)=W(I)
    1 CONTINUE
      DO 2 I=1,96
    2 XX(I)=0.5*(XI(I)+1.)
      XX(97)=1.0
      EXPON=1.0
      DO 3 I=1,96
      YI=2.*(0.5*(1.+XI(I)))**EXPON-1.
      WI(I)=WI(I)/(1.+XI(I))*(1.+YI)*EXPON
      XI(I)=YI
      XX(I)=0.5*(1.+YI)
    3 CONTINUE
      RETURN
      END

      real*8 function cheavy(i,z,eps)
      implicit real*8(a-h,o-z)

c     this function returns the values of C_g(z,Q^2) 
c     and the deriv. wrt log Q^2. Here eps=m^2/Q^2.
c     If i=1  C_g for F2.  If i=2 deriv of C_g for F2 
c     If i=3  C_g for FL.  If i=4 (1-m2/Q2)*beta*Clq(massless)

      if(i.gt.4) stop
      z1=1.-z
      z2=z*z
      eps2=eps*eps
      beta2=1.-4.*eps*z/z1
      if(beta2.lt.0.) go to 10
      beta=dsqrt(beta2)
      a=z2+z1*z1
      b=4.*z*(1.-3.*z)
      c=-8.*z2
      aa=8.*z*z1-1.
      bb=-4.*z*z1
      arg=(1.+beta)/(1.-beta)
      fac=dlog(arg)
      cf=4./3.
      go to (1,2,3,4) i
    1 cheavy=(a+b*eps+c*eps2)*fac+(aa+bb*eps)*beta
      return
    2 cheavy=(-b*eps-2.*c*eps2)*fac+(a+b*eps+c*eps2)/beta
     .      +(-bb*eps)*beta +(aa+bb*eps)*2.*z*eps/z1/beta
      return
    3 cheavy=-bb*beta+c*eps*fac
      return
    4 cheavy=(1.-eps)*beta*4.*cf*z
      return
   10 print 99
   99 format(1x,'x > x0')
      stop
      end