subroutine polint(xa,ya,n,x,y,dy)

	implicit none

	integer::n
	real::dy,x,y,xa(n),ya(n)
	integer,parameter::NMAX=10

!   Given arrays xa and ya, each of length n, and given a value x, 
!   this routine returns a value y, and an error estimate dy. If 
!   P(x) is the polynomial of degree N-1 such that P(xa(i))=ya(i)
!   ,i=1..n, the the returned value y=P(x).

        integer::i,m,ns
	real::den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)

	ns=1
	dif=abs(x-xa(1))

	do i=1,n
	  dift=abs(x-xa(i))
	  if(dift.lt.dif)then
	    ns=i
	    dif=dift
          endif
	  c(i)=ya(i)
	  d(i)=ya(i)
        enddo

	y=ya(ns)
	ns=ns-1
	do m=1,n-1
	  do i=1,n-m
	    ho=xa(i)-x
	    hp=xa(i+m)-x
	    w=c(i+1)-d(i)
	    den=ho-hp
	    if(den.eq.0.) stop
	    den=w/den
	    d(i)=hp*den
	    c(i)=ho*den
          enddo
          if(2*ns.lt.n-m)then
	    dy=c(ns+1)
          else
	    dy=d(ns)
	    ns=ns-1
          endif
	  y=y+dy
        enddo

	return
	end subroutine polint