integration

Below is a function that (i) interpolates an array using a Piecewise Cubic Hermitian Interpolator (using DPCHEZ), and then (ii) integrates the interpolation (using DPCHQA). The required external methods are provided by the PCHIP module.

This method is quite accurate, but not as fast as e.g. the trapezoid method below.

```
function integrate(x, y, a, b) result(r)
!! This function constructs a piecewise cubic Hermitian interpolation of an array y(x) based on
!! discrete numerical data, and subsequently evaluates the integral of the interpolation in the
!! range (a,b). Note that the mesh spacing of x does not necessarily have to be uniform.
real(wp), intent(in) :: x(:) !! Variable x
real(wp), intent(in) :: y(size(x)) !! Function y(x)
real(wp), intent(in) :: a !! Left endpoint
real(wp), intent(in) :: b !! Right endpoint
real(wp) :: r !! Integral ∫y(x)·dx
real(wp), external :: dpchqa
real(wp) :: d(size(x))
integer :: err
! Create a PCHIP interpolation of the input data
call dpchez(size(x), x, y, d, .false., 0, 0, err)
! Integrate the interpolation in the provided range
r = dpchqa(size(x), x, y, d, a, b, err)
end function
```

Below is a simple function for numerically calculating the integral of an an array using the trapezoid method. Since it is expressed using whole-array operations, a good compiler should be able to vectorize it automatically, rendering it very fast.

```
pure function integrate(x, y) result(r)
!! Calculates the integral of an array y with respect to x using the trapezoid
!! approximation. Note that the mesh spacing of x does not have to be uniform.
real(wp), intent(in) :: x(:) !! Variable x
real(wp), intent(in) :: y(size(x)) !! Function y(x)
real(wp) :: r !! Integral ∫y(x)·dx
! Integrate using the trapezoidal rule
associate(n => size(x))
r = sum((y(1+1:n-0) + y(1+0:n-1))*(x(1+1:n-0) - x(1+0:n-1)))/2
end associate
end function
```