I am reading the book Numerical Recipes in Fortran 77: The Art of Scientific Computing (Second Edition) and I came across some methods for numerical integration of 1D functions. More specifically the Gauss-Laguerre, Gauss Hermite, and Gauss Jacobi weights and abscissas appealed to me. I need them because I have some ill-behaved integrands that I need to integrate numerically. As far as I can tell there is no python implementation of these routines, the closest to these methods was this http://docs.scipy.org/doc/numpy/reference/generated/numpy.polynomial.hermite.hermgauss.htmlbut I feel that documentation is a bit incomplete. I'm not sure what to do with this method and slightly puzzled why is there no standard method for these types of quadrature in python.EDIT:I have implemented and compared scipy quadrature with Gauss-Hermite Quadrature on the example problem:\begin{align}I = \int_{-\infty}^{\infty} \left( i\zeta(t) + \overline{i\zeta(t)} \right) \exp(i\omega t) \, \mathrm{d}t\end{align}here $i=\sqrt{-1}$ and $\zeta(t)$ is a complex function of $t$. The overbear denotes complex conjugate. $\zeta(t)$ is not a standard function, it is something I calculated and the samples can be found at http://speedy.sh/SX6JU/seData.npy as a numpy array. Each row in the array is of form [t real(zeta) imag(zeta)] i.e. the sampled real and imaginary values of the function at time t. import numpy as npimport matplotlib.pyplot as pltfrom scipy import integrateT1 =  np.load('seData.npy') #: [Time, x1, x2] from downloaded filedef integrand(t, omega):    x1 = np.interp(t, T1[:,0], T1[:,1])      x2 = np.interp(t, T1[:,0], T1[:,2])        z0 = x1+x2*1j    zeta = z0/np.sqrt(z0**2+1)    return (1j*zeta+np.conjugate(1j*zeta))*np.exp(1j*omega*t) def integral(omega):    def real_func(x,omega):        return np.real(integrand(x,omega))    def imag_func(x,t):        return np.imag(integrand(x,omega))    real_integral = integrate.quad(real_func, -80, 80, args=(omega))    imag_integral = integrate.quad(imag_func, -80, 80, args=(omega))       return real_integral[0] + 1j*imag_integral[0]vintegral = np.vectorize(integral)omega = np.linspace(-30, 30, 601)  I = vintegral(omega)Im not sure what the correct result should be but the noisiness present in the result should not be there I believe. Instead I tried a Gauss-Hermite routine: ti_her, wi_her, = np.polynomial.hermite.hermgauss(90)vintegrand = np.vectorize(integrand)def integral_hermite(omega):    real_integral = np.dot(vintegrand(ti_her, omega).real, wi_her)    imag_integral = np.dot(vintegrand(ti_her, omega).imag, wi_her)    return real_integral + 1j*imag_integralvintegral_hermite = np.vectorize(integral_hermite)I2 = vintegral_hermite(omega)I believe this is more in line with what I should be getting , however I am not sure I have done the correct thing. It seems the order of Gauss-hermite is very high and at low orders (90) it looks like my result shows some type of periodicity/aliasing. 
Is there a Gauss-Laguerre integration routine in Python?
python;interpolation;quadrature