Legacy #1: Temperature fit of CO2 spectrum

Quickly fit an experimental spectrum with a one-temperature model, directly from SpectrumFactory, with fit_legacy()

The method requires a fitting model. An example model is provided in radis.tools.fitting : LTEModel(). Other models can be used; such as in the multi-temperature fit example

More advanced tools for interactive fitting of multi-dimensional, multi-slabs spectra can be found in fitroom. Finally, the GPU-accelerated example shows how to obtain real-time interactive spectra.

from radis import SpectrumFactory, load_spec

Here we get an experimental spectrum from RADIS test cases. Use your own instead.

from radis.test.utils import getTestFile, setup_test_line_databases

setup_test_line_databases()

# fit range
wlmin = 4167
wlmax = 4180

s_exp = (
    load_spec(getTestFile("CO2_measured_spectrum_4-5um.spec"))
    .crop(wlmin, wlmax, "nm")
    .normalize()
    .sort()
    .offset(-0.2, "nm")
)

Customize the LTEModel() for our case: we add a slit (non fittable parameter) and normalize it

from radis.tools.fitting import LTEModel


def LTEModel_withslitnorm(factory, fit_parameters, fixed_parameters):
    s = LTEModel(factory, fit_parameters, fixed_parameters)
    # we could also have added a fittable parameter, such as an offset,
    # or made the slit width a fittable parameter.
    # ... any parameter in model_input will be fitted.
    s.offset(0, "nm")  # or we could have used a fittable parameter below :
    # s.offset(model_input["offset"], 'nm')

    # Alternative: with a wavelength offset
    # WARNING: very sensitive parameter
    # copy_parameters = fit_parameters.copy()
    # offset = copy_parameters["offset"]
    # copy_parameters.pop("offset")
    # s = LTEModel(factory, copy_parameters, fixed_parameters)
    # s.offset(offset, 'nm')

    # Comment:
    # We could also have made the slit width a fittable parameter.
    # ... any parameter in model_input will be fitted.
    # Here we simply employ a fixed slit.
    s.apply_slit(1.4, "nm")
    return s.take("radiance").normalize()

using fit_legacy()

import astropy.units as u

sf = SpectrumFactory(
    wlmin * u.nm,
    wlmax * u.nm,
    molecule="CO2",
    wstep=0.001,  # cm-1
    pressure=1 * 1e-3,  # bar
    cutoff=1e-25,
    isotope="1,2",
    path_length=10,  # cm-1
    mole_fraction=1,
    truncation=1,  # cm-1
    verbose=0,
)
sf.warnings["MissingSelfBroadeningWarning"] = "ignore"
sf.warnings["HighTemperatureWarning"] = "ignore"
sf.load_databank(
    "HITRAN-CO2-TEST"
)  # see 'fetch_databank' below for a more general application
# sf.fetch_databank("hitemp") #use "hitemp" or another database

s_best, best = sf.fit_legacy(
    s_exp.take("radiance"),
    model=LTEModel_withslitnorm,
    fit_parameters={
        "Tgas": 1450,
        # "offset": 0
    },
    bounds={
        "Tgas": [300, 2000],
        # "offset": [-1, 1],
    },
    plot=True,
    solver_options={
        "maxiter": 15,  # 👈 increase to let the fit converge
        "ftol": 1e-15,
        # "gtol": 1e-10,
        # "eps":1e-5
    },
    verbose=2,
)
# plot_diff(s_exp, s_best)  # , show_ruler=True)

• 
• 

/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/develop/radis/misc/log.py:51: UserWarning: Reference databank (2391.79-2399.14cm-1) has 0 lines in range (2391.69-2399.15cm-1) for isotope 2. Change your range or isotope options
  warn(msg)
0.04s - Loaded database
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/develop/radis/misc/warning.py:434: LinestrengthCutoffWarning: Estimated error after discarding lines is large: 0.07%. Consider reducing cutoff
  warnings.warn(WarningType(message))

Now starting the fitting process:
---------------------------------

/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/develop/radis/tools/fitting.py:480: OptimizeWarning: Unknown solver options: disp
  best = minimize(
Tgas=1150.0, Residual: 0.0105 🏆
Tgas=1170.0, Residual: 0.0098 🏆
Tgas=1150.0, Residual: 0.0105
Tgas=1170.0, Residual: 0.0098 🏆
Tgas=2000.0, Residual: 0.0128
Tgas=1980.0, Residual: 0.0128
Tgas=1396.0, Residual: 0.0034 🏆
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1396.0, Residual: 0.0034
Tgas=1416.0, Residual: 0.0031 🏆
Tgas=1397.0, Residual: 0.0034
Tgas=1417.0, Residual: 0.0031 🏆
Tgas=1402.0, Residual: 0.0033
Tgas=1422.0, Residual: 0.0030 🏆
Tgas=1420.0, Residual: 0.0030
Tgas=1440.0, Residual: 0.0028 🏆
Tgas=1459.0, Residual: 0.0027 🏆
Tgas=1479.0, Residual: 0.0028
Init ['Tgas'] = [1150.]['']
Final ['Tgas'] = [1459.]['']
  message: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
  success: True
   status: 0
      fun: 0.0027304260508692734
        x: [ 1.459e+03]
      nit: 4
      jac: [ 3.330e-06]
     nfev: 32
     njev: 16
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
Best ['Tgas'] = [1458.71148816][''] reached at iteration 32/32