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.0006387233734130859s.

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.0002295970916748047s.
0.17s - Loaded database

Commence fitting process for non-LTE spectrum!
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/develop/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.23956 (init = 6000)
    Trot:           3213.32471 (init = 4000)
    mole_fraction:  0.06307478 (init = 0.05)

Successfully finished the fitting process in 5.624126434326172s.
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/develop/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.6884210236199499e-10), np.float64(1.6884210440381976e-10), np.float64(1.6884210060033468e-10), np.float64(1.688421067136535e-10), np.float64(3.7385607833419127e-10), np.float64(3.738560776891107e-10), np.float64(3.7385607874976396e-10), np.float64(3.7385607573061175e-10), np.float64(1.1128815368154832e-10), np.float64(1.1128815386673341e-10), np.float64(1.1128815415450849e-10), np.float64(1.1128815290701178e-10), np.float64(1.0954119407604127e-10), np.float64(1.095411945790136e-10), np.float64(1.0954119426916211e-10), np.float64(1.0954119410139886e-10), np.float64(1.0815468217516548e-10), np.float64(1.081546827629771e-10), np.float64(1.0815468230825383e-10), np.float64(1.0815468237197641e-10), np.float64(1.0308407699815074e-10), np.float64(1.0308407794778126e-10), np.float64(1.0308407689940763e-10), np.float64(1.0308407780252369e-10), np.float64(2.5622137633580595e-10), np.float64(2.5622137978645347e-10), np.float64(2.5622137609875424e-10), np.float64(2.5622137566526895e-10), np.float64(9.533545352732892e-11), np.float64(9.533545478807875e-11), np.float64(9.533545329633652e-11), np.float64(9.533545460776691e-11), np.float64(8.484278666606685e-11), np.float64(8.484278803539322e-11), np.float64(8.484278643134863e-11), np.float64(8.484278763673576e-11), np.float64(2.9421603805650496e-10), np.float64(2.942160348713639e-10), np.float64(2.942160387245155e-10), np.float64(2.942160369909142e-10), np.float64(7.545779416178807e-11), np.float64(7.545779518317105e-11), np.float64(7.54577940909344e-11), np.float64(7.545779462082892e-11), np.float64(6.932655577026052e-11), np.float64(6.932655439919306e-11), np.float64(6.932655641941991e-11), np.float64(6.932655441366059e-11), np.float64(7.949295401768293e-11), np.float64(7.949295541142762e-11), np.float64(7.949295374974318e-11), np.float64(7.949295492727656e-11), np.float64(6.814984180974783e-11), np.float64(6.814984127328906e-11), np.float64(6.814984222613549e-11), np.float64(6.814984098176355e-11), np.float64(6.719214701336819e-11), np.float64(6.71921466208809e-11), np.float64(6.719214730868505e-11), np.float64(6.719214635144555e-11), np.float64(6.549161959461511e-11), np.float64(6.5491619581456e-11), np.float64(6.549161947592724e-11), np.float64(6.549161940246052e-11), np.float64(6.525776498973532e-11), np.float64(6.525776500094872e-11), np.float64(6.525776485832874e-11), np.float64(6.525776484073401e-11), np.float64(6.423320863605763e-11), np.float64(6.423320849770837e-11), np.float64(6.423320864991242e-11), np.float64(6.423320853886716e-11), np.float64(6.418825491123364e-11), np.float64(6.418825498656769e-11), np.float64(6.418825488978427e-11), np.float64(6.4188254948713e-11), np.float64(6.418046440995909e-11), np.float64(6.418046440155933e-11), np.float64(6.418046441523032e-11), np.float64(6.418046440501189e-11), np.float64(6.418033839618138e-11), np.float64(6.418033839454073e-11), np.float64(6.418033839742568e-11), np.float64(6.418033839463649e-11), np.float64(6.41803259770893e-11), np.float64(6.418032597676107e-11), np.float64(6.418032597741024e-11), np.float64(6.418032597662923e-11), np.float64(6.418032367204677e-11), np.float64(6.418032367208248e-11), np.float64(6.418032367202962e-11), np.float64(6.418032367205845e-11), np.float64(6.418032366726754e-11), np.float64(6.418032366727056e-11), np.float64(6.418032366726585e-11), np.float64(6.418032366726905e-11), np.float64(6.418032366726754e-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.589018638822, 4855.661711084132, 0.012435519381435901]
[5105.5890428942585, 4855.661711084132, 0.012435519381435901]
[5105.589018638822, 4855.661735829848, 0.012435519381435901]
[5105.589018638822, 4855.661711084132, 0.01243551971142224]
[5818.1902298659425, 4187.474366077638, 0.04066587315545442]
[5818.190251108292, 4187.474366077638, 0.04066587315545442]
[5818.1902298659425, 4187.4743908815235, 0.04066587315545442]
[5818.1902298659425, 4187.474366077638, 0.04066587364666456]
[5807.164852841079, 4147.480386703035, 0.041997575235356965]
[5807.164874151452, 4147.480386703035, 0.041997575235356965]
[5807.164852841079, 4147.480411453246, 0.041997575235356965]
[5807.164852841079, 4147.480386703035, 0.041997575728911545]
[5757.35968971065, 4111.453794611598, 0.04251231198035729]
[5757.359711318627, 4111.453794611598, 0.04251231198035729]
[5757.35968971065, 4111.453819307815, 0.04251231198035729]
[5757.35968971065, 4111.453794611598, 0.042512312474718955]
[5551.640563642928, 3968.2486534622326, 0.04457882939917948]
[5551.640586323431, 3968.2486534622326, 0.04457882939917948]
[5551.640563642928, 3968.2486778901684, 0.04457882939917948]
[5551.640563642928, 3968.2486534622326, 0.04457882989623189]
[2363.249376768358, 2028.034909542598, 0.08815003377117044]
[2363.2493637903126, 2028.034909542598, 0.08815003377117044]
[2363.249376768358, 2028.0349058091201, 0.08815003377117044]
[2363.249376768358, 2028.034909542598, 0.08815003409436928]
[5270.46034928257, 3795.364739489556, 0.04688423442172476]
[5270.460373065738, 3795.364739489556, 0.04688423442172476]
[5270.46034928257, 3795.3647634759877, 0.04688423442172476]
[5270.46034928257, 3795.364739489556, 0.04688423492075302]
[4949.513775661849, 3610.0316685439625, 0.049418538866642214]
[4949.513800254404, 3610.0316685439625, 0.049418538866642214]
[4949.513775661849, 3610.0316919062316, 0.049418538866642214]
[4949.513775661849, 3610.0316685439625, 0.0494185393666084]
[2440.7362736686623, 2341.042395445625, 0.07075972420068198]
[2440.7362878441095, 2341.042395445625, 0.07075972420068198]
[2440.7362736686623, 2341.0424080507905, 0.07075972420068198]
[2440.7362736686623, 2341.042395445625, 0.07075972465554832]
[4677.684454051518, 3460.7172196365955, 0.051517385561693434]
[4677.684478988293, 3460.7172196365955, 0.051517385561693434]
[4677.684454051518, 3460.7172423739903, 0.051517385561693434]
[4677.684454051518, 3460.7172196365955, 0.05151738606146314]
[4198.1387128777515, 3209.915463716623, 0.05518699198830375]
[4198.138737694841, 3209.915463716623, 0.05518699198830375]
[4198.1387128777515, 3209.9154851308285, 0.05518699198830375]
[4198.1387128777515, 3209.915463716623, 0.05518699248560598]
[4786.937523568357, 3491.189273237789, 0.05073182833636444]
[4786.937548403146, 3491.189273237789, 0.05073182833636444]
[4786.937523568357, 3491.1892961120125, 0.05073182833636444]
[4786.937523568357, 3491.189273237789, 0.05073182883631088]
[4316.412724865164, 3265.0317793454296, 0.054292006852435765]
[4316.412749797664, 3265.0317793454296, 0.054292006852435765]
[4316.412724865164, 3265.031801082158, 0.054292006852435765]
[4316.412724865164, 3265.0317793454296, 0.05429200735059023]
[4284.404410148909, 3228.845832072341, 0.05493288528636168]
[4284.404435055772, 3228.845832072341, 0.05493288528636168]
[4284.404410148909, 3228.845853599457, 0.05493288528636168]
[4284.404410148909, 3228.845832072341, 0.05493288578392239]
[4145.561159800911, 3098.522362436569, 0.057620449144262814]
[4145.561184548382, 3098.522362436569, 0.057620449144262814]
[4145.561159800911, 3098.5223831388885, 0.057620449144262814]
[4145.561159800911, 3098.522362436569, 0.05762044963842157]
[4126.632854965879, 3101.825177880715, 0.058242179053020315]
[4126.632879685501, 3101.825177880715, 0.058242179053020315]
[4126.632854965879, 3101.8251986053556, 0.058242179053020315]
[4126.632854965879, 3101.825177880715, 0.058242179546180176]
[4017.8316356185815, 3184.58345168471, 0.062491130889201156]
[4017.8316601492024, 3184.58345168471, 0.062491130889201156]
[4017.8316356185815, 3184.5834729442477, 0.062491130889201156]
[4017.8316356185815, 3184.58345168471, 0.06249113137334697]
[4030.630557947752, 3216.413561474998, 0.06286712005276221]
[4030.6305825031836, 3216.413561474998, 0.06286712005276221]
[4030.630557947752, 3216.413582928217, 0.06286712005276221]
[4030.630557947752, 3216.413561474998, 0.06286712053592235]
[4019.0353342818726, 3214.7623678999807, 0.06308862304274042]
[4019.0353588148564, 3214.7623678999807, 0.06308862304274042]
[4019.0353342818726, 3214.7623893433115, 0.06308862304274042]
[4019.0353342818726, 3214.7623678999807, 0.06308862352530523]
[4019.7015733526864, 3213.618119420504, 0.06306124378520396]
[4019.7015978869754, 3213.618119420504, 0.06306124378520396]
[4019.7015733526864, 3213.6181408569723, 0.06306124378520396]
[4019.7015733526864, 3213.618119420504, 0.06306124426784294]
[4019.465244631916, 3213.397877778115, 0.06306791802089486]
[4019.4652691657425, 3213.397877778115, 0.06306791802089486]
[4019.465244631916, 3213.3978992132616, 0.06306791802089486]
[4019.465244631916, 3213.397877778115, 0.06306791850351577]
[4019.2484811180757, 3213.320478275058, 0.06307448425482146]
[4019.2485056514774, 3213.320478275058, 0.06307448425482146]
[4019.2484811180757, 3213.3204997097396, 0.06307448425482146]
[4019.2484811180757, 3213.320478275058, 0.0630744847374246]
[4019.239555043352, 3213.3247117710193, 0.06307478481611177]
[4019.239579576736, 3213.3247117710193, 0.06307478481611177]
[4019.239555043352, 3213.3247332057267, 0.06307478481611177]
[4019.239555043352, 3213.3247117710193, 0.06307478529871409]
[4019.239555043352, 3213.3247117710193, 0.06307478481611177]

Total fitting time:
5.624126434326172 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")