home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / polint.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  1.1 KB  |  58 lines
  1. #
  2. # Routine for polynomial interpolation from "Numerical Recipes in C",
  3. # by Press et. al.
  4. #
  5. # This routine is used by qromb.ic
  6. #
  7. # mws, 7/92
  8. #
  9.  
  10. #
  11. # polint
  12. #
  13. # Given arrays xa[1..n] and ya[1..n], and given a value x, this routine
  14. # returns a value y, and an error estimate dy. If P(x) is the polynomial of
  15. # degree n-1 s.t. P(xa[j]) = ya[j] for j = 1..n, then the returned value is
  16. # y = P(x).
  17. #
  18. array polint_c[1]    # will be resized as required
  19. array polint_d[1]
  20. #
  21. func polint(@xa,@ya,n,x,*y,*dy) = {
  22.     local ii,m,den,dif,dift,ho,hp,w,ns
  23.  
  24.     ns = 1
  25.     dif=abs(x-xa[1])
  26.  
  27.     resize(polint_c,n)
  28.     resize(polint_d,n)
  29.  
  30.     for (ii=1;ii<=n;ii+=1) {
  31.         if ( (dift=abs(x-xa[ii])) < dif) {
  32.             ns=ii
  33.             dif=dift
  34.         }
  35.         polint_c[ii]=ya[ii]
  36.         polint_d[ii]=ya[ii]
  37.     }
  38.     y=ya[ns]
  39.     ns -= 1
  40.     for (m=1;m<n;m+=1) {
  41.         for (ii=1;ii<=n-m;ii+=1) {
  42.             ho=xa[ii]-x
  43.             hp=xa[ii+m]-x
  44.             w=polint_c[ii+1]-polint_d[ii]
  45.             # error can occur only if two input xa's are identical
  46.             if ( (den=ho-hp) == 0) error(0)
  47.             den=w/den
  48.             polint_d[ii]=hp*den
  49.             polint_c[ii]=ho*den
  50.         }
  51.         if (2*ns < (n-m)) dy = polint_c[ns+1] else {
  52.             dy = polint_d[ns]
  53.             ns -= 1
  54.         }
  55.         y += dy
  56.     }
  57. }
  58.