*******************************************************************
*	18/12/08									*
*	Interpolation subroutine in Fortran 77 for calculation 	*
*	of skewed distributions	from using the Shuvaev transform.  	*
*											*	
*	See arXiv:0812.3558 [hep-ph] for details				*
*	Version 1.0									*
*	Written by CJN								*
*											*
*	pdf=1 is for CTEQ6.6M partons requires grid cteq6.6m.dat	*
*	pdf=2 is for CTEQ6L partons requires grid cteq6l.dat		*
*	pdf=3 is for MRST2004NLO requires grid mrstnlo.dat		*
*	pdf=4 is for MRST2004LO requires grid mrstlo.dat		*
*	pdf=5 is for MSTW2008NLO partons requires grid mstwnlo.dat	*
*	pdf=6 is for MSTW2008LO partons requires grid mstwlo.dat	*
*******************************************************************

      SUBROUTINE intpol(x,xi,qsq,pdf,resq,resqb,resg)
C Gives the result of the shuvaev transform on quarks (resq), antiquarks (resqb) and gluons (resg)
C for diagonal partons at a given x, xi, qsq (GeV^2).
C The result for the quarks and antiquarks is (x+xi)*F(x,xi)
      IMPLICIT NONE
      INTEGER nq,i,j,k,fl,init,xlpt,rlpt,qlpt,nx,nr,pdf,pt
	PARAMETER (nx=85,nq=20,nr=43)
      REAL*8 x,qsq,q,xmin,xmax,ximin,ximax,qsqmin,qsqmax,xd,xi
      REAL*8 xx(nx),qqsq(nq),rt(nr),xxl(nx),qsql(nq),rtl(nr)
      REAL*8 resq,resqb,resg,r,ipol(3),skd(3,nx,nr,nq)
      REAL*8 xb(2),rb(2),qb(2),f(3,2,2,2)
      DATA xmin,xmax,ximin,ximax,qsqmin,qsqmax/4.d-5,1.d0,0.d0,1.d0,
     &1.25d0,80.d0/
      COMMON/gval/xx,qqsq,rt
      DATA init/0/
      SAVE skd,xxl,rtl,qsql
      IF (x.LT.xmin.OR.x.GT.xmax) then
		print 12, x
      	stop
      ELSEIF (qsq.LT.qsqmin.OR.qsq.GT.qsqmax) then
		print 13, qsq
      	stop
      ELSEIF (x.LT.xi) then
		print 14
      	stop
      ENDIF
	IF (x.GT.0.1d0) print 15,x
      IF (init.EQ.0) then
      IF (pdf.EQ.1) OPEN(10,file='cteq6.6m.dat',status='old')
	IF (pdf.EQ.2) OPEN(10,file='cteq6l.dat',status='old')
      IF (pdf.EQ.3) OPEN(10,file='mrstnlo.dat',status='old')
	IF (pdf.EQ.4) OPEN(10,file='mrstlo.dat',status='old')
      IF (pdf.EQ.5) OPEN(10,file='mstwnlo.dat',status='old')
	IF (pdf.EQ.6) OPEN(10,file='mstwlo.dat',status='old')
	READ (10,*)
	READ (10,*)
	READ (10,*)
      DO 20 i=1,nx
      DO 20 j=1,nr
      DO 20 k=1,nq
      	READ (10,11) skd(1,i,j,k),skd(2,i,j,k),skd(3,i,j,k)
20    CONTINUE
      CLOSE(10)
      DO i=1,nx
	   xxl(i)=DLOG(xx(i))
	ENDDO
	DO i=1,nr-1
	   rtl(i)=-DLOG(rt(i))
	ENDDO
      DO i=1,nq
	   qsql(i)=DLOG(DSQRT(qqsq(i)))
     	ENDDO
      ENDIF
      init=1
      xd=DLOG(x)
      IF (xi/x.GT.rt(nr-1)) r=DLOG(xi/x)
      IF (xi/x.LE.rt(nr-1)) r=xi/x
      q=DLOG(DSQRT(qsq))
      xlpt=pt(xxl,nx,xd)
	IF (xi/x.GT.rt(nr-1)) rlpt=pt(rtl,nr-1,-r)
      IF (xi/x.LE.rt(nr-1)) rlpt=nr-1
      qlpt=pt(qsql,nq,q)
      DO i=1,2
      	  xb(i)=xxl(xlpt+i-1)
      	  qb(i)=qsql(qlpt+i-1)
      ENDDO
      IF (xi/x.GT.rt(nr-1)) then
		rb(1)=-rtl(rlpt)
		rb(2)=-rtl(rlpt+1)
	ELSEIF (xi/x.LE.rt(nr-1)) then
		rb(1)=rt(rlpt)
		rb(2)=rt(rlpt+1)
	ENDIF
      DO 40 fl=1,3
      DO 30 i=1,2
      DO 30 j=1,2 
      DO 30 k=1,2 
      f(fl,i,j,k)=skd(fl,xlpt+i-1,rlpt+j-1,qlpt+k-1)
30    CONTINUE      
      ipol(fl)=(f(fl,1,1,1)*(xb(2)-xd)*(rb(2)-r)*(qb(2)-q)+
     &	   f(fl,2,1,1)*(xd-xb(1))*(rb(2)-r)*(qb(2)-q)+
     &         f(fl,1,2,1)*(xb(2)-xd)*(r-rb(1))*(qb(2)-q)+
     &	   f(fl,2,2,1)*(xd-xb(1))*(r-rb(1))*(qb(2)-q)+
     &         f(fl,1,1,2)*(xb(2)-xd)*(rb(2)-r)*(q-qb(1))+
     &         f(fl,2,1,2)*(xd-xb(1))*(rb(2)-r)*(q-qb(1))+
     &         f(fl,1,2,2)*(xb(2)-xd)*(r-rb(1))*(q-qb(1))+
     &	   f(fl,2,2,2)*(xd-xb(1))*(r-rb(1))*(q-qb(1)))/
     &         ((xb(2)-xb(1))*(rb(2)-rb(1))*(qb(2)-qb(1)))
40    CONTINUE
      resq=ipol(1)
      resqb=ipol(2)
      resg=ipol(3)     
11    FORMAT(4E15.6,I3.1)
12    FORMAT('  WARNING:  x value is out of range,  ','x= ',e10.5)
13    FORMAT('  WARNING:  q^2 value is out of range,  ','q2= ',e10.5)
14	FORMAT('  WARNING:  x<xi')
15	FORMAT('  WARNING:  Routine called with large x,  ','x= ',e10.5)
      RETURN
      END

	BLOCK DATA
	IMPLICIT NONE
	INTEGER nx,nq,nr
	PARAMETER (nx=85,nq=20,nr=43)
	REAL*8 xx(nx),qqsq(nq),rt(nr)
      COMMON/gval/xx,qqsq,rt
	DATA xx/0.4000E-04,0.4500E-04,0.5000E-04,0.5500E-04,
     &0.6000E-04,0.6500E-04,0.7000E-04,0.7500E-04,0.8000E-04,
     &0.8500E-04,0.9000E-04,0.9500E-04,0.1000E-03,0.1500E-03,
     &0.2000E-03,0.2500E-03,0.3000E-03,0.3500E-03,0.4000E-03,
     &0.4500E-03,0.5000E-03,0.5500E-03,0.6000E-03,0.6500E-03,
     &0.7000E-03,0.7500E-03,0.8000E-03,0.8500E-03,0.9000E-03,
     &0.9500E-03,0.1000E-02,0.1500E-02,0.2000E-02,0.2500E-02,
     &0.3000E-02,0.3500E-02,0.4000E-02,0.4500E-02,0.5000E-02,
     &0.5500E-02,0.6000E-02,0.6500E-02,0.7000E-02,0.7500E-02,
     &0.8000E-02,0.8500E-02,0.9000E-02,0.9500E-02,0.1000E-01,
     &0.1500E-01,0.2000E-01,0.2500E-01,0.3000E-01,0.3500E-01,
     &0.4000E-01,0.4500E-01,0.5000E-01,0.5500E-01,0.6000E-01,
     &0.6500E-01,0.7000E-01,0.7500E-01,0.8000E-01,0.8500E-01,
     &0.9000E-01,0.9500E-01,0.1000E+00,0.1500E+00,0.2000E+00,
     &0.2500E+00,0.3000E+00,0.3500E+00,0.4000E+00,0.4500E+00,
     &0.5000E+00,0.5500E+00,0.6000E+00,0.6500E+00,0.7000E+00,
     &0.7500E+00,0.8000E+00,0.8500E+00,0.9000E+00,0.9500E+00,
     &0.1000E+01/
	DATA rt/0.1000E+01,0.9950E+00,0.9900E+00,0.9850E+00,
     &0.9800E+00,0.9650E+00,0.9500E+00,0.9250E+00,0.9000E+00,
     &0.8750E+00,0.8500E+00,0.8250E+00,0.8000E+00,0.7750E+00,
     &0.7500E+00,0.7250E+00,0.7000E+00,0.6500E+00,0.6000E+00,
     &0.5500E+00,0.5000E+00,0.4500E+00,0.4000E+00,0.3500E+00,
     &0.3000E+00,0.2500E+00,0.2000E+00,0.1500E+00,0.1000E+00,
     &0.7500E-01,0.5000E-01,0.3000E-01,0.1000E-01,0.7500E-02,
     &0.5000E-02,0.3000E-02,0.1000E-02,0.5500E-03,0.1000E-03,
     &0.5500E-04,0.1000E-04,0.5000E-05,0.0000E+00/
      DATA qqsq/1.25d0,1.5d0,2.d0,2.5d0,3.2d0,4.d0,5.d0,6.4d0,8.d0,
     &        10.d0,12.d0,18.d0,20.d0,25.d0,30.d0,40.d0,50.d0,60.d0,
     &	  70.d0,80.d0/
      SAVE
	END

      INTEGER FUNCTION pt(xx,nx,x)
	INTEGER nx,upp,low,mid
      REAL*8 x,xx(nx)
      IF (x.EQ.xx(1)) then
      	pt=1
		GOTO 10
      ELSEIF (x.EQ.xx(nx)) then
      	pt=nx-1  
		GOTO 10
	ENDIF
      upp=nx+1
      low=0
11    IF ((upp-low).LE.1) GOTO 12
      mid=(upp+low)/2
      IF (x.GE.xx(mid)) then
      	low=mid
      ELSE
      	upp=mid
      ENDIF
      GOTO 11
12    pt=low
10    RETURN
      END