Chi-square minimization for multiple observed spectra
Forward modeling of multiple spectra for the same object by minimizing the chi-square statistics between the data and atmospheric models.
[1]:
import seda # import the seda package
import importlib
import numpy as np
import matplotlib.pyplot as plt
import os
from matplotlib.ticker import MultipleLocator, FormatStrFormatter, AutoMinorLocator, StrMethodFormatter, NullFormatter
from astropy.io import fits, ascii
SEDA v0.5.6.dev3 package imported
Read the spectra
As an example here, let’s read the near-infrared IRTF/SpeX and the mid-infrared Spitzer/IRS spectra for the T8 (~750 K) brown dwarf 2MASS J04151954-0935066 in Burgasser et al. (2004) and Suárez & Metchev (2022), respectively.
Read SpeX spectrum:
[2]:
# path to the seda package
path_seda = os.path.dirname(os.path.dirname(seda.__file__))
# SpeX spectrum
SpeX_name = path_seda+'/docs/notebooks/data/0415-0935_IRTF_SpeX.dat'
SpeX = ascii.read(SpeX_name)
wl_SpeX = SpeX['wl(um)'] # um
flux_SpeX = SpeX['flux(erg/s/cm2/A)'] # erg/s/cm2/A
eflux_SpeX = SpeX['eflux(erg/s/cm2/A)'] # erg/s/cm2/A
Read IRS spectrum:
[3]:
IRS_name = path_seda+'/docs/notebooks/data/0415-0935_IRS_spectrum.dat'
IRS = ascii.read(IRS_name)
wl_IRS = IRS['wl(um)'] # in um
flux_IRS = IRS['flux(Jy)'] # in Jy
eflux_IRS = IRS['eflux(Jy)'] # in Jy
# convert IRS fluxes from Jy to erg/s/cm2/A
out_convert_flux = seda.synthetic_photometry.convert_flux(wl=wl_IRS, flux=flux_IRS, eflux=eflux_IRS,
unit_in='Jy', unit_out='erg/s/cm2/A')
flux_IRS = out_convert_flux['flux_out'] # in erg/s/cm2/A
eflux_IRS = out_convert_flux['eflux_out'] # in erg/s/cm2/A
Plot SED to check everything looks okay:
[4]:
fig, ax = plt.subplots()
plt.plot(wl_SpeX, flux_SpeX, label='SpeX spectrum')
plt.plot(wl_IRS, flux_IRS, label='IRS spectrum')
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.yaxis.set_minor_locator(AutoMinorLocator())
plt.xscale('log')
plt.yscale('log')
ax.xaxis.set_major_formatter(StrMethodFormatter('{x:.0f}'))
ax.grid(True, which='both', color='gainsboro', linewidth=0.5, alpha=0.5)
ax.legend()
plt.xlabel(r'$\lambda\ (\mu$m)', size=12)
plt.ylabel(r'$F_\lambda\ ($erg s$^{-1}$ cm$^{-2}$ $\AA^{-1}$)', size=12)
plt.show()
Load input data
Look at the input parameters here.
For any SEDA function, we can also see the function description directly on the notebook with the command help(), e.g.:
help(seda.input_parameters.InputData)
Define lists with the input spectra:
[5]:
# wavelenghts
wl_spectra = [wl_SpeX, wl_IRS] # in um
# fluxes
flux_spectra = [flux_SpeX, flux_IRS] # in erg/s/cm2/A
# flux uncertainties
eflux_spectra = [eflux_SpeX, eflux_IRS] # in erg/s/cm2/A
# specify flux units
flux_unit = 'erg/s/cm2/A'
# resolution of each input spectrum (used to convolve the model spectra)
res = [100, 100]
# distance to the target (optional and used to derive a radius)
distance = 5.71 # pc (parallax=175.2+-1.7; Dupuy-Liu2012)
edistance = 0.06 # pc
# load all the input data parameters
my_data = seda.input_parameters.InputData(wl_spectra=wl_spectra, flux_spectra=flux_spectra,
eflux_spectra=eflux_spectra, flux_unit=flux_unit, res=res,
distance=distance, edistance=edistance)
Input data loaded successfully:
2 spectra
Download (if not yet) the atmospheric models you want to use.
Use the commands below to see the available atmospheric models, the links to download them, and other relevant information from models. You can read more about the model here.
Also consider this tutorial to explore the free parameters in the models and their coverage.
[6]:
# available atmospheric models
seda.models.Models().available_models
[6]:
['BT-Settl',
'ATMO2020',
'Sonora_Elf_Owl',
'SM08',
'Sonora_Bobcat',
'Sonora_Diamondback',
'Sonora_Cholla',
'LB23']
[7]:
# some parameters of interest from a selected model
model = 'Sonora_Elf_Owl'
print(seda.models.Models(model).ref) # reference
print(seda.models.Models(model).ADS) # link to paper
print(seda.models.Models(model).download) # link to download the models
seda.models.Models(model).params_unique # coverage of free parameters in the grid
Mukherjee et al. (2024)
https://ui.adsabs.harvard.edu/abs/2024ApJ...963...73M/abstract
['https://zenodo.org/records/10385987', 'https://zenodo.org/records/10385821', 'https://zenodo.org/records/10381250']
[7]:
{'Teff': array([ 275., 300., 325., 350., 375., 400., 425., 450., 475.,
500., 525., 550., 575., 600., 650., 700., 750., 800.,
850., 900., 950., 1000., 1100., 1200., 1300., 1400., 1500.,
1600., 1700., 1800., 1900., 2000., 2100., 2200., 2300., 2400.]),
'logg': array([3. , 3.25, 3.5 , 3.75, 4. , 4.25, 4.5 , 4.75, 5. , 5.25, 5.5 ]),
'logKzz': array([2., 4., 7., 8., 9.]),
'Z': array([-1. , -0.5, 0. , 0.5, 0.7, 1. ]),
'CtoO': array([0.5, 1. , 1.5, 2. , 2.5])}
Load model grid options
Look at the input parameters here.
[8]:
# select the atmospheric models of interest
model = 'Sonora_Elf_Owl'
# path to the directory or directories containing the model spectra
# (update it to your own path)
my_path = '/home/gsuarez/TRABAJO/MODELS/atmosphere_models/Sonora_Elf_Owl/spectra/'
model_dir = [my_path+'output_700.0_800.0/',
my_path+'output_850.0_950.0/',
]
# set parameter ranges to select a grid subset
# when a free parameter range is not specified, the whole grid range will be considered
params_ranges = {
'Teff': [700, 900], # Teff range
'logg': [4.0, 5.0] # logg range
}
# load model options
my_model = seda.input_parameters.ModelOptions(model=model, model_dir=model_dir,
params_ranges=params_ranges)
Model options loaded successfully
Tip: If you plan to model several spectra from the same instrument (same resolution), you can save the convolved model spectra to reuse them and do subsequent fits much faster. For this, set the parameter path_save_spectra_conv in seda.ModelOptions above to a folder path where you want to store the convolved spectra. Once the spectra are stored, the next time you run the code just replace model_dir by the path you used in path_save_spectra_conv and set
skip_convolution=True to avoid the model convolution. This was implemented thanks to this issue.
Load chi-square fit options
Look at the input parameters here.
Consider the default full wavelength range of each input spectrum for the fits. Otherwise, we can use the parameter fit_wl_range to set different fit ranges.
[11]:
# choose a filename (optional) to save the sampling results as a pickle file
# it is convenient to set a non-default name when running
# the code several times in the same folder to avoid overwriting results
chi2_pickle_file = f'{model}_chi2_minimization_multiple_spectra.pickle'
chi2_table_file = f'{model}_chi2_minimization_multiple_spectra.dat'
# load chi-square fit options
my_chi2 = seda.input_parameters.Chi2Options(my_data=my_data, my_model=my_model,
chi2_pickle_file=chi2_pickle_file,
chi2_table_file=chi2_table_file)
3000 model spectra selected with:
Teff range = [700, 900]
logg range = [4.0, 5.0]
Chi-square fit options loaded successfully
elapsed time: 14.1 min
Run chi-square fit
[14]:
out_chi2 = seda.chi2_fit.chi2(my_chi2=my_chi2)
Running chi-square fitting...
chi-square minimization results saved successfully
Chi-square fit ran successfully
elapsed time: 7.0 s
Plot results
SED with the best three fits.
The pickle file generated by seda.chi2_fit and stored with the name my_chi2.chi2_pickle_file is the input file to make plots. We can provide the name by either using my_chi2.chi2_pickle_file (if my_chi2 is in memory) or just typing it.
Note that Elf Owl models do not cover wavelenghts longer than 15 microns.
[5]:
# plot wavelengths in logarithmic scale
chi2_pickle_file = 'Sonora_Elf_Owl_chi2_minimization_multiple_spectra.pickle'
fig, ax = seda.plots.plot_chi2_fit(chi2_pickle_file, N_best_fits=3, xlog=True,
yrange=[1e-20, 1e-15], ori_res=True)
Zoom in on the SpeX spectrum:
[6]:
# plot fluxes in linear scale
fig, ax = seda.plots.plot_chi2_fit(chi2_pickle_file, N_best_fits=3,
xrange=[wl_SpeX.min(), wl_SpeX.max()],
ylog=False, save=False)
Zoom in on the IRS spectrum:
[7]:
# plot fluxes in linear scale
fig, ax = seda.plots.plot_chi2_fit(chi2_pickle_file, N_best_fits=3,
xrange=[wl_IRS.min(), wl_IRS.max()],
ylog=False)
Reduced chi-square as a function of wavelength for the best three fits.
[8]:
fig, ax = seda.plots.plot_chi2_red(chi2_pickle_file, N_best_fits=3, xlog=True)
