C--   G.W. 02/12/2008 Example program for standalone structure function
C--   code corresponding to exactly that used in the MSTW 2008 fits.
C--   Compile with:
C--    g77 example_sf.f mstw2008structurefunctions.f mstwpdf.f alphaS.f

      PROGRAM structurefunctions
      IMPLICIT NONE
      INTEGER iord,ipn,ieigen,neigen,ix,iq
      PARAMETER(neigen=20)
      COMMON/iordCommon/iord
      DOUBLE PRECISION ALPHAS,x,q,GetOnePDF,xg,
     &     f2,f2c,f2b,fl,flc,flb,
     &     f2_p,f2c_p,f2b_p,fl_p,flc_p,flb_p,
     &     f2_m,f2c_m,f2b_m,fl_m,flc_m,flb_m,
     &     f2_max,f2c_max,f2b_max,fl_max,flc_max,flb_max,
     &     f2_min,f2c_min,f2b_min,fl_min,flc_min,flb_min
      INTEGER alphaSorder,alphaSnfmax
      DOUBLE PRECISION distance,tolerance,
     &     mCharm,mBottom,alphaSQ0,alphaSMZ
      COMMON/mstwCommon/distance,tolerance,
     &     mCharm,mBottom,alphaSQ0,alphaSMZ,alphaSorder,alphaSnfmax
      INTEGER ISET
      CHARACTER prefix*50,cl*4,pdfset*50
      COMMON/mstwfiles/iset,prefix,cl

      CALL WATE96               ! needed for Gaussian integration

      IORD = 2                  ! 0 = LO, 1 = NLO, 2 = NNLO
      IPN = 1                   ! 1 = proton, 2 = neutron

      IF (IORD.EQ.0) THEN
         pdfset = "mstw2008lo"
      ELSE IF (IORD.EQ.1) THEN
         pdfset = "mstw2008nlo"
      ELSE IF (IORD.EQ.2) THEN
         pdfset = "mstw2008nnlo"
      END IF
      prefix = "Grids/"//pdfset
      cl = '68cl'
C      cl = '90cl'
      iset = 0                  ! central fit
C--   Heavy quark masses and alphaS values
C--   are stored in a COMMON block.
      xg = GetOnePDF(prefix,iset,0.1d0,10d0,0) ! dummy call
      WRITE(6,*) "mCharm = ",mCharm,", mBottom = ",mBottom
      WRITE(6,*) "alphaS(Q0) = ",alphaSQ0,", alphaS(MZ) = ",
     &     alphaSMZ,", alphaSorder = ",alphaSorder,
     &     ", alphaSnfmax = ",alphaSnfmax

C--   Call the initialisation routine with alpha_S(Q_0).
      CALL INITALPHAS(alphaSorder,1.D0,1.D0,alphaSQ0,
     &     mCharm,mBottom,1.D10)
C--   Check calculated value of alpha_S(M_Z) matches stored value.
      WRITE(6,'(" alphaS(MZ) = ",F7.5," = ",F7.5)')
     &     ALPHAS(91.1876D0),alphaSMZ

C--   Calculate F_2 vs. x in bins of Q2.
C      IORD = 2 ! Temporary: change order of coefficient functions, not PDFs.
      DO iq = 1, 4
         IF (iq.EQ.1) q = sqrt(2d0)
         IF (iq.EQ.2) q = sqrt(5d0)
         IF (iq.EQ.3) q = sqrt(20d0)
         IF (iq.EQ.4) q = sqrt(100d0)
         IF (iq.EQ.1) OPEN(UNIT=66,FILE="f2_"//
     &        pdfset(1:len_trim(pdfset))//"_2GeV2.dat")
         IF (iq.EQ.2) OPEN(UNIT=66,FILE="f2_"//
     &        pdfset(1:len_trim(pdfset))//"_5GeV2.dat")
         IF (iq.EQ.3) OPEN(UNIT=66,FILE="f2_"//
     &        pdfset(1:len_trim(pdfset))//"_20GeV2.dat")
         IF (iq.EQ.4) OPEN(UNIT=66,FILE="f2_"//
     &        pdfset(1:len_trim(pdfset))//"_100GeV2.dat")
         DO ix = 1, 50
            x = 10.d0**(log10(1.D-5) + (ix-1.D0)/(50.D0-1.D0)*
     &           (log10(0.9D0)-log10(1.D-5)))

C--   Calculate F2,F2c,F2b,FL,FLc,FLb with central PDFs.
         iset = 0
         CALL MSTWNC(x,q,ipn,f2,f2c,f2b,fl,flc,flb)
         f2_max=0d0; f2c_max=0d0; f2b_max=0d0
         fl_max=0d0; flc_max=0d0; flb_max=0d0
         f2_min=0d0; f2c_min=0d0; f2b_min=0d0
         fl_min=0d0; flc_min=0d0; flb_min=0d0
         DO ieigen = 1, neigen
            iset = 2*ieigen-1   ! copy to COMMON block
            CALL MSTWNC(x,q,ipn,f2_p,f2c_p,f2b_p,fl_p,flc_p,flb_p)
            iset = 2*ieigen     ! copy to COMMON block
            CALL MSTWNC(x,q,ipn,f2_m,f2c_m,f2b_m,fl_m,flc_m,flb_m)
C--   Add contribution to maximum F2,F2c,F2b,FL,FLc,FLb.
            f2_max = f2_max + (max(f2_p-f2,f2_m-f2,0d0))**2
            f2c_max = f2c_max + (max(f2c_p-f2c,f2c_m-f2c,0d0))**2
            f2b_max = f2b_max + (max(f2b_p-f2b,f2b_m-f2b,0d0))**2
            fl_max = fl_max + (max(fl_p-fl,fl_m-fl,0d0))**2
            flc_max = flc_max + (max(flc_p-flc,flc_m-flc,0d0))**2
            flb_max = flb_max + (max(flb_p-flb,flb_m-flb,0d0))**2
C--   Add contribution to minimum F2,F2c,F2b,FL,FLc,FLb.
            f2_min = f2_min + (min(f2_p-f2,f2_m-f2,0d0))**2
            f2c_min = f2c_min + (min(f2c_p-f2c,f2c_m-f2c,0d0))**2
            f2b_min = f2b_min + (min(f2b_p-f2b,f2b_m-f2b,0d0))**2
            fl_min = fl_min + (min(fl_p-fl,fl_m-fl,0d0))**2
            flc_min = flc_min + (min(flc_p-flc,flc_m-flc,0d0))**2
            flb_min = flb_min + (min(flb_p-flb,flb_m-flb,0d0))**2
         END DO
         WRITE(6,*) x,q**2,f2,sqrt(f2_max),sqrt(f2_min)
         WRITE(66,*) x,q**2,f2,sqrt(f2_max),sqrt(f2_min)
      END DO ! ix
      CLOSE(UNIT=66)
      END DO ! iq

      STOP
      END