      subroutine polint(xa,ya,x,y,dy)
      use types
      IMPLICIT NONE

      real(SP), dimension(:), intent(in) :: xa,ya
      real(SP), intent(in)::x
      real(SP), intent(out) ::y,dy
      integer(I4B):: n, NMAX, i, m,ns
      integer(I4B), dimension(1)::tns
      parameter(NMAX=10)
      real(SP), dimension(size(xa))::c,d,den,ho
      real(SP)::dif,dift,hp,w



      if(size(xa,dim=1).ne.size(ya,dim=1)) then
      write(6,*) 'polint: xa and ya have different sizes'
      stop
      endif
            
      n=size(xa,dim=1)
      c=ya
      d=ya
      ho=xa-x
      tns=minloc(abs(x-xa))
      ns=tns(1)
      y=ya(ns)
      ns=ns-1
	do m=1,n-1
		den(1:n-m)=ho(1:n-m)-ho(1+m:n)
		if (any(den(1:n-m) == 0.0)) then
		write(6,*)	'polint: calculation failure'
                stop
                endif
		den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
		d(1:n-m)=ho(1+m:n)*den(1:n-m)
		c(1:n-m)=ho(1:n-m)*den(1:n-m)
		if (2*ns < n-m) then
			dy=c(ns+1)
		else
			dy=d(ns)
			ns=ns-1
		end if
		y=y+dy
	end do
	END SUBROUTINE polint
      
      




