SUBROUTINE polint(xa,ya,n,x,y,dy) 
IMPLICIT NONE 
INTEGER(4),PARAMETER :: NMAX=10
INTEGER(4) :: n,i,m,ns
REAL(8) :: dy,x,y,xa(n),ya(n) 
REAL(8) :: den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)

!Largest anticipated value of n. 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(xai) = ya_i, i = 1, . . . , n, then the returned value y = P(x). 
 
ns=1 
dif=abs(x-xa(1))

do i=1,n 
   dift=abs(x-xa(i)) 
   if (dift<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==0.)pause 'failure in polint' 
      den=w/den 
      d(i)=hp*den 
      c(i)=ho*den 
   enddo
   if (2*ns<n-m)then 
      dy=c(ns+1) 
   else 
      dy=d(ns) 
      ns=ns-1 
   endif 
   y=y+dy 
enddo

return 
END