subroutine fitpdf(fscale,pdfset,Np)
	implicit double precision (a-h,o-z)
        COMMON /fordyi/ tau
 	INTEGER porder,Np,iset
	double precision tau,fscale
        character*60 pdfset,prefix
	dimension cc(5)

	dimension xx(NP) 
	dimension yglu(Np), yupv(Np), ydnv(Np)
	dimension yusea(Np), ydsea(Np), ystr(Np)
	dimension ychm(Np), ybot(Np)
	dimension xsel(3), ysel(3)
	external SetCtq6,Ctq6Pdf

CCCC    USE THIS ONLY FOR CTEQ PDFS =================================

        if (pdfset.eq.'CTEQ66') then        
           Iset = 400
           call SetCtq6(Iset)
        endif 

CCCC ================================================================


	do i=1,Np
	   xx(i) = tau + (1-tau)*i*(i-1d0)/Np/(Np-1)
	   xi1=xx(i)

	   q=fscale

CCCC     Call PDF set ===============================================

           if (pdfset.eq.'MRSTLO') then
              mode =1       
              call mrstlo(xi1,q,mode,upv,dnv,usea,dsea,str,chm,
     &          bot,glu)  
           elseif (pdfset.eq.'MRST2001') then
              mode =1       
              call mrst2001(xi1,q,mode,upv,dnv,usea,dsea,str,chm,
     &          bot,glu)
           elseif (pdfset.eq.'MRSTNNLO') then
              mode=1
              call mrstnnlo(xi1,q,mode,upv,dnv,usea,dsea,str,chm,
     &          bot,glu)
           elseif (pdfset.eq.'MRST2006NNLO') then 
              prefix = "mrst2006/Grids/mrst2006nnlo" ! prefix for the grid files
              iset = 0
              call GetAllPDFs6(prefix,iset,xi1,q,
     &           upv,dnv,usea,dsea,str,sbar,chm,bot,glu)
           elseif (pdfset.eq.'MSTW2008LO') then
              prefix = "mstw2008/Grids/mstw2008lo"
              iset=0
              call GetAllPDFs(prefix,iset,xi1,q,
     &           upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot)
           elseif (pdfset.eq.'MSTW2008NLO') then  
              prefix = "mstw2008/Grids/mstw2008nlo"      
              iset=0
              call GetAllPDFs(prefix,iset,xi1,q,
     &           upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot)
           elseif (pdfset.eq.'MSTW2008NNLO') then  
              prefix = "mstw2008/Grids/mstw2008nnlo"      
              iset=0
              call GetAllPDFs(prefix,iset,xi1,q,
     &           upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot)
           elseif (pdfset.eq.'CTEQ66') then
              upv = Ctq6Pdf(1,xi1,q)*xi1
              dnv = Ctq6Pdf(2,xi1,q)*xi1
              usea = Ctq6Pdf(-1,xi1,q)*xi1
              dsea = Ctq6Pdf(-2,xi1,q)*xi1
              str = Ctq6Pdf(3,xi1,q)*xi1
              chm = Ctq6Pdf(4,xi1,q)*xi1
              bot = Ctq6Pdf(5,xi1,q)*xi1
              glu = Ctq6Pdf(0,xi1,q)*xi1

              upv = upv-usea
              dnv = dnv-dsea
           endif

CCC====================================================================

	   yupv(i) = upv
	   ydnv(i) = dnv
	   yusea(i) = usea
	   ydsea(i) = dsea
	   ystr(i) = str
	   ychm(i) = chm
	   ybot(i) = bot
	   yglu(i) = glu
	enddo	

	hh = Np*(Np-1)/(1-tau)
	
        open(21,FILE='GL.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='UPV.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='DNV.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='USEA.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='DSEA.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='STR.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
	open(21,FILE='CHM.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)
        open(21,FILE='BTM.dat',STATUS='NEW')
	write(21,*) hh	
	close(21)

	do i=1,Np 
	   if (i==1) then
	      i1=i
	      i2=i+1
	      i3=i+3
	   elseif (i==Np) then 
	      i1=i-2
	      i2=i-1
	      i3=i
	   else
	      i1=i-1
	      i2=i
	      i3=i+1
	   endif 
	   
	   xa1 = xx(i1)
	   xa2 = xx(i2)
	   xa3 = xx(i3)
	   
	   	   
	   ya1 =yglu(i1)
	   ya2 =yglu(i2)
	   ya3 =yglu(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='GL.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)

	   	   
	   ya1 =yupv(i1)
	   ya2 =yupv(i2)
	   ya3 =yupv(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='UPV.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   
	   ya1 =ydnv(i1)
	   ya2 =ydnv(i2)
	   ya3 =ydnv(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='DNV.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   ya1 =yusea(i1)
	   ya2 =yusea(i2)
	   ya3 =yusea(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='USEA.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   
	   ya1 =ydsea(i1)
	   ya2 =ydsea(i2)
	   ya3 =ydsea(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='DSEA.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   ya1 =ystr(i1)
	   ya2 =ystr(i2)
	   ya3 =ystr(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='STR.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   
	   ya1 =ychm(i1)
	   ya2 =ychm(i2)
	   ya3 =ychm(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='CHM.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   ya1 =ybot(i1)
	   ya2 =ybot(i2)
	   ya3 =ybot(i3)
	   
	   call givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	   
	   open(21,FILE='BTM.dat',STATUS='OLD',ACCESS='append')		
	   write(21,*) xx(i), a, b, c
	   close(21)
	   
	   
	enddo
	
	end


	subroutine givecoeff(xa1,xa2,xa3,ya1,ya2,ya3, a, b, c)
	implicit double precision (a-h,o-z)
	
	den=(xa1-xa2)*(xa2-xa3)*(xa3-xa1)
	a = -ya1*(xa2-xa3)/den
	a = a-ya2*(xa3-xa1)/den
	a = a-ya3*(xa1-xa2)/den
      
	b =  ya1*(xa2*xa2-xa3*xa3)/den
	b = b+ya2*(xa3*xa3-xa1*xa1)/den
	b = b+ya3*(xa1*xa1-xa2*xa2)/den
      
	c = -ya1*(xa2-xa3)*xa2*xa3/den
	c = c-ya2*(xa3-xa1)*xa1*xa3/den
	c = c-ya3*(xa1-xa2)*xa1*xa2/den
	return
	end 




c       gluon pdf
	
	double precision function retpdfg(x)
	implicit double precision(a-h,o-z)    
	COMMON /fordyi/ tau
	double precision tau,lambda
	dimension cc(5)
	COMMON /gparam/ xxg, AAG, BBG, CCG, hhg
	dimension xxg(1500), AAG(1500), BBG(1500), CCG(1500)
	
	hh=hhg
	
	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))
	
	L =rr 
	
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfg = CCG(L)+BBG(L)*x+AAG(L)*x**2
	else
	   retpdfg = 0.0d0
	endif

	end


c     up-valence pdf

      double precision function retpdfuv(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /uvparam/ xxuv, AAUV, BBUV, CCUV, hhuv
      dimension xxuv(1500), AAUV(1500), BBUV(1500), CCUV(1500)

	hh=hhuv

	lambda = (x-tau)*hh

	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfuv = CCUV(L)+BBUV(L)*x+AAUV(L)*x**2
	else
	   retpdfuv = 0.0d0
	endif

	end



c     down-valence pdf

      double precision function retpdfdv(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /dvparam/ xxdv, AADV, BBDV, CCDV, hhdv
      dimension xxdv(1500), AADV(1500), BBDV(1500), CCDV(1500)

	hh=hhdv

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfdv = CCDV(L)+BBDV(L)*x+AADV(L)*x**2
	else
	   retpdfdv = 0.0d0
	endif

	end



c     up-sea pdf

      double precision function retpdfus(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /usparam/ xxus, AAUS, BBUS, CCUS, hhus
      dimension xxus(1500), AAUS(1500), BBUS(1500), CCUS(1500)

	hh=hhus

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfus = CCUS(L)+BBUS(L)*x+AAUS(L)*x**2
	else
	   retpdfus = 0.0d0
	endif

	end



c     down-sea pdf

      double precision function retpdfds(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /dsparam/ xxds, AADS, BBDS, CCDS, hhds
      dimension xxds(1500), AADS(1500), BBDS(1500), CCDS(1500)

	hh=hhds

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfds = CCDS(L)+BBDS(L)*x+AADS(L)*x**2
	else
	   retpdfds = 0.0d0
	endif

	end


c     strange-sea pdf

      double precision function retpdfstr(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /strparam/ xxstr, AASTR, BBSTR, CCSTR, hhstr
      dimension xxstr(1500), AASTR(1500), BBSTR(1500), CCSTR(1500)

	hh=hhstr

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfstr = CCSTR(L)+BBSTR(L)*x+AASTR(L)*x**2
	else
	   retpdfstr = 0.0d0
	endif

	end


c     charm-sea pdf

      double precision function retpdfchm(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /chmparam/ xxchm, AACHM, BBCHM, CCCHM, hhchm
      dimension xxchm(1500), AACHM(1500), BBCHM(1500), CCCHM(1500)

	hh=hhchm

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfchm = CCCHM(L)+BBCHM(L)*x+AACHM(L)*x**2
	else
	   retpdfchm = 0.0d0
	endif

	end



c     bottom-sea pdf

      double precision function retpdfbtm(x)
      implicit double precision(a-h,o-z)    
      COMMON /fordyi/ tau
      double precision tau,lambda
      dimension cc(5)
      COMMON /btmparam/ xxbtm, AABTM, BBBTM, CCBTM, hhbtm
      dimension xxbtm(1500), AABTM(1500), BBBTM(1500), CCBTM(1500)

	hh=hhbtm

	lambda = (x-tau)*hh
	rr = 0.5d0*(1.0d0+dsqrt(dabs(1+4d0*lambda)))

	L =rr 
	if ((L.gt.0).and.(L.lt.1501)) then
	   retpdfbtm = CCBTM(L)+BBBTM(L)*x+AABTM(L)*x**2
	else
	   retpdfbtm = 0.0d0
	endif

	end