The objective is to fit a theoretical spectrum shape, which takes the form of a Lorentzian to a point cloud.
In theory
The Lorentzian is defined by :
(1)\[F(\omega) = \frac{S}{\pi} * \frac{\gamma}{(\omega-\omega_m)^2 + \gamma^2}\]
Where \(\omega\) is the emission frequency; \(\omega_m\) is the “average” frequency, corresponding to the centre of the Lorentzian; \(S = \int_{-\infty}^{\infty} F(\omega) d\omega\) is the total area under the spectrum curve (equivalent to the total energy); \(\gamma\) is the broadening coefficient of the gas under consideration.
In order to perform this theoretical fit, we need to transform this function so that it takes the form of a polynomial. We therefore take its inverse:
(2)\[ \begin{align}\begin{aligned}F(\omega)^{-1} &= \frac{\pi}{S} \frac{(\omega-\omega_m)^2 + \gamma^2}{\gamma}\\& = \frac{\pi}{S\gamma}\omega^2 - \frac{2\pi\omega_m}{S\gamma}\omega + \frac{\pi(\omega_m^2 + \gamma^2)}{S\gamma}\end{aligned}\end{align} \]
For the program, the values of omega are likely to be too large. We then define :
(3)\[ \begin{align}\begin{aligned}x &= \frac{\omega - \bar{\omega}}{\sqrt{\overline{\omega^2} - \bar{\omega}^2}}\\& \equiv \frac{\omega - \bar{\omega}}{\sigma}\end{aligned}\end{align} \]
This gives us :
(4)\[F(\omega)^{-1} = \frac{\pi\gamma^2}{S\gamma}x^2 + \frac{2\pi\sigma}{S\gamma}(\bar{\omega}-\omega_m)x + \frac{\pi}{S\gamma}(\gamma^2 + (\bar{\omega} - \omega_m)^2)\]
Which can be rewritten as:
(5)\[F(\omega)^{-1} = a_2x^2 + a_1x + a_0\]
The objective is to find the combination \({a_0, a_1, a_2}\) so as to minimize the coefficients
(6)\[E(a_0,a_1,a_2) = \sum_{i=1}^{N} W_i (a_2 x_i^2 + a_1 x_i + a_0 - y_i)^2\]
To find a minimum, we therefore look for the points where the derivative cancels, which gives us:
(7)\[\begin{split}\begin{cases}
E_{a_0}' &= 0 \\
E_{a_1}' &= 0 \\
E_{a_2}' &= 0
\end{cases}
\rightarrow
\begin{cases}
\sum_{i=1}^N W_i 2(a_2 x_i^2 + a_1 x_i + a_0 - y_i) &= 0 \\
\sum_{i=1}^N W_i 2(a_2 x_i^2 + a_1 x_i + a_0 - y_i) x_i &= 0 \\
\sum_{i=1}^N W_i 2(a_2 x_i^2 + a_1 x_i + a_0 - y_i) x_i^2 &= 0
\end{cases}\end{split}\]
Which when decomposed gives us:
(8)\[\begin{split}\begin{cases}
a_2 \sum_{i=1}^N W_i x_i^2 + a_1 \sum_{i=1}^N W_i x_i + a_0 \sum_{i=1}^N W_i \\
a_2 \sum_{i=1}^N W_i x_i^3 + a_1 \sum_{i=1}^N W_i x_i^2 + a_0 \sum_{i=1}^N W_i x_i \\
a_2 \sum_{i=1}^N W_i x_i^4 + a_1 \sum_{i=1}^N W_i x_i^3 + a_0 \sum_{i=1}^N W_i x_i^2
\end{cases}
=
\begin{cases}
\sum_{i=1}^N W_i y_i \\
\sum_{i=1}^N W_i y_i x_i \\
\sum_{i=1}^N W_i y_i x_i^2
\end{cases}\end{split}\]
If we divide everything by \(N\), we get the average of all terms:
(9)\[\begin{split}\begin{cases}
a_2 \langle W \rangle \langle x^2 \rangle &+ a_1 \langle W \rangle \langle x \rangle &+ a_0 \langle W \rangle \\
a_2 \langle W \rangle \langle x^3 \rangle &+ a_1 \langle W \rangle \langle x^2 \rangle &+ a_0 \langle W \rangle \langle x \rangle \\
a_2 \langle W \rangle \langle x^4 \rangle &+ a_1 \langle W \rangle \langle x^3 \rangle &+ a_0 \langle W \rangle \langle x^2 \rangle
\end{cases}
=
\begin{cases}
\langle W \rangle \langle y \rangle \\
\langle W \rangle \langle y x \rangle \\
\langle W \rangle \langle y x^2 \rangle
\end{cases}\end{split}\]
Thus, we can divide everything by \(\langle W \rangle\) and write this system as a matrix equation:
(10)\[\begin{split}\begin{pmatrix}
1 & \langle x \rangle & \langle x^2 \rangle \\
\langle x \rangle & \langle x^2 \rangle & \langle x^3 \rangle \\
\langle x^2 \rangle & \langle x^3 \rangle & \langle x^4 \rangle
\end{pmatrix}.
\begin{pmatrix}
a_0 \\
a_1 \\
a_2
\end{pmatrix}=
\begin{pmatrix}
\langle y \rangle \\
\langle y x \rangle \\
\langle y x^2 \rangle
\end{pmatrix}\end{split}\]
So we can determine the coefficients \(a_0, a_1, a_2\)
(11)\[\begin{split}a_0 = \frac{
\begin{vmatrix}
\langle y \rangle & \langle x \rangle & \langle x^2 \rangle \\
\langle yx \rangle & \langle x^2 \rangle & \langle x^3 \rangle \\
\langle yx^2 \rangle & \langle x^3 \rangle & \langle x^4 \rangle
\end{vmatrix}
}{
\begin{vmatrix}
1 & \langle x \rangle & \langle x^2 \rangle \\
\langle x \rangle & \langle x^2 \rangle & \langle x^3 \rangle \\
\langle x^2 \rangle & \langle x^3 \rangle & \langle x^4 \rangle
\end{vmatrix}
}\end{split}\]
(12)\[\begin{split}a_1 &= \frac{
\begin{vmatrix}
\langle 1 \rangle & \langle y \rangle & \langle x^2 \rangle \\
\langle x \rangle & \langle yx \rangle & \langle x^3 \rangle \\
\langle x^2 \rangle & \langle yx^2 \rangle & \langle x^4 \rangle
\end{vmatrix}
}{
\begin{vmatrix}
1 & \langle x \rangle & \langle x^2 \rangle \\
\langle x \rangle & \langle x^2 \rangle & \langle x^3 \rangle \\
\langle x^2 \rangle & \langle x^3 \rangle & \langle x^4 \rangle
\end{vmatrix}
}\end{split}\]
(13)\[\begin{split}a_2 &= \frac{
\begin{vmatrix}
\langle 1 \rangle & \langle x \rangle & \langle y \rangle \\
\langle x \rangle & \langle x^2 \rangle & \langle yx \rangle \\
\langle x^2 \rangle & \langle x^3 \rangle & \langle yx^2 \rangle
\end{vmatrix}
}{
\begin{vmatrix}
1 & \langle x \rangle & \langle x^2 \rangle \\
\langle x \rangle & \langle x^2 \rangle & \langle x^3 \rangle \\
\langle x^2 \rangle & \langle x^3 \rangle & \langle x^4 \rangle
\end{vmatrix}
}\end{split}\]
Once these coefficients have been calculated, the values of \(S\), \(gamma\) and \(\bar{\omega}\) can be found. Their expression can be recovered from the equations (4) and (5)
(14)\[\bar{\omega} = \langle w \rangle - \sigma \frac{a_1}{2 a_2}\]
(15)\[\gamma = \sigma \sqrt{\frac{a_0}{a_2} - \frac{a_1^2}{4 a_2^2}}\]
(16)\[S = \frac{\pi \sigma}{\sqrt{a_0 a_2} - \frac{a_1^2}{4}}\]
With these 3 parameters, we then have a complete description of \(F(\omega)\) as initially described in the equation (1)
In practice
The program reads a file containing on each line a value of \(F(\omega)\). We know that these values are given for a \(\omega\) starting at 2280 and each line increments \(\omega\) by 0.01.
With this information, we have the coordinates of each point. However, we want to invert the coordinate \(y\) because we are interested in \(F(\omega)^{-1}\) in order to have a function in the form of a polynomial.
Also, we calculate for each value of \(\omega\), the associated value of \(x\).
At this point, we have a new set of coordinates that can be exploited. We then calculate the average values of the different elements (\(x\), \(y\) but also their square etc.). The calculation of this average is done via a dedicated routine.
Once this is done, we can now calculate the coefficients \(a_0\), \(a_1\) and \(a_2\).
From this we can calculate the values of \(\bar{\omega}\), \(\gamma\) and \(S\).
Using the first spectrum corresponding to a pressure of 1 atm, we obtain the following result:
However, we notice that the curve we obtain does not totally coincide with our initial data. Indeed, here, we have given an equal weight to each of the experimental values. However, the high values are more significant than the low values because the noise becomes relatively negligible. So, to avoid trying to adjust the curve to the ambient noise, and thus having this kind of error, we give a weight proportional to the squared intensity of each measurement.
In this way, we obtain the following result:
The operation can be repeated for the other spectra, thus for other pressures, which gives us:
