Multiple species¶
This example demonstrates how to combine multiple species in one go for the simulation.
To do that, it simply concatenates the databases for two or more species.
To ensure we retain the information which line corresponds with which species, we add a species column with this information.
For simplicity of this interactive example, the spectra are simulated based on the linestrength only, inherently assuming an equal amount for each selected species.
Furthermore, it assumes equilibration between rotational or vibrational states, such that a single $T_{rot}$ and $T_{vib}$ can be used to describe the spectrum.
In [1]:
Copied!
from pathlib import Path
import Moose
import lmfit
import pandas as pd
import numpy as np
import panel as pn
import panel.widgets as pnw
from importlib import resources
import plotly.graph_objects as go
pn.extension('plotly')
species = [p.stem for p in Path(resources.files("Moose").joinpath("data")).glob('*.db')]
params = lmfit.create_params(**Moose.default_params)
params['A'].vary = False
params['b'].vary = False
params['wl_pad'].value = 200
params['T_vib'].value= 4000
def load_species_db(name, *args):
df = Moose.query_DB(name,*args)
df['Species'] = name
return df.sort_values('air_wavelength')
dbs = pd.concat([load_species_db(n,(300,500)) for n in ['N2CB', 'C2_swan']])
x = np.linspace(300,500,2000)
y = Moose.model_for_fit(x,**{p:params[p].value for p in params}, sim_db=dbs)
from pathlib import Path
import Moose
import lmfit
import pandas as pd
import numpy as np
import panel as pn
import panel.widgets as pnw
from importlib import resources
import plotly.graph_objects as go
pn.extension('plotly')
species = [p.stem for p in Path(resources.files("Moose").joinpath("data")).glob('*.db')]
params = lmfit.create_params(**Moose.default_params)
params['A'].vary = False
params['b'].vary = False
params['wl_pad'].value = 200
params['T_vib'].value= 4000
def load_species_db(name, *args):
df = Moose.query_DB(name,*args)
df['Species'] = name
return df.sort_values('air_wavelength')
dbs = pd.concat([load_species_db(n,(300,500)) for n in ['N2CB', 'C2_swan']])
x = np.linspace(300,500,2000)
y = Moose.model_for_fit(x,**{p:params[p].value for p in params}, sim_db=dbs)
In [2]:
Copied!
fig = go.Figure()
fig.add_scatter(x=x, y=y, name='Spectrum')
fig.update_xaxes(title_text='λ (nm)')
fig.update_yaxes(title_text='I (a.u.)')
select_db = pnw.MultiChoice(value=['N2CB','C2_swan',],options = species, name='Select transitions', max_items=4)
select_range = pnw.RangeSlider(name='λ (nm)', start=200, end=700, value=(300,500))
pop_fracs = [pnw.FloatInput(name=f'Pop. species {i}', value=1,start=0,end=100, step=0.1) for i in range(1,5)]
status_msg = pnw.TextInput(name='Log', disabled=True)
input_widgets = {p:pnw.FloatSlider(name=p, value=params[p].value, start=params[p].min, end=params[p].max, step=params[p].value/10) for p in params if params[p].vary is True}
def cb_update_spectrum(event):
x = np.linspace(*select_range.value,2000)
if event.new != []:
dbs = pd.concat([load_species_db(n, select_range.value) for n in select_db.value])
y = np.zeros(x.shape[0])
for i, (n,g) in enumerate(dbs.groupby('Species')):
y+=Moose.model_for_fit(x,**{k:v.value for k,v in input_widgets.items()}, **{p:params[p].value for p in params if p not in input_widgets.keys()}, sim_db=g)*pop_fracs[i].value
with fig.batch_update():
fig.data[0].x = x
fig.data[0].y=y
else:
fig.data[0].x = x
fig.data[0].y = np.zeros(x.shape[0])
status_msg.value = "Please select at least one species"
for wid in list(input_widgets.values())+[select_range,select_db]+pop_fracs:
if hasattr(wid,'value_throttled'):
to_watch = 'value_throttled'
else:
to_watch = 'value'
wid.param.watch(cb_update_spectrum,to_watch)
layout = pn.Column(pn.pane.Plotly(fig, sizing_mode='stretch_width'),pn.FlexBox(select_db, select_range, status_msg), pn.FlexBox(*input_widgets.values()), pn.Row(*pop_fracs))
layout
fig = go.Figure()
fig.add_scatter(x=x, y=y, name='Spectrum')
fig.update_xaxes(title_text='λ (nm)')
fig.update_yaxes(title_text='I (a.u.)')
select_db = pnw.MultiChoice(value=['N2CB','C2_swan',],options = species, name='Select transitions', max_items=4)
select_range = pnw.RangeSlider(name='λ (nm)', start=200, end=700, value=(300,500))
pop_fracs = [pnw.FloatInput(name=f'Pop. species {i}', value=1,start=0,end=100, step=0.1) for i in range(1,5)]
status_msg = pnw.TextInput(name='Log', disabled=True)
input_widgets = {p:pnw.FloatSlider(name=p, value=params[p].value, start=params[p].min, end=params[p].max, step=params[p].value/10) for p in params if params[p].vary is True}
def cb_update_spectrum(event):
x = np.linspace(*select_range.value,2000)
if event.new != []:
dbs = pd.concat([load_species_db(n, select_range.value) for n in select_db.value])
y = np.zeros(x.shape[0])
for i, (n,g) in enumerate(dbs.groupby('Species')):
y+=Moose.model_for_fit(x,**{k:v.value for k,v in input_widgets.items()}, **{p:params[p].value for p in params if p not in input_widgets.keys()}, sim_db=g)*pop_fracs[i].value
with fig.batch_update():
fig.data[0].x = x
fig.data[0].y=y
else:
fig.data[0].x = x
fig.data[0].y = np.zeros(x.shape[0])
status_msg.value = "Please select at least one species"
for wid in list(input_widgets.values())+[select_range,select_db]+pop_fracs:
if hasattr(wid,'value_throttled'):
to_watch = 'value_throttled'
else:
to_watch = 'value'
wid.param.watch(cb_update_spectrum,to_watch)
layout = pn.Column(pn.pane.Plotly(fig, sizing_mode='stretch_width'),pn.FlexBox(select_db, select_range, status_msg), pn.FlexBox(*input_widgets.values()), pn.Row(*pop_fracs))
layout
Out[2]:
Model and fit the combined spectrum of two species¶
Alternatively, the provided sample file contains emission of the CN violet system and the N$_2$ second positive system.
In the case of the former there are multiple overlapping vibrational bands, but the latter only contains one vibrational band, for the (0,2) transition.
Thus to model the spectrum we allow for a vibrational and rotational temperature for CN, while N$_2$ only is assigned a rotational temperature.
In [3]:
Copied!
def model_two_species(x:np.array, mu:float=0,sigma:float=0.01,gamma:float=0.01,T_rot_1:float=300,T_rot_2:float=300,T_vib:float=300, ratio:float=1, A:float=1,b:float=0,resolution:int=100, wl_pad:float=10,sim_db1=None,sim_db2=None):
"""A model function that combines two line-by-line dataframes for two species.
The respective contributions between the two species are scaled according to the `ratio` argument"""
sticks = np.concatenate([Moose.create_stick_spectrum(300,T_rot_1,sim_db1), Moose.create_stick_spectrum(T_vib,T_rot_2,sim_db2)*np.array([1,ratio])])
equid = Moose.equidistant_mesh(sticks[sticks[:,0].argsort()],wl_pad,resolution)
v = Moose.apply_voigt(equid,sigma, gamma)
matched = Moose.match_spectra((x-mu).reshape(-1,1), v)
matched[:,1] = (matched[:,1]-matched[:,1].min())/(matched[:,1].max()-matched[:,1].min())
return A*matched[:,1]+b
db1 = Moose.query_DB('N2CB', (350,400))
db2 = Moose.query_DB('CNBX', (350,400))
x=np.linspace(370,390,1000)
y=model_two_species(x,0,0.1,0.1,300,300,300,1,sim_db1=db1,sim_db2=db2)
def model_two_species(x:np.array, mu:float=0,sigma:float=0.01,gamma:float=0.01,T_rot_1:float=300,T_rot_2:float=300,T_vib:float=300, ratio:float=1, A:float=1,b:float=0,resolution:int=100, wl_pad:float=10,sim_db1=None,sim_db2=None):
"""A model function that combines two line-by-line dataframes for two species.
The respective contributions between the two species are scaled according to the `ratio` argument"""
sticks = np.concatenate([Moose.create_stick_spectrum(300,T_rot_1,sim_db1), Moose.create_stick_spectrum(T_vib,T_rot_2,sim_db2)*np.array([1,ratio])])
equid = Moose.equidistant_mesh(sticks[sticks[:,0].argsort()],wl_pad,resolution)
v = Moose.apply_voigt(equid,sigma, gamma)
matched = Moose.match_spectra((x-mu).reshape(-1,1), v)
matched[:,1] = (matched[:,1]-matched[:,1].min())/(matched[:,1].max()-matched[:,1].min())
return A*matched[:,1]+b
db1 = Moose.query_DB('N2CB', (350,400))
db2 = Moose.query_DB('CNBX', (350,400))
x=np.linspace(370,390,1000)
y=model_two_species(x,0,0.1,0.1,300,300,300,1,sim_db1=db1,sim_db2=db2)
In [4]:
Copied!
df = pd.read_csv("CN/370-400nm_9ns_calib.txt").dropna()
model = lmfit.Model(model_two_species, sim_db1=db1, sim_db2=db2)
fit_params = lmfit.create_params(
**{k: v for k, v in Moose.default_params.items() if k not in ["T_rot"]}
)
fit_params.add("T_rot_1", value=1000, min=300, max=1e4)
fit_params.add("T_rot_2", value=6000, min=300, max=1e4)
fit_params.add("ratio", value=1, min=0, max=100)
fig_fit = go.Figure()
fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum")
fig_fit.add_scatter(
x=df["Wavelength"],
y=model.eval(x=df["Wavelength"].values, params=fit_params),
name="Fit",
)
fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual")
fig_fit.update_xaxes(title_text="λ (nm)")
fig_fit.update_yaxes(title_text="I (a.u.)")
btn_fit = pnw.Button(name="Fit", button_type="primary")
widgets = {
p: pnw.FloatInput(name=p, value=fit_params[p].value)
for p in fit_params
if p not in ["A", "b"]
}
def iter_cb(params, iter, resid, *args, **kwargs):
for p in params:
if p not in ["A", "b"]:
widgets[p].value = params[p].value
if iter % 10 == 0:
with fig_fit.batch_update():
fig_fit.data[2].y = resid
# fig_fit.data[1].y = model_two_species(**params.valuesdict(),**kwargs, sim_db1=db1,sim_db)
fig_fit.data[1].y = model.eval(params=params, x=kwargs.get("x"))
def cb_fit(event):
btn_fit.button_type = "warning"
result = model.fit(
df["I"].values, x=df["Wavelength"].values, params=fit_params, iter_cb=iter_cb
)
btn_fit.button_type = "primary"
iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values)
btn_fit.on_click(cb_fit)
layout = pn.Column(
pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"),
btn_fit,
pn.FlexBox(*widgets.values()),
)
layout
df = pd.read_csv("CN/370-400nm_9ns_calib.txt").dropna()
model = lmfit.Model(model_two_species, sim_db1=db1, sim_db2=db2)
fit_params = lmfit.create_params(
**{k: v for k, v in Moose.default_params.items() if k not in ["T_rot"]}
)
fit_params.add("T_rot_1", value=1000, min=300, max=1e4)
fit_params.add("T_rot_2", value=6000, min=300, max=1e4)
fit_params.add("ratio", value=1, min=0, max=100)
fig_fit = go.Figure()
fig_fit.add_scatter(x=df["Wavelength"], y=df["I"], name="Spectrum")
fig_fit.add_scatter(
x=df["Wavelength"],
y=model.eval(x=df["Wavelength"].values, params=fit_params),
name="Fit",
)
fig_fit.add_scatter(x=df["Wavelength"], y=np.zeros(df.shape[0]), name="Residual")
fig_fit.update_xaxes(title_text="λ (nm)")
fig_fit.update_yaxes(title_text="I (a.u.)")
btn_fit = pnw.Button(name="Fit", button_type="primary")
widgets = {
p: pnw.FloatInput(name=p, value=fit_params[p].value)
for p in fit_params
if p not in ["A", "b"]
}
def iter_cb(params, iter, resid, *args, **kwargs):
for p in params:
if p not in ["A", "b"]:
widgets[p].value = params[p].value
if iter % 10 == 0:
with fig_fit.batch_update():
fig_fit.data[2].y = resid
# fig_fit.data[1].y = model_two_species(**params.valuesdict(),**kwargs, sim_db1=db1,sim_db)
fig_fit.data[1].y = model.eval(params=params, x=kwargs.get("x"))
def cb_fit(event):
btn_fit.button_type = "warning"
result = model.fit(
df["I"].values, x=df["Wavelength"].values, params=fit_params, iter_cb=iter_cb
)
btn_fit.button_type = "primary"
iter_cb(result.params, -10, result.residual, x=df["Wavelength"].values)
btn_fit.on_click(cb_fit)
layout = pn.Column(
pn.pane.Plotly(fig_fit, sizing_mode="stretch_width"),
btn_fit,
pn.FlexBox(*widgets.values()),
)
layout
Out[4]:
In [ ]:
Copied!