Summary: We present a code to compute the oscillatory cuspoid canonical integrals and their first order partial derivatives. The algorithm is based on the method of J.N.L. Conno
and P.R. Curtis
[J. Phys. A 15, 1179–1190 (1982; Zbl 0485.65019
)], in which the integration path along the real axis is replaced by a more convenient contour in the complex plane, rendering the oscillatory integrand more amenable to numerical quadrature. Our code is a modern implementation of this method, presented in a modular fashion as a Fortran 90 module containing the relevant subroutines. The code is robust, with the novel feature that the algorithm implements an adaptive contour procedure, choosing contours that avoid the violent oscillatory and exponential natures of the integrand and modifying its choice as necessary. In many cases, the subroutines are efficient and provide results of high accuracy. Plots and results for the Pearcey and swallowtail cuspoid integrals given in previous articles are reproduced and verified. Our subroutines significantly extend the range of application of this earlier work.