      subroutine matrixelem(vars,wlo,Nmax,xorki)
      implicit none
      INTEGER wlo,jj,isetval, Nmax
      double precision vars(2)
      character*60 collider, pdfset
      double precision muf, xorki(Nmax)

      double precision sect1,sect2,sect3,sect4,sect5,sect6,sect7
      double precision x1, x2, tau, x, z, ris, tmp0, tmp1, tmp2, tmp

      double precision L, Lt, Lr, alphas, NF, meas 
      double precision Rplus, Rdelta, Rreg
      
      double precision Lgg,Lqg,Lqqb,Lqq,Lq1q2
      double precision Lggo,Lqgo,Lqqbo,Lqqo,Lq1q2o

      double precision ggdelta, ggplus, ggreg, qgreg, qqbreg, qqreg
      double precision q1q2reg
      common /xsection/ L, Lt, Lr, alphas, NF 
      COMMON /kinem/ tau

      COMMON /lumiB/ collider, pdfset, muf

      external lumif



      Rplus=1d0
      Rdelta =1d0
      Rreg=1d0

C---- Evaluate  functions in the two kinematic points ------ 
      x=vars(1)
      z = 1d0
      meas =1/x
      x1 = x
      x2 = tau/x1/z

      sect1 = ggdelta(z, wlo)



      x=vars(1)
      z = vars(2)
      meas =1/x
      x1 = x
      x2 = tau/x1/z

      sect2 = ggplus(z,wlo)


      sect3 = qgreg(z,wlo)
      sect4 = qqbreg(z,wlo)
      sect5 = qqreg(z,wlo)
      sect6 = q1q2reg(z,wlo)
      sect7 = ggreg(z,wlo)



      do jj=1,Nmax

         ris=0d0
         isetval=jj-1
         
C-------PARTONIC LUMINOSITIES -----------------
         x=vars(1)
         z = 1d0
         meas =1/x
         x1 = x
         x2 = tau/x1/z
         
         call lumif(collider,x1,x2,muf,pdfset, isetval, Lggo,Lqgo,
     &        Lqqbo,Lqqo,Lq1q2o)
         
         
         x=vars(1)
         z = vars(2)
         meas =1/x
         x1 = x
         x2 = tau/x1/z
         
         
         call lumif(collider,x1,x2,muf,pdfset, isetval, Lgg,Lqg,
     &        Lqqb,Lqq,Lq1q2)
         
         
         
C-------------------------------------------------------------
         

C     -------- DELTA FUNCTION TERMS -------------
         
         x=vars(1)
         z = 1d0
         
         meas =1/x
         x1 = x
         x2 = tau/x1/z
         
         tmp = Lggo*meas*sect1
         tmp = Rdelta*tmp
         ris = ris +tmp
         
         
C     -------- PLUS FUNCTION TERMS --------------
         
         tmp1 = Lggo*meas
         
         x=vars(1)
         z = vars(2)
         
         meas =1/x
         x1 = x
         x2 = tau/x1/z
         tmp2 = Lgg*meas
         
         tmp = (tmp2-tmp1)/(1d0-z)
         
         tmp = tmp*sect2
         
         tmp =Rplus*tmp
         ris = ris + tmp
         
         
C     ----------- REGULAR TERMS -----------------
         
         
         x=vars(1)
         z = vars(2)
         meas =1/x
         x1 = x
         x2 = tau/x1/z

         tmp=0d0
         tmp = tmp+Lqg*sect3
         tmp = tmp+Lqqb*sect4
         tmp = tmp+Lqq*sect5
         tmp = tmp+Lq1q2*sect6
         tmp = tmp+Lgg*sect7

         tmp = tmp*meas
         tmp =Rreg*tmp
         ris =ris+tmp

         xorki(jj)=ris
      enddo
      end


      double precision function ggreg(x,wlo)
      implicit none 
      double precision x, xggrN, xggrNN
      double precision tmp
      integer wlo
      if (wlo.EQ.0) then 
         tmp=0d0
      elseif (wlo.EQ.1) then 
         tmp = 1d0*xggrN(x)
      elseif (wlo.EQ.2) then 
         tmp = 1d0*xggrNN(x)
      endif
         
      ggreg = tmp     
      end






      double precision function qgreg(x,wlo)
      implicit none 
      double precision x, xqgN, xqgNN
      double precision tmp
      integer wlo
      
      if (wlo.EQ.0) then 
         tmp=0d0
      elseif (wlo.EQ.1) then 
         tmp = 1d0*xqgN(x)
      elseif (wlo.EQ.2) then 
         tmp = 1d0*xqgNN(x)
      endif
         
      qgreg = tmp     
      end




      double precision function qqbreg(x,wlo)
      implicit none 
      double precision x, xqqbN, xqqbNN
      double precision tmp
      integer wlo
      
      if (wlo.EQ.0) then 
         tmp=0d0
      elseif (wlo.EQ.1) then 
         tmp = 1d0*xqqbN(x)
      elseif (wlo.EQ.2) then 
         tmp = 1d0*xqqbNN(x)
      endif
         qqbreg = tmp     
      end




      double precision function qqreg(x,wlo)
      implicit none 
      double precision x, xqqNN
      double precision tmp
      integer wlo
      
      if (wlo.EQ.0) then 
         tmp=0d0
      elseif (wlo.EQ.1) then 
         tmp = 0d0
      elseif (wlo.EQ.2) then 
         tmp = xqqNN(x)*1d0
      endif
         qqreg = tmp     
      end




      double precision function q1q2reg(x,wlo)
      implicit none 
      double precision x, xqqprime
      double precision tmp
      integer wlo

      
      if (wlo.EQ.0) then 
         tmp=0d0
      elseif (wlo.EQ.1) then 
         tmp = 0d0
      elseif (wlo.EQ.2) then 
         tmp = xqqprime(x)*1d0
      endif
         q1q2reg = tmp     
      end



      double precision function ggdelta(z,wlo)
      implicit none 
      double precision tmp
      double precision z, tmp1, tmp2, tmp3, z4
      double precision L, Lt, Lr, alphas, NF, pi, z2, z3
      integer wlo


      common /xsection/ L, Lt, Lr, alphas, NF 
      common /consts/ pi, z2, z3


      z4=1.082323234d0
      tmp1 = 1d0 
      tmp2 = (-23.D0/6.D0*L+23.D0/6.D0*Lr+6*z2+11.D0/2.D0)
      tmp3 =
     # ((529.D0/48.D0-18*z2)*L**2+(-23*z2-317.D0/12.D0-171.D0/2.D0*z
     #3)*L-529.D0/24.D0*Lr*L+137.D0/24.D0*Lt+349.D0/6.D0*z2-445.D0/12.D0
     #*z3+303.D0/8.D0-54*z2**2+1071.D0/8.D0*z4+529.D0/48.D0*Lr**2+(69.D0
     #/2.D0*z2+875.D0/24.D0)*Lr)
 

      if (wlo.EQ.0) then 
         tmp=tmp1
      elseif (wlo.EQ.1) then 
         tmp = tmp1*1d0+alphas*tmp2*1d0
      elseif (wlo.EQ.2) then 
         tmp = tmp1*1d0+alphas*tmp2*1d0+alphas**2*tmp3*1d0
      endif
      
      ggdelta = tmp

      end


      double precision function ggplus(z,wlo)
      implicit none 
      double precision z,x, LLz, DD0, DD1, DD2, DD3, taux
      double precision L, Lt, Lr, alphas, NF, tmp
      double precision tmp1, tmp2, tmp3, z4, pi, z2, z3

      integer wlo

      common /xsection/ L, Lt, Lr, alphas, NF 

      common /consts/ pi, z2, z3       
      x=z

      z4=1.082323234d0


      LLz = dlog(1d0-z)

      DD0=1d0
      DD1=LLz
      DD2=LLz**2
      DD3=LLz**3

      tmp1=0d0      
      tmp2 = -6*L*DD0+12*DD1

      taux =115.D0/4.D0*L**2+(-349.D0/6.D0+45*z2)*L-69.D0/2.D0*Lr*L
     &     +23*z2-233.D0/9.D0+351.D0/2.D0*z3
      tmp3 =taux*DD0 
     &     +(349.D0/3.D0-90*z2+69*Lr-46*L+36*L**2)*DD1
     &     +(-23-108*L)*DD2 
     &     +72*DD3


      if (wlo.EQ.0) then 
         tmp=tmp1
      elseif (wlo.EQ.1) then 
         tmp = tmp1+alphas*tmp2*1d0
      elseif (wlo.EQ.2) then 
         tmp = tmp1+alphas*tmp2*1d0+alphas**2*tmp3*1d0
      endif

      
      ggplus = tmp
      end
