Note

Download full notebook here

Querying Opacities

PICASO currently comes with single opacity database that includes both the continuum and molecular opacity. The opacity file that comes included with has a 0.3-14 \(\mu\)m opacity grid sampled at R=15000. This is why we consistently bin our spectra to R~100. Contact us with opacity needs.

We chose those use sqlite3 for our database because of it’s 1) user-friendliness, 2) speed, 3) scalability, 4) compatibility with parallel programing. We tried out various other methods as well– json, hdf5, ascii, sqlalchemy– but sqlite3 was truly better for this specific problem.

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

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

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

#this is where your opacity file should be located if you've set your environments correctly
db_filename = os.path.join(os.getenv('picaso_refdata'), 'opacities','opacities.db')

What Species, Pressures and Temperatures are Available

molecules, pt_pairs = opa.molecular_avail(db_filename)

print(molecules)

['AlH', 'C2H2', 'C2H4', 'C2H6', 'CH4', 'CO', 'CO2', 'CaH', 'CrH', 'Cs', 'Fe', 'FeH', 'H2', 'H2O', 'H2S', 'H3+', 'HCN', 'K', 'Li', 'LiCl', 'LiF', 'LiH', 'MgH', 'N2', 'N2O', 'NH3', 'Na', 'O2', 'O3', 'OCS', 'PH3', 'Rb', 'SiO', 'TiH', 'TiO', 'VO']

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 isnt 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-', 'H2CH4', 'H2H', 'H2H2', 'H2He', 'H2N2']

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)