Example #3: non-equilibrium spectrum (Tvib, Trot, x_CO)

With the new fitting module introduced in fit_spectrum() function, and in the example of Tgas fitting using new fitting module, we can see its 1-temperature fitting performance for equilibrium condition.

This example features how new fitting module can fit non-equilibrium spectrum, with multiple fit parameters, such as vibrational/rotational temperatures, or mole fraction, etc.

This is a real fitting case introduced Grimaldi’s thesis: https://doi.org/10.2514/1.T6768 . This case features CO molecule emitting on a wide range of spectrum. This case also includes a user-defined trapezoid slit function, which is accounted by the distortion (dispersion) of the slit, as a result of the influences from experimental parameters and spectrometer dispersion parameters during the experiment.

• 
• 

======================= COMMENCE FITTING PROCESS =======================

Successfully retrieved the experimental data in 0.0007104873657226562s.

Acquired spectral quantity 'radiance' from the spectrum.
NaN values successfully purged.Number of data points left: 1479 points.
Successfully refined the experimental data in 0.0002532005310058594s.
0.19s - Loaded database

Commence fitting process for non-LTE spectrum!
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/warning.py:434: PerformanceWarning: 'gu' was recomputed although 'gp' already in DataFrame. All values are equal
  warnings.warn(WarningType(message))
[[Fit Statistics]]
    # fitting method   = L-BFGS-B
    # function evals   = 96
    # data points      = 1
    # variables        = 3
    chi-square         = 4.1191e-21
    reduced chi-square = 4.1191e-21
    Akaike info crit   = -40.9386489
    Bayesian info crit = -46.9386489
##  Warning: uncertainties could not be estimated:
    this fitting method does not natively calculate uncertainties
    and numdifftools is not installed for lmfit to do this. Use
    `pip install numdifftools` for lmfit to estimate uncertainties
    with this fitting method.
[[Variables]]
    Tvib:           4019.23957 (init = 6000)
    Trot:           3213.32473 (init = 4000)
    mole_fraction:  0.06307478 (init = 0.05)

Successfully finished the fitting process in 7.707149982452393s.
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/misc/curve.py:240: UserWarning: Presence of NaN in curve_divide!
Think about interpolation=2
  warnings.warn(

======================== END OF FITTING PROCESS ========================


Residual history:

[np.float64(1.6884210236199504e-10), np.float64(1.6884210440381979e-10), np.float64(1.6884210060033468e-10), np.float64(1.6884210671365353e-10), np.float64(3.738560781412175e-10), np.float64(3.7385607749613685e-10), np.float64(3.7385607855679e-10), np.float64(3.7385607553763786e-10), np.float64(1.1128815372128413e-10), np.float64(1.1128815390646927e-10), np.float64(1.112881541942443e-10), np.float64(1.1128815294674758e-10), np.float64(1.0954119400145755e-10), np.float64(1.0954119450442979e-10), np.float64(1.0954119419457839e-10), np.float64(1.0954119402681508e-10), np.float64(1.0815468208401315e-10), np.float64(1.0815468267182474e-10), np.float64(1.0815468221710159e-10), np.float64(1.08154682280824e-10), np.float64(1.0308407647672229e-10), np.float64(1.0308407742635263e-10), np.float64(1.0308407637797923e-10), np.float64(1.0308407728109503e-10), np.float64(2.5622069956334375e-10), np.float64(2.562207030140033e-10), np.float64(2.562206993262903e-10), np.float64(2.562206988928061e-10), np.float64(9.533544922549783e-11), np.float64(9.53354504862476e-11), np.float64(9.533544899450552e-11), np.float64(9.533545030593561e-11), np.float64(8.48427735600552e-11), np.float64(8.484277492938125e-11), np.float64(8.484277332533733e-11), np.float64(8.484277453072342e-11), np.float64(2.942162855844573e-10), np.float64(2.9421628239932186e-10), np.float64(2.9421628625246667e-10), np.float64(2.942162845188683e-10), np.float64(7.545778695626082e-11), np.float64(7.545778797764293e-11), np.float64(7.545778688540757e-11), np.float64(7.545778741530068e-11), np.float64(6.932656171835516e-11), np.float64(6.932656034728585e-11), np.float64(6.932656236751499e-11), np.float64(6.932656036175411e-11), np.float64(7.949298283128666e-11), np.float64(7.949298422503194e-11), np.float64(7.949298256334672e-11), np.float64(7.949298374088136e-11), np.float64(6.81498435553155e-11), np.float64(6.81498430188567e-11), np.float64(6.814984397170331e-11), np.float64(6.814984272733105e-11), np.float64(6.719214857166717e-11), np.float64(6.71921481791797e-11), np.float64(6.719214886698425e-11), np.float64(6.719214790974426e-11), np.float64(6.549161955740076e-11), np.float64(6.549161954424184e-11), np.float64(6.549161943871277e-11), np.float64(6.54916193652463e-11), np.float64(6.525776504687181e-11), np.float64(6.525776505808518e-11), np.float64(6.525776491546529e-11), np.float64(6.525776489787045e-11), np.float64(6.423320929399242e-11), np.float64(6.42332091556419e-11), np.float64(6.423320930784754e-11), np.float64(6.423320919680103e-11), np.float64(6.418825495583594e-11), np.float64(6.418825503117023e-11), np.float64(6.418825493438653e-11), np.float64(6.418825499331539e-11), np.float64(6.418046441427161e-11), np.float64(6.418046440587178e-11), np.float64(6.41804644195429e-11), np.float64(6.41804644093244e-11), np.float64(6.418033839618103e-11), np.float64(6.418033839454034e-11), np.float64(6.418033839742535e-11), np.float64(6.418033839463614e-11), np.float64(6.418032597727522e-11), np.float64(6.418032597694697e-11), np.float64(6.418032597759614e-11), np.float64(6.418032597681513e-11), np.float64(6.418032367206344e-11), np.float64(6.418032367209917e-11), np.float64(6.418032367204627e-11), np.float64(6.418032367207511e-11), np.float64(6.418032366726698e-11), np.float64(6.418032366727002e-11), np.float64(6.41803236672653e-11), np.float64(6.418032366726851e-11), np.float64(6.418032366726698e-11)]

Fitted values history:

[6000.0, 4000.0, 0.05]
[6000.00002, 4000.0, 0.05]
[6000.0, 4000.000024494898, 0.05]
[6000.0, 4000.0, 0.0500000005]
[5105.589014486926, 4855.661714738742, 0.012435519427667975]
[5105.589038742363, 4855.661714738742, 0.012435519427667975]
[5105.589014486926, 4855.661739484458, 0.012435519427667975]
[5105.589014486926, 4855.661714738742, 0.012435519757654313]
[5818.190228435902, 4187.474367518324, 0.04066587313976969]
[5818.190249678251, 4187.474367518324, 0.04066587313976969]
[5818.190228435902, 4187.4743923222095, 0.04066587313976969]
[5818.190228435902, 4187.474367518324, 0.040665873630979836]
[5807.164848659225, 4147.48038979007, 0.041997575235936876]
[5807.164869969598, 4147.48038979007, 0.041997575235936876]
[5807.164848659225, 4147.480414540282, 0.041997575235936876]
[5807.164848659225, 4147.48038979007, 0.04199757572949145]
[5757.359686017544, 4111.453799922989, 0.04251231193185258]
[5757.359707625521, 4111.453799922989, 0.04251231193185258]
[5757.359686017544, 4111.453824619207, 0.04251231193185258]
[5757.359686017544, 4111.453799922989, 0.04251231242621424]
[5551.640562063738, 3968.2486675433383, 0.04457882915299982]
[5551.64058474424, 3968.2486675433383, 0.04457882915299982]
[5551.640562063738, 3968.248691971274, 0.04457882915299982]
[5551.640562063738, 3968.2486675433383, 0.04457882965005222]
[2363.2518437831704, 2028.035310886635, 0.0881500560958061]
[2363.2518308050835, 2028.035310886635, 0.0881500560958061]
[2363.2518437831704, 2028.0353071531306, 0.0881500560958061]
[2363.2518437831704, 2028.035310886635, 0.08815005641900468]
[5270.46023593549, 3795.36469150175, 0.046884234996570164]
[5270.4602597186595, 3795.36469150175, 0.046884234996570164]
[5270.46023593549, 3795.364715488182, 0.046884234996570164]
[5270.46023593549, 3795.36469150175, 0.04688423549559842]
[4949.513430682526, 3610.0314990484503, 0.0494185411328201]
[4949.5134552750815, 3610.0314990484503, 0.0494185411328201]
[4949.513430682526, 3610.031522410719, 0.0494185411328201]
[4949.513430682526, 3610.0314990484503, 0.049418541632786285]
[2440.7348269315166, 2341.0416722030936, 0.07075974094128255]
[2440.7348411069424, 2341.0416722030936, 0.07075974094128255]
[2440.7348269315166, 2341.041684808247, 0.07075974094128255]
[2440.7348269315166, 2341.0416722030936, 0.07075974139614881]
[4677.684247042535, 3460.717135614628, 0.05151738666294371]
[4677.684271979313, 3460.717135614628, 0.05151738666294371]
[4677.684247042535, 3460.7171583520226, 0.05151738666294371]
[4677.684247042535, 3460.717135614628, 0.05151738716271342]
[4198.138537235358, 3209.915406743408, 0.05518699273191109]
[4198.138562052449, 3209.915406743408, 0.05518699273191109]
[4198.138537235358, 3209.915428157613, 0.05518699273191109]
[4198.138537235358, 3209.915406743408, 0.055186993229213326]
[4786.938383930233, 3491.189773422493, 0.050731820854976586]
[4786.93840876502, 3491.189773422493, 0.050731820854976586]
[4786.938383930233, 3491.1897962967187, 0.050731820854976586]
[4786.938383930233, 3491.189773422493, 0.050731821354923025]
[4316.412798024692, 3265.031847683382, 0.05429200564275049]
[4316.412822957192, 3265.031847683382, 0.05429200564275049]
[4316.412798024692, 3265.0318694201105, 0.05429200564275049]
[4316.412798024692, 3265.031847683382, 0.05429200614090495]
[4284.404479920046, 3228.845901782738, 0.05493288400739585]
[4284.404504826909, 3228.845901782738, 0.05493288400739585]
[4284.404479920046, 3228.8459233098542, 0.05493288400739585]
[4284.404479920046, 3228.845901782738, 0.054932884504956564]
[4145.561172586466, 3098.5223656113867, 0.05762044917567394]
[4145.561197333936, 3098.5223656113867, 0.05762044917567394]
[4145.561172586466, 3098.522386313707, 0.05762044917567394]
[4145.561172586466, 3098.5223656113867, 0.0576204496698327]
[4126.632862305605, 3101.825182110372, 0.05824217878616823]
[4126.632887025227, 3101.825182110372, 0.05824217878616823]
[4126.632862305605, 3101.8252028350125, 0.05824217878616823]
[4126.632862305605, 3101.825182110372, 0.0582421792793281]
[4017.831626970656, 3184.5833354398346, 0.06249112747735987]
[4017.831651501277, 3184.5833354398346, 0.06249112747735987]
[4017.831626970656, 3184.5833566993715, 0.06249112747735987]
[4017.831626970656, 3184.5833354398346, 0.06249112796150568]
[4030.6305844299595, 3216.4136047497223, 0.06286712013815833]
[4030.6306089853915, 3216.4136047497223, 0.06286712013815833]
[4030.6305844299595, 3216.4136262029415, 0.06286712013815833]
[4030.6305844299595, 3216.4136047497223, 0.06286712062131848]
[4019.035314678171, 3214.7623983051985, 0.06308862400585541]
[4019.035339211154, 3214.7623983051985, 0.06308862400585541]
[4019.035314678171, 3214.7624197485293, 0.06308862400585541]
[4019.035314678171, 3214.7623983051985, 0.06308862448842022]
[4019.701570616155, 3213.6181185448568, 0.06306124382660912]
[4019.701595150444, 3213.6181185448568, 0.06306124382660912]
[4019.701570616155, 3213.618139981325, 0.06306124382660912]
[4019.701570616155, 3213.6181185448568, 0.0630612443092481]
[4019.4652564456005, 3213.397863673256, 0.06306791743856699]
[4019.465280979427, 3213.397863673256, 0.06306791743856699]
[4019.4652564456005, 3213.397885108402, 0.06306791743856699]
[4019.4652564456005, 3213.397863673256, 0.06306791792118792]
[4019.248506097908, 3213.3204797019584, 0.06307448349188255]
[4019.24853063131, 3213.3204797019584, 0.06307448349188255]
[4019.248506097908, 3213.3205011366404, 0.06307448349188255]
[4019.248506097908, 3213.3204797019584, 0.06307448397448569]
[4019.2395650969092, 3213.3247316929546, 0.06307478479743257]
[4019.2395896302933, 3213.3247316929546, 0.06307478479743257]
[4019.2395650969092, 3213.324753127662, 0.06307478479743257]
[4019.2395650969092, 3213.3247316929546, 0.06307478528003489]
[4019.2395650969092, 3213.3247316929546, 0.06307478479743257]

Total fitting time:
7.707149982452393 s

import astropy.units as u
import numpy as np

from radis import load_spec
from radis.test.utils import getTestFile
from radis.tools.new_fitting import fit_spectrum

# The following function and parameters reproduce the slit function in Dr. Corenti Grimaldi's thesis.


def slit_dispersion(w):
    phi = -6.33
    f = 750
    gr = 300
    m = 1
    phi *= -2 * np.pi / 360
    d = 1e-3 / gr
    disp = (
        w
        / (2 * f)
        * (-np.tan(phi) + np.sqrt((2 * d / m / (w * 1e-9) * np.cos(phi)) ** 2 - 1))
    )
    return disp  # nm/mm


slit = 1500  # µm
pitch = 20  # µm
top_slit_um = slit - pitch  # µm
base_slit_um = slit + pitch  # µm
center_slit = 5090
dispersion = slit_dispersion(center_slit)
top_slit_nm = top_slit_um * 1e-3 * dispersion
base_slit_nm = base_slit_um * 1e-3 * dispersion * 1.33


# -------------------- Step 1. Load experimental spectrum -------------------- #


specName = (
    "Corentin_0_100cm_DownSampled_20cm_10pctCO2_1-wc-gw450-gr300-sl1500-acc5000-.spec"
)
s_experimental = load_spec(getTestFile(specName)).offset(-0.6, "nm")


# -------------------- Step 2. Fill ground-truths and data -------------------- #


# Experimental conditions which will be used for spectrum modeling. Basically, these are known ground-truths.
experimental_conditions = {
    "molecule": "CO",  # Molecule ID
    "isotope": "1,2,3",  # Isotope ID, can have multiple at once
    "wmin": 2270
    * u.nm,  # Starting wavelength/wavenumber to be cropped out from the original experimental spectrum.
    "wmax": 2400 * u.nm,  # Ending wavelength/wavenumber for the cropping range.
    "pressure": 1.01325
    * u.bar,  # Total pressure of gas, in "bar" unit by default, but you can also use Astropy units.
    "path_length": 1
    / 195
    * u.cm,  # Experimental path length, in "cm" unit by default, but you can also use Astropy units.
    "slit": {  # Experimental slit. In simple form: "[value] [unit]", i.e. "-0.2 nm". In complex form: a dict with parameters of apply_slit()
        "slit_function": (top_slit_nm, base_slit_nm),
        "unit": "nm",
        "shape": "trapezoidal",
        "center_wavespace": center_slit,
        "slit_dispersion": slit_dispersion,
        "inplace": False,
    },
    "cutoff": 0,  # (RADIS native) Discard linestrengths that are lower that this to reduce calculation time, in cm-1.
    "databank": "hitemp",  # Databank used for calculation. Must be stated.
}

# List of parameters to be fitted, accompanied by its initial value
fit_parameters = {
    "Tvib": 6000,  # Vibrational temperature, in K.
    "Trot": 4000,  # Rotational temperature, in K.
    "mole_fraction": 0.05,  # Species mole fraction, from 0 to 1.
}

# List of bounding ranges applied for those fit parameters above.
# You can skip this step and let it use default bounding ranges, but this is not recommended.
# Bounding range must be at format [<lower bound>, <upper bound>].
bounding_ranges = {
    "Tvib": [2000, 7000],
    "Trot": [2000, 7000],
    "mole_fraction": [0, 0.1],
}

# Fitting pipeline setups.
fit_properties = {
    "method": "lbfgsb",  # Preferred fitting method. By default, "leastsq".
    "fit_var": "radiance",  # Spectral quantity to be extracted for fitting process, such as "radiance", "absorbance", etc.
    "normalize": False,  # Either applying normalization on both spectra or not.
    "max_loop": 300,  # Max number of loops allowed. By default, 200.
    "tol": 1e-20,  # Fitting tolerance, only applicable for "lbfgsb" method.
}

"""

For the fitting method, you can try one among 17 different fitting methods and algorithms of LMFIT,
introduced in `LMFIT method list <https://lmfit.github.io/lmfit-py/fitting.html#choosing-different-fitting-methods>`.

You can see the benchmark result of these algorithms here:
`RADIS Newfitting Algorithm Benchmark <https://github.com/radis/radis-benchmark/blob/master/manual_benchmarks/plot_newfitting_comparison_algorithm.py>`.

"""


# -------------------- Step 3. Run the fitting and retrieve results -------------------- #

# Conduct the fitting process!
s_best, result, log = fit_spectrum(
    s_exp=s_experimental,  # Experimental spectrum.
    fit_params=fit_parameters,  # Fit parameters.
    bounds=bounding_ranges,  # Bounding ranges for those fit parameters.
    model=experimental_conditions,  # Experimental ground-truths conditions.
    pipeline=fit_properties,  # # Fitting pipeline references.
)

# Now investigate the result logs for additional information about what's going on during the fitting process

print("\nResidual history: \n")
print(log["residual"])

print("\nFitted values history: \n")
for fit_val in log["fit_vals"]:
    print(fit_val)

print("\nTotal fitting time: ")
print(log["time_fitting"], end=" s\n")