************************************************************************
*
*     StructureFunctionxQ.f:
*
*     This function returns the value of the StructureFunctions at the
*     scale Q in GeV and for the bjorken variable x using the interpolation
*     of the cached StructureFunctions.
*
************************************************************************
      function StructureFunctionxQ(proc,sf,comp,x,Q)
*
      implicit none
*
      include "../commons/grid.h"
      include "../commons/gridQ.h"
      include "../commons/scales.h"
      include "../commons/StructureFunctionsxQ.h"
      include "../commons/TimeLike.h"
**
*     Input Variables
*
      double precision x,Q
      character*2 proc
      character*2 sf
      character*6 comp
**
*     Internal Variables
*
      integer iproc,isf,icomp
      integer idx,idQ
      integer alpha,tau
      integer lx,ux,lQ,uQ
      integer isg,iQtrans,iQtransp
      integer tQ,ideg
      double precision Q2
      double precision w_int_xQ,wg
      double precision tol
      parameter(tol=1d-10)
**
*     Output Variables
*
      double precision StructureFunctionxQ
*
*     Check that the APFEL evolution has been initialized
*
      if(InCacheSFs.ne."done")then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Impossible to use this function because the"
         write(6,*) "structure functions have not been cached."
         write(6,*) "Call 'CacheStructureFunctionsAPFEL(Q0)'"
         write(6,*) "before calling StructureFunctionxQ."
         write(6,*) "   "
         call exit(-10)
      endif
*
      Q2 = Q * Q
*
*     Check consistency of the input variables
*
      if(proc.ne."EM".and.proc.ne."NC".and.proc.ne."CC")then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Invalid process, proc =",proc
         call exit(-10)
      endif
*
      if(proc.eq."EM") iproc = 0
      if(proc.eq."NC") iproc = 1
      if(proc.eq."CC") iproc = 2
*
      if(sf.ne."F2".and.sf.ne."FL".and.sf.ne."F3")then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Invalid StructureFunction, sf =",sf
         call exit(-10)
      endif
*
      if(sf.eq."F2") isf = 1
      if(sf.eq."FL") isf = 2
      if(sf.eq."F3") isf = 3
*
      if(comp(1:5).ne."light".and.comp(1:5).ne."charm".and.
     1   comp(1:6).ne."bottom".and.comp(1:3).ne."top".and.
     2   comp(1:5).ne."total")then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Invalid component, comp =",comp
         call exit(-10)
      endif
*
      if(comp(1:5).eq."light")  icomp = 3
      if(comp(1:5).eq."charm")  icomp = 4
      if(comp(1:6).eq."bottom") icomp = 5
      if(comp(1:3).eq."top")    icomp = 6
      if(comp(1:5).eq."total")  icomp = 7
*
*     Check that structure functions have been called in the allowed region
*
      if(x.lt.xmin(1)*(1d0-tol).or.x.gt.xmax*(1d0+tol))then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Value of x out of range, x =",x
         call exit(-10)
      endif
*
      if(Q2.lt.Q2min*(1d0-tol).or.Q2.gt.Q2max*(1d0+tol))then
         write(6,*) "In StructureFunctionxQ.f:"
         write(6,*) "Value of Q out of range, Q =",Q
         call exit(-10)
      endif
*
*     Trim away numerical fluctiations
*
      if(x.lt.xmin(1)) x  = xmin(1)
      if(x.gt.xmax)    x  = 1d0
*
      if(Q2.lt.Q2min)  Q2 = Q2min
      if(Q2.gt.Q2max)  Q2 = Q2max
*
*     1) If the value of Q2 happens to fall in the tiny gap
*     between two subsequent subgrids, push it on the first
*     node of the second grid.
*     2) Check that the Q2-subgrid has enough points for the
*     interpolation, otherwise reduce the interpolation degree.
*     3) Determine tQ, i.e. the parameter to avoid to interpolate
*     over the thresholds.
*
      idQ      = inter_degreeQ
      tQ       = 0
      iQtransp = 0
      iQtrans  = nQ(nfin) - 1
      do isg=nfin,nffi
*     1)
         if(Q2.gt.Q2g(iQtrans).and.Q2.lt.Q2g(iQtrans+1))
     1        Q2 = Q2g(iQtrans+1)
*     2)
         if(Q2.ge.Q2g(iQtransp).and.Q2.le.Q2g(iQtrans).and.
     1      nQ(isg).lt.idQ+1) idQ = nQ(isg) - 1
*     3)
         do ideg=2,idQ
            if(Q2.gt.Q2g(iQtrans-ideg+1).and.Q2.le.Q2g(iQtrans))
     1           tQ = ideg - 1
         enddo
         iQtransp = iQtrans + 1
         iQtrans  = iQtrans + nQ(isg+1)
      enddo
*
*     Interpolation
*
      idx = inter_degree(0)
      do alpha=0,nin(0)
         if(xg(0,alpha).gt.x) goto 11
      enddo
 11   lx = alpha - 1
      ux = lx + idx + 1
      do tau=0,nQ2g
         if(Q2g(tau).gt.Q2) goto 12
      enddo
 12   lQ = tau - 1 - tQ
      uQ = lQ + idQ + 1
      StructureFunctionxQ = 0d0
      do alpha=lx,ux
         do tau=lQ,uQ
            wg = w_int_xQ(tQ,idx,idQ,alpha,tau,x,Q2)
            StructureFunctionxQ = StructureFunctionxQ
     1           + wg * SFxQ(iproc,isf,icomp,alpha,tau)
         enddo
      enddo
*
      if(Timelike) StructureFunctionxQ = StructureFunctionxQ / x
*
      return
      end