Scipy.integrate Pseudo-voigt Function, Integral Becomes 0

I am writing a script for fitting peak shapes to spectroscopic data in Python, using Scipy, Numpy and Matplotlib. It can fit multiple peaks at once. The peak profile (for now) is P

Solution 1:

As I'm sure you are aware, the interval (-inf, inf) is pretty big. :) A Gaussian decays very fast, so except for an interval near the peak, a Gaussian is numerically indistinguishable from 0. I suspect quad is simply not seeing your peak.

A simple work-around is to split the integral into two intervals, (-inf, pos) and (pos, inf). (Your function is symmetric about Pos, so you really only need twice the integral over (-inf, pos).)

Here's an example. I don't know if these parameter values are anywhere near the typical values you use, but they illustrate the point.

In [259]: pos = 1500.0

In [260]: amp = 4.0

In [261]: gammal = 0.5

In [262]: fracl = 0# Pure Gaussian

quad thinks the integral is 0:

In [263]: quad(PseudoVoigtFunction, -np.inf, np.inf, args=(pos, amp, gammal, fracl))
Out[263]: (0.0, 0.0)

Instead, integrate over (-inf, pos) and (pos, inf):

In[264]: quad(PseudoVoigtFunction, -np.inf, pos, args=(pos, amp, gammal, fracl))
Out[264]: (1.5053836955785238, 3.616268258191726e-11)

In[265]: quad(PseudoVoigtFunction, pos, np.inf, args=(pos, amp, gammal, fracl))
Out[265]: (1.5053836955785238, 3.616268258191726e-11)

So the integral over (-inf, inf) is approximately 3.010767.