Note

Download full notebook here

Querying Opacities

PICASO currently includes several different flavors of opacity database on Zenodo (incl different Rs, molecules, wavelength ranges).

We chose those use sqlite3 for our database because of it’s 1) user-friendliness, 2) speed, 3) scalability, 4) compatibility with parallel programming. In 2019, we tried out various other methods as well– json, hdf5, ascii, sqlalchemy– but sqlite3 was truly better for this specific problem. Having revisited this problem recently hdf5 has now outpaced sqlite3 in terms of storage and performance. Therefore in future PICASO v5 we will move away from sqlite3 and toward hdf5 which has a much larger user base now.

In this tutorial you wil learn how to quickly grab any opacity data

from picaso import justplotit as jpi
from picaso import justdoit as jdi
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import picaso.opacity_factory as opa
#plotting
jpi.output_notebook()

WARNING: Failed to load Vega spectrum from /data/reference_data/picaso/ref4/stellar_grids/calspec/alpha_lyr_stis_011.fits; Functionality involving Vega will be severely limited: FileNotFoundError(2, 'No such file or directory') [stsynphot.spectrum]

Loading BokehJS ...

General Sqlite3 and how our database is structured

• 3 Tables: header, continuum, and molecular

• header: contains units and the wavenumber grid

• continuum: contains a grid continuum opacity and is temperature dependent

• molecular: contains all molecular opacity and is pressure-temperature dependent.

How to Query the Database

opanxn = jdi.opannection()
#this is where your opacity file should be located if you've set your environments correctly
db_filename = opanxn.db_filename

What version db are you using ?

Note this only works with PICASO Opacities v3 on Zenodo and above

opa.get_metadata_item(db_filename, 'version'),opa.get_metadata_item(db_filename, 'zenodo')

('default_3.3', '10.5281/zenodo.14861730')

What Species, Pressures and Temperatures are Available

molecules, pt_pairs = opa.molecular_avail(db_filename)

print(molecules)

['C2H2', 'CH4', 'CO', 'CO2', 'H2O', 'H2S', 'K', 'Na', 'OCS', 'SO2']

pt_pairs[0:10] #full list of all the pt pairs. the first value is just an index for bookkeeping

[(1, 1e-06, 75.0),
 (2, 3e-06, 75.0),
 (3, 1e-05, 75.0),
 (4, 3e-05, 75.0),
 (5, 0.0001, 75.0),
 (6, 0.0003, 75.0),
 (7, 0.001, 75.0),
 (8, 0.003, 75.0),
 (9, 0.01, 75.0),
 (10, 0.03, 75.0)]

species_to_get = ['H2O']
temperature = list(np.arange(400,1000,10)) #this will find the nearest point for you in the db
pressure = [1]*len(temperature) #in bars

#NOTE: temperature and pressure lists are defined as pairs. The list above is
#grabbing temperatures from 400-1000 all at 1 bar.

data  = opa.get_molecular(db_filename, species_to_get, temperature,pressure)

data.keys()

dict_keys(['H2O', 'wavenumber'])

#note these temperatures might be different from your input
#if there isn't an exact point matching your input
data['H2O'].keys()

dict_keys([400.0, 425.0, 450.0, 475.0, 500.0, 525.0, 550.0, 575.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0])

data['H2O'][400.0].keys() #in this case we had the exact point (1.0 bar)

dict_keys([1.0])

Spectral11  = jpi.pals.Spectral11
f = jpi.figure(y_axis_type='log',height=500,y_range=[1e-29,1e-17],
          y_axis_label='Cross Section (cm2/species)', x_axis_label='Micron')
for T , C in zip(data['H2O'].keys(), Spectral11):
    x,y = jdi.mean_regrid(data['wavenumber'],data['H2O'][T][1.0], R=200)
    f.line(1e4/x,y, color=C,legend_label=str(T))
jpi.show(f)

Get Continuum Opacity

Very similar to molecular opacity, except not pressure dependent

molecules, temperatures = opa.continuum_avail(db_filename)

print(molecules)

['H-bf', 'H-ff', 'H2-', 'H2H2', 'H2He']

Note: The continuum is perhaps the most “hard-coded” aspect of this code. If you want to add new continuum molecules you must also update this part of picaso:

https://github.com/natashabatalha/picaso/blob/8ae68dfb0bd7c876f5f4fe7ee1d2ac00366d15ef/picaso/atmsetup.py#L313

data  = opa.get_continuum(db_filename, ['H2He'], list(temperature))

data.keys()

dict_keys(['H2He', 'wavenumber'])

data['H2He'].keys()

dict_keys([400.0, 425.0, 450.0, 475.0, 500.0, 525.0, 550.0, 575.0, 600.0, 625.0, 650.0, 675.0, 700.0, 725.0, 750.0, 775.0, 800.0, 825.0, 850.0, 875.0, 900.0, 925.0, 950.0, 975.0, 1000.0])

f = jpi.figure(y_axis_type='log',height=300,y_range=[1e-15,1e-5],
          y_axis_label='Cross Section (cm-1 amagat-2)', x_axis_label='Micron')
for T , C in zip(list(data['H2He'].keys())[::2], Spectral11):
    x,y = jdi.mean_regrid(data['wavenumber'],data['H2He'][T], R=200)
    f.line(1e4/x,y, color=C,legend_label=str(T))

jpi.show(f)