Name

---

      TRAPZD

Purpose

---

      Compute the nth stage of refinement of an extended trapezoidal rule.

Explanation

---

      This procedure is called by QSIMP and QTRAP. Algorithm from Numerical
      Recipes, Section 4.2. TRAPZD is meant to be called iteratively from
      a higher level procedure.

Calling Sequence

---

      TRAPZD, func, A, B, S, step, [ _EXTRA = ]

Inputs

---

      func - scalar string giving name of function to be integrated. This
              must be a function of one variable.
      A,B - scalars giving the limits of the integration

Input-output

---

      S - scalar giving the total sum from the previous iterations on
              input and the refined sum after the current iteration on output.
      step - LONG scalar giving the number of points at which to compute the
              function for the current iteration. If step is not defined on
              input, then S is intialized using the average of the endpoints
              of limits of integration.

Optional Input Keywords

---

      Any supplied keywords will be passed to the user function via the
      _EXTRA facility.

Notes

---

      (1) TRAPZD will check for math errors (except for underflow) when
      computing the function at the endpoints, but not on subsequent
      iterations.
      (2) TRAPZD always uses double precision to sum the function values
      but the call to the user-supplied function is double precision only if
      one of the limits A or B is double precision.

Revision History

---

      Written W. Landsman August, 1991
      Always use double precision for TOTAL March, 1996
      Pass keyword to function via _EXTRA facility W. Landsman July 1999
      Don't check for floating underflow W.Landsman April 2008