You can use recursion.
Also from scipy.integrate.simps it says the assumed axis is the last axis.
import numpy
import scipy.integrate
def intNd(c,axes):
''' c is the array
axes is a list of the corresponding coordinates
'''
assert len(c.shape) == len(axes)
assert all([c.shape[i] == axes[i].shape[0]
for i in range(len(axes))])
if len(axes) == 1:
return scipy.integrate.simps(c,axes[0])
else:
return intNd(scipy.integrate.simps(c,axes[-1]),axes[:-1])
You may also consider improving the efficiency for integrating large complex arrays by integrating along the longest dimension first. However. not integrating along the last axis, which is the fastest, may impose some penalty. I did not look into these details myself.