C     ----------------------------------------------------------------
      program myevol
C     ----------------------------------------------------------------
C--   Custom evolution unpolarised LO 
C--   Check result aigainst build-in unpolarised evolution 
C     ----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      
      data iord/1/, nfin/0/, iosp/3/               !LO, VFNS, iospl
      data as0/0.364/, r20/2.D0/                   !alphas 
      data q2c/7.D0/, q2b/25.D0/, q2t/1000.D0/     !thresholds
      data q20/5.D0/                               !q20

      dimension xmin(5),ixwt(5)
      data xmin/1.D-5,0.2D0,0.4D0,0.6D0,0.75D0/    !multiple x grid
      data ixwt/    1,    2,    4,    8,    16/    !xgrid weights
      data nxin/100/                               !x grid pnts

      dimension qq(2),qw(2)
      data qq/2.D0,1.D4/, qw/1.D0,1.D0/, nqin/60/  !mu2 grid

      external func                                !input parton dists
      dimension def(-6:6,12)                       !flavor composition
      data def  /
C--   tb  bb  cb  sb  ub  db   g   d   u   s   c   b   t
C--   -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6  
     + 0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,   !dval
     + 0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,   !uval
     + 0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,   !sval
     + 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,   !dbar
     + 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,   !ubar
     + 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !sbar
     + 0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,   !cval
     + 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,   !cbar
     + 52*0.    /
     
      dimension pdfs1(-6:6), pdfs2(-6:6)
      
      real time1, time2
      
      external myweight

C--   lun = 6 stdout, -6 stdout w/out banner page
      lun    = 6
      lunout = abs(lun) 
      call qcinit(lun,' ')
      call gxmake(xmin,ixwt,5,nxin,nx,iosp)
      call gqmake(qq,qw,2,nqin,nq)

C--   Weight tables
      call readwt(20,'../weights/unpol.wgt',id1,id2,nw,ierr)
      if(ierr.ne.0) then
        call fillwt(1,id1,id2,nw)
        call dmpwgt(1,20,'../weights/unpol.wgt')
      endif
      call readwt(20,'../weights/custom.wgt',id1,id2,nw,ierr)
      if(ierr.ne.0) then
        call fillwc(myweight,id1,id2,nw)
        call dmpwgt(4,20,'../weights/custom.wgt')
      endif 

C--   Evolution      
      call setord(iord)                         !LO, NLO, NNLO
      call setalf(as0,r20)                      !input alphas
      iqc  = iqfrmq(q2c)                        !charm threshold
      iqb  = iqfrmq(q2b)                        !bottom threshold
      iqt  = iqfrmq(q2t)                        !top threshold
      call setcbt(nfin,iqc,iqb,iqt)             !thesholds in the vfns
      iq0  = iqfrmq(q20)                        !starting scale
      call cpu_time(time1)                      !start timing
      call evolfg(1,func,def,iq0,eps)           !unpolarised evolution
      call evolfg(4,func,def,iq0,eps)           !custom evolution 
      
C--   Compare the two evolutions
C--   Pdfs vs xi at q2min and q2max
      xmil = log(xmin(1))
      xmal = 0.D0
      nxxx = 10
      bw   = (xmal-xmil)/nxxx
      dmax = 0.D0
      do ix = 1,nxxx
        xl = xmil + (ix-0.5)*bw
        xx = exp(xl)
        do iq = 1,2
          call fpdfxq(1,xx,qq(iq),pdfs1,1)      !unpolarised
          call fpdfxq(4,xx,qq(iq),pdfs2,1)      !custom evolution
          do id = -6,6
            dif  = pdfs1(id)-pdfs2(id)
            dmax = max(dmax,abs(dif))
          enddo
        enddo  
      enddo
C--   End timing      
      call cpu_time(time2)

C--   This should be zero    
      write(lunout,'(/'' This should be zero dmax = '', E15.5)') dmax 
      
      write(lunout,'(/'' time spent:'',F7.4)') time2-time1         

      end
    
C     ================================================================
      
C     =================================================
      subroutine myweight(w,nw,nwords,idpij,mxord,idum)
C     =================================================

C--   Fill LO Pij

C--   w           (in)   store
C--   nw          (in)   number of words available
C--   nwords      (out)  number of words used < 0 not enough space
C--   idpij       (out)  list of Pij table identifiers
C--   mxord       (out)  maximum perturbative order LO,NLO,NNLO
C--   idum               not used at present        

      implicit double precision (a-h,o-z)
      
      dimension w(*)
      dimension idpij(7,3)
      dimension itypes(4)
      
      external AChi
      external PQQR, PQQS, PQQD                         ! PQQ
      external PQGA                                     ! PQG
      external PGQA                                     ! PGQ
      external PGGA, PGGR, PGGS, PGGD                   ! PGG

C--   S/R name for error message      
      call setUmsg('MyWeight')
C--   Max perturbative order
      mxord = 1
C--   Partition
      itypes(1) = 2
      itypes(2) = 2 
      itypes(3) = 0
      itypes(4) = 0
      call BookTab(w,nw,itypes,nwords)
C--   Not enough space       
      if(nwords.le.0) return
C--   Assign table indices
      idPij(1,1)   =  101                               ! PQQ
      idPij(2,1)   =  201                               ! PQG
      idPij(3,1)   =  102                               ! PGQ
      idPij(4,1)   =  202                               ! PGG
      idPij(5,1)   =  101                               ! PPL
      idPij(6,1)   =  101                               ! PMI
      idPij(7,1)   =  101                               ! PVA
C--   Fill tables      
      call MakeWRS(w, idPij(1,1), PQQR, PQQS, AChi, 0)
      call MakeWtD(w, idPij(1,1), PQQD, AChi)
      call MakeWtA(w, idPij(2,1), PQGA, AChi)
      call MakeWtA(w, idPij(3,1), PGQA, AChi)
      call MakeWtA(w, idPij(4,1), PGGA, AChi)
      call MakeWRS(w, idPij(4,1), PGGR, PGGS, AChi, 0)
      call MakeWtD(w, idPij(4,1), PGGD, AChi)
C--   Done!            
      idum  = 0
      call clrUmsg

      return
      end
      
C==   ================================================================      

C     ====================================
      double precision function Achi(qmu2)
C     ====================================

      implicit double precision (a-h,o-z)
      
      qq    = qmu2    !avoid compiler warning
      Achi  = 1.D0
      
      return
      end
      
C     =========================================
      DOUBLE PRECISION FUNCTION PQQR(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning
      IDUM = NF        !avoid compiler warning 
      PQQR = 4.D0 * ( 1.D0 + X*X ) / 3.D0
 
      RETURN
      END

C     =========================================
      DOUBLE PRECISION FUNCTION PQQS(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning
      IDUM = NF        !avoid compiler warning 
      PQQS = 1.D0 / (1.D0-X)
 
      RETURN
      END 
      
C     =============================================
      double precision function PQQD(x,qmu2,nf)
C     =============================================

      implicit double precision (a-h,o-z)
      
      DDUM = X       !avoid compiler warning
      DDUM = QMU2    !avoid compiler warning
      IDUM = NF      !avoid compiler warning
      PQQD = 2.D0
      
      return
      end 
      
C     =========================================
      DOUBLE PRECISION FUNCTION PQGA(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning
      PQGA = NF * ( X*X + (1.D0-X)*(1.D0-X) )
 
      RETURN
      END

C     =========================================
      DOUBLE PRECISION FUNCTION PGQA(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning
      IDUM = NF        !avoid compiler warning 
      PGQA = 4.D0 * ( 1.D0 + (1.D0-X)*(1.D0-X) ) / ( 3.D0*X )
 
      RETURN
      END

C     =========================================
      DOUBLE PRECISION FUNCTION PGGA(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM  = QMU2      !avoid compiler warning
      IDUM  = NF        !avoid compiler warning       
      ONEMX = 1.D0 - X
      PGGA  = 6.D0 * ( ONEMX/X + X*ONEMX )
 
      RETURN
      END
      
C     =========================================
      DOUBLE PRECISION FUNCTION PGGR(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning
      IDUM = NF        !avoid compiler warning       
      PGGR = 6.D0 * X 
 
      RETURN
      END
      
C     =========================================
      DOUBLE PRECISION FUNCTION PGGS(X,QMU2,NF)
C     =========================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DDUM = QMU2      !avoid compiler warning 
      IDUM = NF        !avoid compiler warning      
      PGGS = 1.D0 / ( 1.D0 - X )
 
      RETURN
      END
      
C     =========================================
      double precision function PGGD(X,QMU2,NF)
C     =========================================

      implicit double precision (a-h,o-z)
      
      XX   = X       !avoid compiler warning
      QQ   = QMU2    !avoid compiler warning
      PGGD = 6.D0*(11.D0/12.D0 - NF/18.D0)
      
      return
      end                                        
                                      
C==   ================================================================

C     ======================================
      double precision function func(ipdf,x)
C     ======================================

      implicit double precision (a-h,o-z)

                     func = 0.D0
      if(ipdf.eq. 0) func = xglu(x)
      if(ipdf.eq. 1) func = xdnv(x)
      if(ipdf.eq. 2) func = xupv(x)
      if(ipdf.eq. 3) func = 0.D0
      if(ipdf.eq. 4) func = xdbar(x)
      if(ipdf.eq. 5) func = xubar(x)
      if(ipdf.eq. 6) func = xsbar(x)
      if(ipdf.eq. 7) func = 0.D0
      if(ipdf.eq. 8) func = 0.D0
      if(ipdf.eq. 9) func = 0.D0
      if(ipdf.eq.10) func = 0.D0
      if(ipdf.eq.11) func = 0.D0
      if(ipdf.eq.12) func = 0.D0

      return
      end
 
C     =================================
      double precision function xupv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data au /5.107200D0/

      xupv = au * x**0.8D0 * (1.D0-x)**3.D0

      return
      end

C     =================================
      double precision function xdnv(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ad /3.064320D0/

      xdnv = ad * x**0.8D0 * (1.D0-x)**4.D0

      return
      end

C     =================================
      double precision function xglu(x)
C     =================================

      implicit double precision (a-h,o-z)

      data ag    /1.7D0 /

      common /msum/ glsum, uvsum, dvsum

      xglu = ag * x**(-0.1D0) * (1.D0-x)**5.D0

      return
      end

C     ==================================
      double precision function xdbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      data adbar /0.1939875D0/

      xdbar = adbar * x**(-0.1D0) * (1.D0-x)**6.D0

      return
      end

C     ==================================
      double precision function xubar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xubar = xdbar(x) * (1.D0-x)

      return
      end

C     ==================================
      double precision function xsbar(x)
C     ==================================

      implicit double precision (a-h,o-z)

      xsbar = 0.2D0 * (xdbar(x)+xubar(x))

      return
      end     







