Note

Download full notebook here

Getting Started : Basic Inputs and Outputs

If you are here then you have already successfully

1. Installed the code

2. Downloaded the necessary reference data

3. Added export picaso_refdata="/path/to/picaso/reference" to ~/.bash_profile

If you have not done these things, please return to Installation Guilde

#picaso
from picaso import justdoit as jdi
from picaso import justplotit as jpi
import numpy as np

jpi.output_notebook()

Loading BokehJS ...

#double check that your reference file path has been set
import os
refdata = os.getenv("picaso_refdata")
print(refdata)
#if you are having trouble setting this you can do it right here in the command line
#os.environ["picaso_refdata"]= add your path here AND COPY AND PASTE ABOVE
#IT WILL NEED TO GO ABOVE YOUR PICASO IMPORT

/data/reference_data/picaso/reference

Connect to Opacity Database

There is a full notebook in the tutorials devoted to learning how to make opacity databases in our format. If you cloned from github, there should also be an opacity database there called opacity.db.

help(jdi.opannection)

Help on function opannection in module picaso.justdoit:

opannection(wave_range=None, filename_db=None, raman_db=None, resample=1, ck_db=None, deq=False, on_fly=False, gases_fly=None, ck=False, verbose=True)
    Sets up database connection to opacities.

    Parameters
    ----------
    wave_range : list of float
        Subset of wavelength range for which to run models for
        Default : None, which pulls entire grid
    filename_db : str
        Filename of opacity database to query from
        Default is none which pulls opacity file that comes with distribution
    raman_db : str
        Filename of raman opacity cross section
        Default is none which pulls opacity file that comes with distribution
    resample : int
        Default=1 (no resampling) PROCEED WITH CAUTION!!!!!This will resample your opacites.
        This effectively takes opacity[::BINS] depending on what the
        sampling requested is. Consult your local theorist before
        using this.
    ck_db : str
        ASCII filename of ck file
    verbose : bool
        Error message to warn users about resampling. Can be turned off by supplying
        verbose=False

opannection has a few default parameters. A few notes:

1. wave_range can be used to select a subset of the full opacity database you are using

2. filename_db is a sqlite database that should have been downloaded and stored with your reference data (see Installation Documentation)

3. raman_db is a small file that should already be in your reference folder from Github.

opacity = jdi.opannection(wave_range=[0.3,1]) #lets just use all defaults

Load blank slate

start_case = jdi.inputs()

In order to run the code we need (at the minimum) specific info about the:

• phase angle

• planet : gravity

• star : temperature, metallicity, gravity

• atmosphere : P-T profile, chemical composition, cloud properties (discussed later)

Additionally, we have some optional parameters that will be discussed in later notebooks:

• approx : approximations to rapidly compute reflected light

• disco : number of individual facets you want to compute before integrating to 1D spectrum (think disco ball)

• opacities : keeping track of opacity files used

• test_mode : used to do benchmark testing

Set Planet & Star Properties

#phase angle
start_case.phase_angle(0) #radians

#define gravity
start_case.gravity(gravity=25, gravity_unit=jdi.u.Unit('m/(s**2)')) #any astropy units available

#define star
start_case.star(opacity, 5000,0,4.0) #opacity db, pysynphot database, temp, metallicity, logg

Set Atmospheric Composition

There are different options for setting atmospheric composition.

1) Specifying a file path to model run
2) Give arbitrary pressure, temperature and composition directly as a dictionary input

Option 1) Specify file path

Below, I am loading in a profile path for Jupiter that should be included in your reference data

print(jdi.jupiter_pt()) #should return the path to your reference data

/data/reference_data/picaso/reference/base_cases/jupiter.pt

start_case.atmosphere(filename=jdi.jupiter_pt(), delim_whitespace=True)

/home/nbatalh1/codes/picaso/picaso/justdoit.py:1742: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
  df = pd.read_csv(filename, **pd_kwargs)

File format

1. Must have pressure(bars), temperature(Kelvin), and case sensitive molecule names (e.g. TiO, Na, H2O, etc) for mixing ratios (in no particular order)

2. Can specify any necessary key word arguments for pd.read_csv at the end

PICASO will auto-compute mixing ratios, determine what CIA is neceesary and compute mean molecular weight based on these headers. Take at the preloaded example below

#to give you an idea
comp_file = jdi.pd.read_csv(jdi.jupiter_pt(), delim_whitespace=True)
#see example below
comp_file.head()

/tmp/ipykernel_107841/325953466.py:2: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
  comp_file = jdi.pd.read_csv(jdi.jupiter_pt(), delim_whitespace=True)

	pressure	temperature	e-	H2	H	H+	H-	VO	TiO	CO2	He	H2O	CH4	CO	NH3	N2	PH3
0	0.000001	150.87	4.500000e-38	0.837	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	0.163	0.000069	0.000466	4.500000e-38	0.000137	5.420000e-17	4.500000e-38
1	0.000001	149.68	4.500000e-38	0.837	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	0.163	0.000035	0.000466	4.500000e-38	0.000137	1.580000e-17	4.500000e-38
2	0.000002	148.40	4.500000e-38	0.837	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	0.163	0.000017	0.000466	4.500000e-38	0.000137	4.370000e-18	4.500000e-38
3	0.000003	147.00	4.500000e-38	0.837	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	0.163	0.000008	0.000466	4.500000e-38	0.000137	1.140000e-18	4.500000e-38
4	0.000004	145.46	4.500000e-38	0.837	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	4.500000e-38	0.163	0.000004	0.000466	4.500000e-38	0.000137	2.740000e-19	4.500000e-38

Create 1D Reflected Light Spectrum

Let’s create our first spectrum of Jupiter’s reflected light at full phase

df = start_case.spectrum(opacity)

Checkout out what was returned (Note this is a change in v1.0)

df.keys()

dict_keys(['wavenumber', 'albedo', 'bond_albedo', 'fpfs_reflected'])

Regrid Opacities to Constant Resolution

wno, alb, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected']
wno, alb = jdi.mean_regrid(wno, alb , R=150)

jpi.show(jpi.spectrum(wno, alb, plot_width=500,x_range=[0.3,1]))

FpFs is the relative flux of the planet and star or :

\(\frac{f_p}{f_s} = a_\lambda \left( \frac{ r_p}{a} \right) ^2\)

where \(a\) is the semi-major axis. You may have noticed that we did not supply a radius or semi-major axis in the above code. Therefore, if you print out \(\frac{f_p}{f_s}\) you will see this:

fpfs

['Semi-major axis not supplied. If you want fpfs, add it to `star` function. ',
 'Planet Radius not supplied. If you want fpfs, add it to `gravity` function with a mass.']

Let’s add this to the star function so we can get the relative flux as well..

start_case.star(opacity, 5000,0,4.0,semi_major=1, semi_major_unit=jdi.u.Unit('au'))
start_case.gravity(radius=1, radius_unit=jdi.u.Unit('R_jup'),
                   mass = 1, mass_unit=jdi.u.Unit('M_jup'))
df = start_case.spectrum(opacity)
wno, alb, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected']

fpfs

array([9.56281071e-10, 7.27078706e-10, 7.39534725e-10, ...,
       1.44657210e-07, 1.44654860e-07, 1.44651937e-07])

Option 2) Arbitrary PT and Chemistry

Sometimes for testing (or for atmospheres that we don’t fully understand) an isothermal, well-mixed profile is sufficient. If we don’t want to load in a full profile, we can give it a simple DataFrame with the info we need.

start_case.atmosphere( df = jdi.pd.DataFrame({'pressure':np.logspace(-6,2,60),
                                                 'temperature':np.logspace(-6,2,60)*0+200,
                                                 "H2":np.logspace(-6,2,60)*0+0.837,
                                                 "He":np.logspace(-6,2,60)*0+0.163,
                                                 "CH4":np.logspace(-6,2,60)*0+0.000466})
                     )

df = start_case.spectrum(opacity)
wno_ch4, alb_ch4, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected']
wno_ch4, alb_ch4 = jdi.mean_regrid(wno_ch4, alb_ch4 , R=150)

jpi.show(jpi.spectrum(wno_ch4, alb_ch4, plot_width=500,x_range=[0.3,1]))

See how the plot above is much easier to interpret than the one with the full set of molecular input. Here we can clearly see the effects of methane opacity, raman and rayleigh scattering (the next notebook will include a tutorial for more diagnostic plotting)

Diagnostic help: Sometimes it helps to exclude molecule to see how it is influencing the spectrum

Take a look below

start_case.atmosphere(filename=jdi.jupiter_pt(), exclude_mol='H2O', delim_whitespace=True)

df = start_case.spectrum(opacity)
wno_nowater, alb_nowater, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected']
wno_nowater, alb_nowater= jdi.mean_regrid(wno_ch4, alb_ch4 , R=150)
fig = jpi.spectrum(wno, alb, plot_width=500)
fig.line(1e4/wno_nowater, alb_nowater, line_width=2, color='red')
jpi.show(fig)

/home/nbatalh1/codes/picaso/picaso/justdoit.py:1742: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
  df = pd.read_csv(filename, **pd_kwargs)