Opacity Factory¶
PICASO currently comes with single opacity database that includes both the continuum and molecular opacity. The opacity file is already included in the reference/ directory. You may have to install Github Large File System if you notice that opacity.db is downloading as an empty file.
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
[1]:
import warnings
warnings.filterwarnings('ignore')
import io
import numpy as np
import os
import picaso.opacity_factory as opa
#plotting
from bokeh.io import output_notebook
output_notebook()
from bokeh.plotting import show,figure
General Sqlite3 and how our database is structured¶
3 Tables:
header,continuum, andmolecularheader: 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¶
[2]:
#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','opacity.db')
What Species, Pressures and Temperatures are Available¶
[3]:
molecules, pt_pairs = opa.molecular_avail(db_filename)
[4]:
print(molecules)
['C2' 'C2H2' 'CH4' 'CO' 'CO2' 'CaH' 'FeH' 'H2O' 'H2S' 'HCN' 'K' 'MgH' 'N2'
'N2O' 'NH3' 'NO' 'Na' 'O2' 'O3' 'OCS' 'OH' 'PH3' 'SH' 'SO2' 'TiO' 'VO']
[5]:
pt_pairs[0:10] #full list of all the pt pairs. the first value is just an index for bookkeeping
[5]:
[(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)]
[6]:
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)
[7]:
data.keys()
[7]:
dict_keys(['H2O', 'wavenumber'])
[8]:
#note these temperatures might be different from your input
#if there isnt an exact point matching your input
data['H2O'].keys()
[8]:
dict_keys([400.0, 425.0, 450.0, 475.0, 500.0, 575.0, 650.0, 725.0, 800.0, 900.0, 1000.0])
[9]:
data['H2O'][400.0].keys() #in this case we had the exact point (1.0 bar)
[9]:
dict_keys([1.0])
[10]:
from bokeh.palettes import Spectral11
f = figure(y_axis_type='log',height=300,y_range=[1e-29,1e-20],
y_axis_label='Cross Section (cm2/species)', x_axis_label='Micron')
for T , C in zip(data['H2O'].keys(), Spectral11):
f.line(1e4/data['wavenumber'],data['H2O'][T][1.0], color=C,legend_label=str(T))
show(f)
Get Continuum Opacity¶
Very similar to molecular opacity, except not pressure dependent
[11]:
molecules, temperatures = opa.continuum_avail(db_filename)
[12]:
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:
[13]:
data = opa.get_continuum(db_filename, ['H2H2'], list(temperature))
[14]:
data.keys()
[14]:
dict_keys(['H2H2', 'wavenumber'])
[15]:
data['H2H2'].keys()
[15]:
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])
[16]:
from bokeh.palettes import Spectral11
f = figure(y_axis_type='log',height=300,y_range=[1e-35,1e-5],
y_axis_label='Cross Section (cm-1 amagat-2)', x_axis_label='Micron')
for T , C in zip(list(data['H2H2'].keys())[::2], Spectral11):
f.line(1e4/data['wavenumber'],data['H2H2'][T], color=C,legend_label=str(T))
show(f)