Free retrieval with petitRADTRANSΒΆ

This tutorial contains only a code for running a free parameter retrieval on the data of ROXs 42 Bb. It is required to manually install petitRADTRANS and recommended to run the retrieval on a cluster, in particular with scattering=True. More details will follow.

Make sure to comment anything Database related out when running run_multinest because writing to the HDF5 database is not possible with multiprocessing. Therefore, first add the companion data to the HDF5 database, then run AtmosphericRetrieval and run_multinest on a cluster, and then create all the plots while commenting out the retrieval part.

[1]:
import species

species.SpeciesInit()

database = species.Database()

database.add_object(object_name='ROXs 42 Bb',
                    distance=(135., 0.),
                    app_mag={'Keck/NIRC2.J': (16.91, 0.11),
                             'Keck/NIRC2.H': (15.8, 0.05),
                             'Keck/NIRC2.Ks': (15.01, 0.06),
                             'Keck/NIRC2.Lp': (13.97, 0.06),
                             'Keck/NIRC2.NB405': (13.90, 0.12),
                             'Keck/NIRC2.Ms': (14.01, 0.23)},
                    spectrum={'NIFS': ('data/roxs42bb_nifs_j_calib.dat', None, 5000.),
                              'OSIRIS': ('data/roxs42bb_osiris_h_calib.dat', None, 3000.),
                              'SINFONI': ('data/roxs42bb_sinfoni_k_calib.dat', None, 1500.)},
                    deredden=None)

retrieve = species.AtmosphericRetrieval(object_name='ROXs 42 Bb',
                                        line_species=['CO_all_iso', 'H2O', 'CH4', 'NH3', 'CO2', 'Na', 'K', 'TiO', 'VO', 'FeH', 'H2S'],
                                        cloud_species=['MgSiO3(c)_cd', 'Fe(c)_cd', 'Al2O3(c)_cd'],
                                        scattering=True,
                                        output_folder='multinest',
                                        wavel_range=(1.1, 2.46),
                                        inc_spec=['NIFS', 'OSIRIS', 'SINFONI'],
                                        inc_phot=False,
                                        pressure_grid='smaller',
                                        weights=None)

retrieve.run_multinest(bounds={'logg': (2., 5.5),
                               'c_o_ratio': (0.1, 1.),
                               'metallicity': (-3., 3.),
                               'radius': (1., 5.),
                               'OSIRIS': ((0.8, 1.2), None, None),
                               'fsed': (0., 10.),
                               'log_kzz': (4., 12.),
                               'log_tau_cloud': (-2., 1.),
                               'fe_mgsio3_ratio': (-2., 2.),
                               'al2o3_mgsio3_ratio': (-2., 2.),
                               'ism_ext': (1.7, 1.7)},
                       chemistry='equilibrium',
                       quenching=None,
                       pt_profile='molliere',
                       n_live_points=1000,
                       resume=True,
                       plotting=False,
                       pt_smooth=0.)

database.add_retrieval(tag='roxs42bb',
                       output_folder='multinest',
                       inc_teff=False)

database.get_retrieval_teff(tag='roxs42bb',
                            random=30)

species.plot_posterior(tag='roxs42bb',
                       offset=(-0.3, -0.35),
                       vmr=True,
                       inc_mass=True,
                       inc_pt_param=False,
                       output='plot/posterior.pdf')

samples, radtrans = database.get_retrieval_spectra(tag='roxs42bb',
                                                   random=30,
                                                   wavel_range=(0.5, 6.),
                                                   spec_res=2000.)

species.plot_pt_profile(tag='roxs42bb',
                        random=100,
                        xlim=(0., 6000.),
                        offset=(-0.07, -0.14),
                        output='plot/pt_profile_grains.pdf',
                        radtrans=radtrans,
                        extra_axis='grains')

species.plot_opacities(tag='roxs42bb',
                       offset=(-0.1, -0.14),
                       output='plot/opacities.pdf',
                       radtrans=radtrans)

species.plot_clouds(tag='roxs42bb',
                    offset=(-0.1, -0.15),
                    output='plot/clouds.pdf',
                    radtrans=radtrans,
                    composition='MgSiO3')

best = database.get_probable_sample(tag='roxs42bb')

objectbox = database.get_object('ROXs 42 Bb',
                                inc_phot=False)

objectbox = species.update_spectra(objectbox, best)

residuals = species.get_residuals(datatype='model',
                                  spectrum='petitradtrans',
                                  parameters=best,
                                  objectbox=objectbox,
                                  inc_phot=False,
                                  inc_spec=True,
                                  radtrans=radtrans)

modelbox = radtrans.get_model(model_param=best,
                              spec_res=2000.,
                              plot_contribution='plot/contribution.pdf')

no_clouds = best.copy()
no_clouds['log_tau_cloud'] = -100.
model_no_clouds = radtrans.get_model(no_clouds)

species.plot_spectrum(boxes=[samples, modelbox, model_no_clouds, objectbox],
                      filters=None,
                      plot_kwargs=[{'ls': '-', 'lw': 0.1, 'color': 'gray'},
                                   {'ls': '-', 'lw': 0.5, 'color': 'black'},
                                   {'ls': '--', 'lw': 0.3, 'color': 'black'},
                                   {'NIFS': {'marker': 'o', 'ms': 2., 'color': 'tab:green', 'ls': 'none', 'alpha': 0.2, 'mew': 0., 'label': 'NIFS'},
                                    'OSIRIS': {'marker': 'o', 'ms': 2., 'color': 'tab:blue', 'ls': 'none', 'alpha': 0.2, 'mew': 0., 'label': 'OSIRIS'},
                                    'SINFONI': {'marker': 'o', 'ms': 2., 'color': 'tab:orange', 'ls': 'none', 'alpha': 0.2, 'mew': 0., 'label': 'SINFONI'}}],
                      residuals=residuals,
                      xlim=(1.1, 2.5),
                      ylim=(0.15e-16, 1.15e-15),
                      ylim_res=(-5., 5.),
                      scale=('linear', 'linear'),
                      offset=(-0.6, -0.05),
                      figsize=(8, 2.5),
                      legend=[{'loc': 'upper right', 'fontsize': 8.}, {'loc': 'lower left', 'fontsize': 8.}],
                      output='plot/spectrum.pdf')

# Plot the SED over broader wavelength range

modelbox = radtrans.get_model(model_param=best,
                              spec_res=200.,
                              plot_contribution='plot/contribution.pdf')

objectbox = database.get_object('ROXs 42 Bb',
                                inc_phot=True)

for item in ['NIFS', 'OSIRIS', 'SINFONI']:
    objectbox.spectrum[item] = (objectbox.spectrum[item][0][::50, ],
                                objectbox.spectrum[item][1],
                                objectbox.spectrum[item][2],
                                objectbox.spectrum[item][3])

residuals = species.get_residuals(datatype='model',
                                  spectrum='petitradtrans',
                                  parameters=best,
                                  objectbox=objectbox,
                                  inc_phot=True,
                                  inc_spec=True,
                                  radtrans=radtrans)

synphot = species.multi_photometry(datatype='model',
                                   spectrum='petitradtrans',
                                   filters=objectbox.filters,
                                   parameters=best,
                                   radtrans=radtrans)

object_kwargs = {'NIFS': {'marker': 'o', 'ms': 2.5, 'color': 'tab:green', 'ls': 'none', 'alpha': 0.5, 'mew': 0., 'label': 'Gemini/NIFS'},
                 'OSIRIS': {'marker': 'o', 'ms': 2.5, 'color': 'tab:blue', 'ls': 'none', 'alpha': 0.5, 'mew': 0., 'label': 'Keck/OSIRIS'},
                 'SINFONI': {'marker': 'o', 'ms': 2.5, 'color': 'tab:orange', 'ls': 'none', 'alpha': 0.5, 'mew': 0., 'label': 'VLT/SINFONI'}}

for item in objectbox.magnitude:
    object_kwargs[item] = {'marker': 's', 'ms': 3.5, 'color': 'black'}

species.plot_spectrum(boxes=[samples, modelbox, objectbox, synphot],
                      filters=None,
                      plot_kwargs=[{'ls': '-', 'lw': 0.1, 'color': 'gray'},
                                   {'ls': '-', 'lw': 0.5, 'color': 'black'},
                                   object_kwargs,
                                   None],
                      residuals=residuals,
                      xlim=(0.8, 5.),
                      ylim=(0.15e-15, 1.9e-15),
                      ylim_res=(-7., 7.),
                      scale=('linear', 'linear'),
                      offset=(-0.6, -0.04),
                      figsize=(8, 2.5),
                      legend=[None, {'loc': 'upper right', 'fontsize': 10.}],
                      quantity='flux',
                      output='plot/full_sed.pdf')
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_1007/3466150118.py in <module>
----> 1 import species
      2
      3 species.SpeciesInit()
      4
      5 database = species.Database()

ModuleNotFoundError: No module named 'species'