Scipy.integrate Псевдо-фойгтовая функция, интеграл становится 0 ⇐ Python

[Anonymous]

Я пишу скрипт для подгонки формы пиков к спектроскопическим данным на Python, используя Scipy, Numpy и Matplotlib. Он может соответствовать нескольким вершинам одновременно. Пиковый профиль (на данный момент) — это псевдо-Фойгт, который представляет собой линейную комбинацию гауссова (так называемого нормального) и лоренцева (также известного как Коши) распределения.

У меня есть переключатель, с помощью которого я могу либо позволить программному обеспечению оптимизировать вклад гауссова и лоренциана, либо установить для него фиксированное значение (где 0 = чистый гауссиан и 1 = чистый лоренциан). Работает так, как должно, отображение подобранных пиков выглядит так, как ожидалось. Проблема начинается, когда я пытаюсь вычислить интегралы пиков, используя

Код: Выделить всё

scipy.integrate

.

So far I tried scipy.integrate.quad, scipy.integrate.quadrature, scipy.integrate.fixed_quad and scipy.integrate.romberg. The integral becomes something like

Код: Выделить всё

1.73476E-34

(not always the same number) when the peak is pure Gaussian, even for peaks that clearly have larger areas than neighbouring peaks that are not pure Gaussian but return a finite integral of something on the order of 10 to 1000. Here's what the relevant parts look like:

Код: Выделить всё

# Function defining the peak functions for plotting and integration
# WavNr: Wave number, the x-axis over which shall be integrated
# Pos: Peak center position
# Amp: Amplitude of the peak
# GammaL: Gamma parameter of the Lorentzian distribution
# FracL: Fraction of Lorentzian distribution
def PseudoVoigtFunction(WavNr, Pos, Amp, GammaL, FracL):
SigmaG = GammaL / np.sqrt(2*np.log(2)) # Calculate the sigma parameter  for the Gaussian distribution from GammaL (coupled in Pseudo-Voigt)
LorentzPart = Amp * (GammaL**2 / ((WavNr - Pos)**2 + GammaL**2)) # Lorentzian distribution
GaussPart = Amp * np.exp( -((WavNr - Pos)/SigmaG)**2) # Gaussian distribution
Fit = FracL * LorentzPart + (1 - FracL) * GaussPart # Linear combination of the two parts (or distributions)
return Fit

This is the called by the plot function via:

Код: Выделить всё

Fit = PseudoVoigtFunction(WavNr, Pos, Amp, GammaL, FracL)

which works fine. It is also called by the integrator via:

Код: Выделить всё

PeakArea, PeakAreaError = integrate.quad(PseudoVoigtFunction, -np.inf, np.inf, args=(Pos, Amp, GammaL, FracL))

or any of the other variants supplied by scipy.integrate, all with the same result that if FracL = 0, then PeakArea = (almost) 0.

I'm sure the problem is just me being too stupid to figure out how scipy.integrate works with slightly more complicated functions than I can find examples for. Hoping someone sees the obvious error that I don't. Two days of searching stackoverflow and Scipy Docs and rearranging and completely rewriting my code have taken me nowhere. I'm suspecting that the args in the scipy.integrate are somehow related to the problem but for all I could find out, they seem to be arranged correctly.

Thanks in advance,
Os


Источник: https://stackoverflow.com/questions/376 ... -becomes-0