custom `Python3` code which exploits a similar approach to that of the `HFS` method implemented in CLASS. Differently from the latter, our code uses the column density ($N$) and the excitation temperature ($T_{\ensuremath{\mathrm{ex}}}$) as fit parameters, whereas the total opacity of the transition ($\tau$) is regarded as a derived quantity. To this aim, the routine makes use of the physical line parameters and implements the computation of the rotational partition function for any given $T_{\ensuremath{\mathrm{ex}}}$. The set of the fit parameters also includes the line systemic velocity ($\varv_0$) in the local standard of rest, and the full-width half maximum of the lines (FWHM, $\Delta\varv$). Briefly, for a given input the total line opacity is first calculated through: $$\label{eq:tau}
        \tau = \sqrt{\frac{\ln 2}{16\pi^3}}\:
        \frac{c^3  A_{ul} g_u N}{\nu^3 Q_{\ensuremath{\mathrm{rot}}}(T_{\ensuremath{\mathrm{ex}}}) \Delta\varv} 
        \exp\left(\frac{-E_{\ensuremath{\mathrm{low}}}}{k_{\ensuremath{\mathrm{B}}} T_\text{ex}}\right) 
        \left[1 - \exp\left(\frac{-h\nu_0}{k_{\ensuremath{\mathrm{B}}} T_\text{ex}} \right)\right] \,,$$

where $\nu_0$ is the rest frequency of the (unsplit) rotational transition, $c$ is the speed of light, $Q_{\ensuremath{\mathrm{rot}}}(T_{\ensuremath{\mathrm{ex}}})$ is the rotational partition function at $T_{\ensuremath{\mathrm{ex}}}$, $g_u$ is the degeneracy of the upper level, $A_{ul}$ is the Einstein coefficient for spontaneous emission, $E_{\ensuremath{\mathrm{low}}}$ is the energy of the lower state, $h$ is the Planck constant, and $k_{\ensuremath{\mathrm{B}}}$ is the Boltzmann constant. Then, assuming the local 