c
c     This subroutine is used to perform an interpolation/extrapolation
c     of various components after an inversion process has been performed.
c     The subroutine performes the interpolation/extrapolation based on
c     a polynomial of degree N (where N represents the # of points passed
c     to this subroutine).  The code was taken from "NUMERICAL RECIPES 
c     The Art of Scientific Computing" Cambridge University Press 1986.
c
c     This routine is generic and can be used for other applications
c     other than the inversion routine. It needs to be put in the math
c     library
c
c     Program output is y.  Note you can also output error dy, it is calculated
c     here but not passed anywhere.  
c
c     A.N.Hawkins 7/21/95 
c
      subroutine polint
     &                 (xa,ya,n,x,y)
	  
	  
	  real x, y, dif, dift,ho,hp,dy,den,w  
	  integer i ,m,n ,ns,nmax
      parameter (nmax=10)
	  real xa(n), ya(n), c(nmax), d(nmax)
	  ns=1 
	  dif=ABS(x-xa(1))
	  do 11 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)
11    continue
	  y=ya(ns)
	  ns=ns-1

	  do 13 m=1,n-1
	    do 12 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.)pause
		  den=w/den
		  d(i)=hp*den
		  c(i)=ho*den
12      continue
      if (2*ns.lt.n-m)then
	    dy=c(ns+1) 
	  else
	    dy=d(ns) 
		ns=ns-1
	  endif
	  y=y+dy
13    continue
      return
	  end