home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / qromb.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  1.1 KB  |  51 lines
  1. #
  2. # qromb
  3. #
  4. # Performs numerical integration of a function, using Romberg's method.
  5. #
  6. # Error 102 indicates that QTRAP_JMAX is too low to get desired accuracy.
  7. #
  8. # This routine is from "Numerical Recipes in C" by Press et. al.
  9. # It uses trapzd(), defined in trapzd.ic, and polint(), defined in polint.ic.
  10. #
  11. # mws, 7/92. 
  12. #
  13. silent
  14. echo "Reading polint..."
  15. read "polint.ic"
  16. echo "Reading trapzd..."
  17. read "trapzd.ic"
  18. echo "Defining qromb..."
  19.  
  20. QROMB_EPS = 1e-6
  21. QROMB_JMAX = 20
  22. QROMB_JMAXP = QROMB_JMAX+1
  23. QROMB_K = 5
  24. array qromb_s[QROMB_JMAXP+1]
  25. array qromb_h[QROMB_JMAXP+1]
  26. array qromb_c[QROMB_K]
  27. array qromb_d[QROMB_K]
  28.  
  29. func qromb(~f,a,b) = {
  30.     local ss,dss,k,j
  31.  
  32.     qromb_h[1]=1
  33.     for (j=1;j<=QROMB_JMAX;j+=1) {
  34.         qromb_s[j]=trapzd(f,a,b,j);
  35.         if (j >= QROMB_K) {
  36.             for (k = 1; k <= QROMB_K; k+=1) {
  37.                 qromb_c[k] = qromb_h[j-QROMB_K+k]
  38.                 qromb_d[k] = qromb_s[j-QROMB_K+k]
  39.             }
  40.             polint(qromb_c,qromb_d,QROMB_K,0,ss,dss)
  41.             if (abs(dss) < QROMB_EPS*abs(ss)) return ss
  42.         }
  43.         qromb_s[j+1]=qromb_s[j]
  44.         qromb_h[j+1]=0.25*qromb_h[j]
  45.     }
  46.     error(102)    # too many iterations required
  47. }
  48.  
  49. echo "Done."
  50. verbose
  51.