Note

Download full notebook here

Running the Code

In this tutorial you will learn:

  1. How to run the code!

You should already be familiar with:

  1. How to pick condensates to run in a cloud model

  2. What is \(f_{sed}\) and what number is right for my model?

  3. What are the chemical limitations in the code

  4. How to compute initial Mie scattering grid

[1]:

from bokeh.io import output_notebook
from bokeh.plotting import show, figure
from bokeh.palettes import Colorblind
output_notebook()
import numpy as np
import pandas as pd
import astropy.units as u

#cloud code
import virga.justdoit as jdi
Loading BokehJS ...

Simple Isothermal/Constant \(K_{z}\) Example

Let’s take the isothermal example from the first tutorial to begin gaining intuition for using the code

[2]:

pressure = np.logspace(-5,3,30) #simple isotherml PT profile (kelvin)
temperature = np.zeros(30)+1600
kz = np.zeros(30)+ 1e10 #constant kz profile

metallicity = 1 #atmospheric metallicity relative to Solar
mean_molecular_weight = 2.2 # atmospheric mean molecular weight

#get pyeddy recommendation for which gases to run
recommended_gases = jdi.recommend_gas(pressure, temperature,
                                     metallicity,mean_molecular_weight)

print(recommended_gases)

['CaAl12O19', 'CaTiO3', 'Cr', 'Fe', 'Mg2SiO4', 'MgSiO3', 'SiO2', 'TiO2']

Define the Atmosphere Class

[3]:

#set cloud species to condense
#let's take the recommended gases at face value for now
sum_planet = jdi.Atmosphere(['MgSiO3'],fsed=0.1,mh=metallicity,
                 mmw = mean_molecular_weight)

Set gravity and the P/T/Kz profile

[4]:

#set the planet gravity
sum_planet.gravity(gravity=357.00, gravity_unit=u.Unit('cm/(s**2)'))

#PT
sum_planet.ptk(df = pd.DataFrame({'pressure':pressure, 'temperature':temperature,
                           'kz':kz}))

[5]:

#you can also read in a file using pandas.read_csv()
pd.DataFrame({'pressure':pressure, 'temperature':temperature,
                           'kz':kz}).to_csv('test.csv')
sum_planet.ptk(filename='test.csv',usecols = [1,2,3])

The minimum values for kz

There is a parameter called kz_min in the ptk function. The default value is kz_min=1e5. This number isn’t necessarily a hard and fast rule, but it does avoid numerical instabilities. Another thing to note is that this excludes negative kz values. The code will alert you if you are tripping the kz_min reset.

[6]:

#set the planet gravity
sum_planet.gravity(gravity=357.00, gravity_unit=u.Unit('cm/(s**2)'))

#This is what happens when you trip up the minimum KZ value
sum_planet.ptk(df = pd.DataFrame({'pressure':pressure, 'temperature':temperature,
                           'kz':kz*-1}))

Run the code

[7]:

#directory where mieff files are
mieff_directory = '/data/virga/'
sum_planet.ptk(filename='test.csv',usecols = [1,2,3])
#start with the simplest output
#get total opacity, single scattering, asymmetry, and individual optical depths
df_out = sum_planet.compute(directory=mieff_directory)

Not doing sublayer as cloud deck at the bottom of pressure grid

[8]:

#get full dictionary output
all_out = sum_planet.compute(directory=mieff_directory, as_dict=True)

Not doing sublayer as cloud deck at the bottom of pressure grid

Exploring dict Output

There are many things to explore in the dict output. In the next tutorial, we will reproduce some of the most common plots used in papers to analyze cloud runs.

[9]:

#see what outputs exist in the dictionary
all_out.keys()

[9]:

dict_keys(['pressure', 'pressure_unit', 'temperature', 'temperature_unit', 'wave', 'wave_unit', 'condensate_mmr', 'cond_plus_gas_mmr', 'mean_particle_r', 'droplet_eff_r', 'r_units', 'column_density', 'column_density_unit', 'opd_per_layer', 'single_scattering', 'asymmetry', 'opd_by_gas', 'condensibles', 'scalar_inputs', 'fsed', 'altitude', 'layer_thickness', 'z_unit', 'mixing_length', 'mixing_length_unit', 'kz', 'kz_unit', 'scale_height', 'cloud_deck'])

Hopefully the names of the dict elements are self explanatory. If not, here is a brief description of each. nl=number of layers, nw=number of wavelengths, ng=number of gases.

  • condensate_mmr(nl x ng): Mean mass mixing ratio (concentration) of the condensate. See qc in A&M 01 Eqn. 4 and Eqn. 8.

  • cond_plus_gas_mmr(nl x ng): Mean mass mixing ratio of the consensate plus the vapor. See qt in A&M 01 Eqn. 7.

  • mean_particle_r(nl x ng): Geometric mean particle radius. See r_g in A&M 01 Eqn. 13.

  • droplet_eff_r(nl x ng): Effetive (area-weighted) droplet radius. See r_eff in A&M 01 Eqn. 17

  • column_density(nl x ng): Total number concentration of particles PER LAYER. See N in A&M 01 Eqn. 14

  • opd_per_layer(nl x nw): Total extinction PER layer. This includes all gases.

  • single_scattering(nl x nw): Total single scattering albedo

  • asymmetry(nl x nw): Total asymmetry

  • opd_by_gas(nl x ng): Optical depth for conservative geometric scatteres separated by gas. E.g. Fig 7 in Morley+2012

Additionally, scalar input contains some of the original input

[10]:

all_out['scalar_inputs']

[10]:

{'mh': 1, 'mmw': 2.2, 'sig': 2.0, 'nrad': 60, 'rmin': 1e-08}

[ ]: