function fintr(x,f,xf,npf,nintr)
c==========================================================
c fintr.f: polinomial interpolation
c----------------------------------------------------------
c x - argument for which interpolation is needed
c f(npf),xf(npf)- table function and argument
c nintr - number of points for interpolation
c
c comment: quite old f77 program but still runs well
c==========================================================
      implicit real*8 (a-h,o-z),integer * 4(i-n)
      dimension f(npf),xf(npf)
c
      if(x.gt.xf(1)) go to 10
        fintr=f(1)
        go to 100
   10 continue
      if(x.lt.xf(npf)) go to 20
        fintr=f(npf)
        go to 100
   20 continue
      if(npf.gt.nintr) go to 25
        i1=1
        i2=npf
        nn=npf
        go to 60
   25 continue
      do 30 i=2,npf
        if(.not. (xf(i-1).le.x  .and.  xf(i).ge.x)) go to 30
          i1=i-1
          i2=i
          go to 35
   30 continue
   35 continue
      id=  int((dfloat(nintr)+0.5)/2.)-1
      i1=i1-id
      i2=i1+nintr-1
cl----------- interpolation is asymmetric at the end
      if(i1.ge.1) go to 40
        i1=1
   40 continue
      if(i2.le.npf) go to 50
        i2=npf
   50 continue
      nn=i2-i1+1
   60 continue
      call polint(f(i1),xf(i1),nn,x,r)
      fintr=r
  100 return
      end

         subroutine polint(f,a,n,x,r)
c----------------------------------------------------------
c service program for fintr
c----------------------------------------------------------
      implicit real*8 (a-h,o-z), integer *4 (i-n)
      dimension f(n),a(n)
      r=0.0
      do 1 j=1,n
      al=1.0
      do 2 i=1,n
         if (i-j) 3,2,3
 3       al=al*(x-a(i))/(a(j)-a(i))
 2    continue
 1    r=r+al*f(j)
      return
      end