simpsn.c 1.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081
  1. /* simpsn.c */
  2. /* simpsn.c
  3. * Numerical integration of function tabulated
  4. * at equally spaced arguments
  5. */
  6. /* Coefficients for Cote integration formulas */
  7. /* Note: these numbers were computed using 40-decimal precision. */
  8. #define NCOTE 8
  9. /* 6th order formula */
  10. /*
  11. static double simcon[] =
  12. {
  13. 4.88095238095238095E-2,
  14. 2.57142857142857142857E-1,
  15. 3.2142857142857142857E-2,
  16. 3.2380952380952380952E-1,
  17. };
  18. */
  19. /* 8th order formula */
  20. static double simcon[] =
  21. {
  22. 3.488536155202821869E-2,
  23. 2.076895943562610229E-1,
  24. -3.27336860670194003527E-2,
  25. 3.7022927689594356261E-1,
  26. -1.6014109347442680776E-1,
  27. };
  28. /* 10th order formula */
  29. /*
  30. static double simcon[] =
  31. {
  32. 2.68341483619261397039E-2,
  33. 1.77535941424830313719E-1,
  34. -8.1043570626903960237E-2,
  35. 4.5494628827962161295E-1,
  36. -4.3515512265512265512E-1,
  37. 7.1376463043129709796E-1,
  38. };
  39. */
  40. /* simpsn.c 2 */
  41. /* 20th order formula */
  42. /*
  43. static double simcon[] =
  44. {
  45. 1.182527324903160319E-2,
  46. 1.14137717644606974987E-1,
  47. -2.36478370511426964E-1,
  48. 1.20618689348187566E+0,
  49. -3.7710317267153304677E+0,
  50. 1.03367982199398011435E+1,
  51. -2.270881584397951229796E+1,
  52. 4.1828057422193554603E+1,
  53. -6.4075279490154004651555E+1,
  54. 8.279728347247285172085E+1,
  55. -9.0005367135242894657916E+1,
  56. };
  57. */
  58. /* simpsn.c 3 */
  59. double simpsn( f, delta )
  60. double f[]; /* tabulated function */
  61. double delta; /* spacing of arguments */
  62. {
  63. extern double simcon[];
  64. double ans;
  65. int i;
  66. ans = simcon[NCOTE/2] * f[NCOTE/2];
  67. for( i=0; i < NCOTE/2; i++ )
  68. ans += simcon[i] * ( f[i] + f[NCOTE-i] );
  69. return( ans * delta * NCOTE );
  70. }