The Code

Submodules

picaso.atmsetup module

class picaso.atmsetup.ATMSETUP(config)

Bases: object

Reads in default source configuration from JSON and creates a full atmosphere class.

• Gets PT profile

• Computes mean molecular weight, density, column density, mixing ratios

• Gets cloud profile

• Gets stellar profile

add_warnings(warn)

Accumulate warnings from ATMSETUP

as_dict()

Get output into picklable dict format

calc_PT(logPc, T, logKir, logg1)

Calculates parameterized PT profile from Guillot. This isntance is here primary for the retrieval scheme, so this can be updated

Parameters:

• T (float) –

• logKir (float) –

• logg1 (float) –

• logPc (float) –

disect(g, t)

This disects the 3d input to a 1d input which is a function of a single gangle and tangle. This makes it possible for us to get the opacities at each facet before we go into the flux calculation

Parameters:

• g (int) – Gauss angle index

• t (int) – Tchebyshev angle index

get_altitude(p_reference=1, constant_gravity=False)

Calculates z and gravity

Parameters:

• p_reference (float) – Reference pressure for radius. (bars)

• constant_gravity (bool) – creates zero altitude dependence in the height

get_clouds(wno)

Get cloud properties from .cld input returned from eddysed. The eddysed file should have the following specifications

1. Have the following column names (any order) opd g0 w0

2. Be white space delimeted

3. Has to have values for pressure levels (N) and wwavenumbers (M). The row order should go:

level wavenumber opd w0 g0 1. 1. … . . 1. 2. … . . 1. 3. … . . . . … . . . . … . . 1. M. … . . 2. 1. … . . . . … . . N. . … . .

Warning

The order of the rows is very important because each column will be transformed into a matrix that has the size: [nlayer,nwave].

Parameters:

wno (array) – Array in ascending order of wavenumbers. This is used to regrid the cloud output

Todo

• Take regridding out of here and add it to justdoit

get_column_density()

Calculates the column desntiy based on TP profile: LAYER unit = g/cm2 = pressure/gravity =dyne/cm2/gravity = (g*cm/s2)/cm2/(cm/s2)=g/cm2

get_constants()

This function gets all conversion factors based on input units and all constants. Goal here is to get everything to cgs and make sure we are only call lengths of things once (e.g. number of wavelength points, number of layers etc)

Some stuff will be added to this in other routines (e.g. number of layers)

get_density()

Calculates density of atmospheres used on TP profile: LEVEL units of cm-3

get_mmw()

Returns the mean molecular weight of the atmosphere

get_needed_continuum(available_ray_mol, available_continuum)

This will define which molecules are needed for the continuum opacities. THis is based on temperature and molecules. This is terrible code I

Parameters:

• available_ray_mol (list of str) – list of available rayleigh molecules

• available_continuum (list of str) – list of available continuum molecules

get_profile()

Get profile from file or from user input file or direct pandas dataframe. If PT profile is not given, use parameterization

Currently only needs inputs from config file

Todo

• Add regridding to this by having users be able to set a different nlevel than the input cloud code is

get_profile_3d()

Get profile from file or from user input file or direct pandas dataframe. If PT profile is not given, use parameterization

Currently only needs inputs from config file

Parameters:

• lat (array) – Array of latitudes (radians)

• lon (array) – Array of longitudes (radians)

Todo

• Add regridding to this by having users be able to set a different nlevel than the input cloud code is

get_surf_reflect(nwno)

Gets the surface reflectivity from input

get_weights(molecule)

Automatically gets mean molecular weights of any molecule. Requires that user inputs case sensitive molecules i.e. TiO instead of TIO.

Parameters:

molecule (str or list) – any molecule string e.g. “H2O”, “CH4” etc “TiO” or [“H2O”, ‘CH4’] Isotopologues should be specified with a dash separating the elemtents. E.g. 12C-16O2 or 13C-16O2. If you only want to specify one then you would have 12C-O2.

Returns:

molecule as keys with molecular weights as value

Return type:

dict

picaso.atmsetup.convert_to_simple(iso_name)

Converts iso name (e.g. 13C-16O2 to CO2) Returns same name if not (e.g. CO2 gives CO2)

picaso.atmsetup.separate_molecule_name(molecule_name)

used for separating molecules For example, CO2 becomes “C” “O2”

picaso.atmsetup.separate_string_number(string)

used for separating numbers from molecules For example, CO2 becomes “C” “O2” in separate_molecule_name then this function turns it into [[‘C’],[‘O’,’2’]]

picaso.build_3d_input module

picaso.build_3d_input.lon_lat_to_cartesian(lon, lat, R=1)

Calculates lon, lat coordinates of a point on a sphere with radius R

Parameters:

• lon (array) – Longitude

• lat (array) – Latitude

• R (float) – radius of sphere

Return type:

x, y, z

picaso.build_3d_input.make_3d_cld_input(ng, nt, phase_angle, input_file, output_file, lat_range=None, lon_range=None, rand_coverage=1)

Discontinued Program to create 3d CLOUD input. Used to feed GCM input into disco ball.

Parameters:

• ng (int) – Number of gauss angles

• nt (int) – Number of Tchebysehv angles

• phase_angle (float) – Geometry of phase angle

• input_file (str) – CLD input file you want to create the 3d grid from. Currently takes in a single 1d CLD profile. Input file must be structured with the following columns in the correct order: [‘lvl’, ‘wv’,’opd’,’g0’,’w0’,’sigma’] input file

• output_file (str) – Output file location

• lat_range (list) – (Optional)Range of latitudes to exclude in the CLD file

• lon_range (list) – (Optional)Range of longitudes to exclud in the CLD file

• rand_coverage (float) – (Optional)Fractional cloud coverage. rand_cov=1 introduces full cloud coverage. rand_cov=0 has no coverage.

• **kwargs (dict) – Key word arguments used for panads.read_csv

Return type:

Creates output file. No other returns.

Examples

Basic use of creating a 3D file:

>>> phase_angle = 0
>>> number_gauss_angles = 2
>>> number_tcheby_angles = 3
>>> phase_angle = 0
>>> input_file=jdi.jupiter_cld()
>>> output_file='jupiter3D.hdf5'
>>>make_3d_pt_input(number_gauss_angles, number_tcheby_angles,
>>>         phase_angle, input_file,output_file)

picaso.build_3d_input.make_3d_pt_input(ng, nt, phase_angle, input_file, output_file, **kwargs)

Discontinued Program to create 3d PT input. Used to feed GCM input into disco ball.

Parameters:

• ng (int) – Number of gauss angles

• nt (int) – Number of Tchebysehv angles

• phase_angle (float) – Geometry of phase angle

• input_file (str) – PT input file you want to create the 3d grid from. Currently takes in a single 1d PT profile. Must have labeled, whitespace delimeted columns. In the columns, must have temperature, pressure, and at least some molecular species (which are ALL case-sensitive)

• output_file (str) – Output file location

Return type:

Creates output file. No other returns.

Examples

Basic use of creating a 3D file:

>>> phase_angle = 0
>>> number_gauss_angles = 2
>>> number_tcheby_angles = 3
>>> phase_angle = 0
>>> input_file=jdi.jupiter_pt()
>>> output_file='jupiter3D.hdf5'
>>>make_3d_pt_input(number_gauss_angles, number_tcheby_angles,
>>>         phase_angle, input_file,output_file)

Can also use **kwargs to control how pandas reads in input file

>>>make_3d_pt_input(number_gauss_angles, number_tcheby_angles, >>> phase_angle, input_file,output_file,usecols=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,20,21])

picaso.build_3d_input.rebin_mitgcm_cld(ng, nt, phase_angle, input_file, output_file, names=['i', 'j', 'opd', 'g0', 'w0'])

Rebin post processed GCM cloud grid to a smaller grid. This function creates a pickle file that can be directly input into clouds_3d.

Parameters:

• ng (int) – Number of gauss angles

• nt (int) – Number of Tchebysehv angles

• phase_angle (float) – Geometry of phase angle (radians)

• input_file (str) – PT input file you want to create the 3d grid from. Currently takes in a single 1d PT profile. Must have labeled, whitespace delimeted columns. In the columns, must have temperature, pressure, and at least some molecular species (which are ALL case-sensitive)

• output_file (str) – output file location

• names (list,str,optional) – List of names for the post processed file. The Default value is what generallly comes out of A&M code.. But double check before running!!

picaso.build_3d_input.rebin_mitgcm_pt(ng, nt, phase_angle, input_file, output_file, p_unit='Pa', kzz_unit='m*m/s', run_chem=False, MH=None, CtoO=None)

Rebin GCM grid to a smaller grid. This function creates a pickle file that can be directly input into atmosphere_3d.

Parameters:

• ng (int) – Number of gauss angles

• nt (int) – Number of Tchebysehv angles

• phase_angle (float) – Geometry of phase angle (radians)

• input_file (str) – PT input file you want to create the 3d grid from. Currently takes in a single 1d PT profile. Must have labeled, whitespace delimeted columns. In the columns, must have temperature, pressure, and at least some molecular species (which are ALL case-sensitive)

• output_file (str) – output file location

• p_unit (str,optional) – Pressure Unit (default Pascal)

• kzz_unit (str,optional) – Kzz Unit (default is m*m/s)

• run_chem (bool , optional) – (Optional) runs chemistry and adds to json file

• run_chem – (Optional) MIT gcm output doesn’t usually come with chemistry in the files If this is the case then you can post process chemistry by setting this to True.

• MH (float , optional) – (Optional) This is only used if run_chem is set to True. It is the Metallicity of the planet in NON log units. MH = 1 is solar.

• CtoO (float , optional) – (Optional) This is the C/O ratio of the planet also in NON log units. CtoO = 0.55 is solar.

picaso.build_3d_input.regrid_xarray(xarray_original, num_gangle=None, num_tangle=None, phase_angle=None, latitude=None, longitude=None)

This regrid functionality uses the xarray and xesmf method. It requires to additional packages that are not required by the rest of PICASO (xesmf). See tutorial for additional tips on how to structure your xarray input.

Parameters:

• xarray_original (xarray) – xarray of all 3d input you would like regridded

• num_gangle (int) – Number of gauss angles (not needed if directly supply lat and lons)

• num_tangle (int) – Number of cheychev angles (not needed if directly supply lat and lons)

• phase_angle (int) – Radians phase angle (not needed if directly supply lat and lons)

• latitude (array) – array of latitudes in degrees (not needed if supplied num_gangle, num_tangle, phase_angle)

• longitude (array) – array of longitudes in degrees (not needed if supplied num_gangle, num_tangle, phase_angle)

• Output –

• ------ – Dictionary of regridded data and new latitude longitude grid

picaso.build_3d_input.variable_pressure_example()

picaso.disco module

picaso.disco.compress_disco(nwno, cos_theta, xint_at_top, gweight, tweight, F0PI)

Last step in albedo code. Integrates over phase angle based on the Gaussian-Chebychev weights in geometry.json

Parameters:

• nwno (int) – Number of wavenumbers

• cos_theta (float) – Cosine of phase angle

• xint_at_top (ndarray of floats) – Planetary intensity at the top of the atmosphere with dimensions (ng, nt, nwno)

• gweight (ndarray of floats) – Gaussian weights for integration

• tweight (ndarray of floats) – Chebychev weights for integration

• F0PI (ndarray of floats) – Stellar flux

picaso.disco.compress_thermal(nwno, flux_at_top, gweight, tweight)

Last step in albedo code. Integrates over phase angle based on the Gaussian-Chebychev weights in disco.get_angles_1d or 3d

Parameters:

• nwno (int) – Number of wavenumbers

• flux_at_top (ndarray of floats) – Thermal Flux at the top of the atmosphere with dimensions (ng, nt, nwno) or could also be (ng,nt,nlayer,nwno)

• gweight (ndarray of floats) – Gaussian weights for integration

• tweight (ndarray of floats) – Chebychev weights for integration

picaso.disco.compute_disco(ng, nt, gangle, tangle, phase_angle)

Computes ubar0, the incident angle, and ubar1, the outgoing angle from the chebyshev angles in geometry.json

Parameters:

• ng (int) – Number of gauss angles

• nt (int) – Number of tchebyshev angles

• gangle (float) – Gaussian Angles

• tangle (float) – Chevychev Angles

• phase_angle (float) – Planetary phase angle

Returns:

• ubar0 – the incident angles

• ubar1 – the outgoing angles

• cos_theta – Cosine of the phase angle

picaso.disco.get_angles_1d(ngauss)

Computes integration over half sphere. ngauss=5 is 5 angles over half a sphere.

These are taken from : Abramowitz and Stegun Table 25.8 https://books.google.com/books?id=32HJ1orVXDEC&pg=PA876&lpg=PA876&dq=Abramowitz+and+Stegun+Table+25.8&source=bl&ots=_eKu6k5SG8&sig=ACfU3U0_piSAi8aDHOVd0itf_QSAFxcFLQ&hl=en&sa=X&ved=2ahUKEwj5neLR2fDoAhXemHIEHUSVCjEQ6AEwD3oECAwQKQ#v=onepage&q=Abramowitz%20and%20Stegun%20Table%2025.8&f=false

Parameters:

ngauss (int) – Can only be 5, 6, 7, 8

picaso.disco.get_angles_3d(num_gangle, num_tangle)

Computes angles for full disco ball

Parameters:

• num_gangles (int) – Number of Gauss angles

• num_tangles (int) – Number of Tchebychev angles desired

Returns:

Gauss Angles,Gauss Weights,Tchebyshev Angles,Tchebyshev weights

Return type:

np.ndarray, np.ndarray, np.ndarray, np.ndarray

picaso.elements module

Properties of the chemical elements.

Each chemical element is represented as an object instance. Physicochemical and descriptive properties of the elements are stored as instance attributes.

Author:

Christoph Gohlke

Version:

2015.01.29

Requirements

• CPython 2.7 or 3.4

References

1. http://physics.nist.gov/PhysRefData/Compositions/

2. http://physics.nist.gov/PhysRefData/IonEnergy/tblNew.html

3. http://en.wikipedia.org/wiki/%(element.name)s

4. http://www.miranda.org/~jkominek/elements/elements.db

Examples

>>> from elements import ELEMENTS
>>> len(ELEMENTS)
109
>>> str(ELEMENTS[109])
'Meitnerium'
>>> ele = ELEMENTS['C']
>>> ele.number, ele.symbol, ele.name, ele.eleconfig
(6, 'C', 'Carbon', '[He] 2s2 2p2')
>>> ele.eleconfig_dict
{(1, 's'): 2, (2, 'p'): 2, (2, 's'): 2}
>>> sum(ele.mass for ele in ELEMENTS)
14659.1115599
>>> for ele in ELEMENTS:
...     ele.validate()
...     ele = eval(repr(ele))

picaso.fluxes module

picaso.fluxes.blackbody(t, w)

Blackbody flux in cgs units in per unit wavelength (cm)

Parameters:

• t (array,float) – Temperature (K)

• w (array, float) – Wavelength (cm)

Return type:

ndarray with shape ntemp x numwave in units of erg/cm/s2/cm

picaso.fluxes.blackbody_climate_deprecate(wave, temp, bb, y2, tp, tmin, tmax)

picaso.fluxes.blackbody_integrated(T, wave, dwave)

This computes the total energey per wavenumber bin needed for the climate calculation Note that this is different than the raw flux at an isolated wavenumber. Therefore this function is different than the blackbody function in picaso.fluxes which computes blackbody in raw cgs units.

Parameters:

• T (float, array) – temperature in Kelvin

• wave (float, array) – wavenumber in cm-1

• dwave (float, array) – Wavenumber bins in cm-1

Returns:

num temperatures by num wavenumbers units of ergs/cm*2/s/cm-1 for integrated bins ()

Return type:

array

picaso.fluxes.calculate_flux(F, G, X)

Calculate fluxes

picaso.fluxes.chapman(pressure, pm, hratio)

Computes Chapman function for use in tidal routine

Parameters

pressurefloat

pressure

pmfloat

Some pressure (?)

hratiofloat

Ratio of Scale Height over Chapman Scale Height

Returns

Chapman function

picaso.fluxes.did_grad_cp(t, p, t_table, p_table, grad, cp, calc_type)

Parameters:

• t (float) – Temperature value

• p (float) – Pressure value

• t_table (array) – array of Temperature values with 53 entries

• p_table (array) – array of Pressure value with 26 entries

• grad (array) – array of gradients of dimension 53*26

• cp (array) – array of cp of dimension 53*26

• calc_type (int) – not used to make compatible with nopython.

Returns:

grad_x,cp_x

Return type:

float

picaso.fluxes.get_kzz(pressure, temp, grav, mmw, tidal, flux_net_ir_layer, flux_plus_ir_attop, t_table, p_table, grad, cp, calc_type, nstr)

picaso.fluxes.get_reflected_1d(nlevel, wno, nwno, numg, numt, dtau, tau, w0, cosb, gcos2, ftau_cld, ftau_ray, dtau_og, tau_og, w0_og, cosb_og, surf_reflect, ubar0, ubar1, cos_theta, F0PI, single_phase, multi_phase, frac_a, frac_b, frac_c, constant_back, constant_forward, get_toa_intensity=1, get_lvl_flux=0, toon_coefficients=0, b_top=0)

Computes toon fluxes given tau and everything is 1 dimensional. This is the exact same function as get_flux_geom_3d but is kept separately so we don’t have to do unecessary indexing for fast retrievals. :param nlevel: Number of levels in the model :type nlevel: int :param wno: Wave number grid in cm -1 :type wno: array of float :param nwno: Number of wave points :type nwno: int :param numg: Number of Gauss angles :type numg: int :param numt: Number of Chebyshev angles :type numt: int :param DTAU: This is the opacity contained within each individual layer (defined at midpoints of “levels”)

WITHOUT D-Eddington Correction Dimensions=# layer by # wave

Parameters:

• TAU (ndarray of float) – This is the cumulative summed opacity WITHOUT D-Eddington Correction Dimensions=# level by # wave

• W0 (ndarray of float) – This is the single scattering albedo, from scattering, clouds, raman, etc WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• COSB (ndarray of float) – This is the asymmetry factor WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• GCOS2 (ndarray of float) – Parameter that allows us to directly include Rayleigh scattering = 0.5*tau_rayleigh/(tau_rayleigh + tau_cloud)

• ftau_cld (ndarray of float) – Fraction of cloud extinction to total = tau_cloud/(tau_rayleigh + tau_cloud)

• ftau_ray (ndarray of float) – Fraction of rayleigh extinction to total = tau_rayleigh/(tau_rayleigh + tau_cloud)

• dtau_og (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT the delta eddington correction, if it was specified by user Dimensions=# layer by # wave

• tau_og (ndarray of float) – This is the cumulative summed opacity WITHOUT the delta eddington correction, if it was specified by user Dimensions=# level by # wave

• w0_og (ndarray of float) – Same as w0 but WITHOUT the delta eddington correction, if it was specified by user

• cosb_og (ndarray of float) – Same as cosbar buth WITHOUT the delta eddington correction, if it was specified by user

• surf_reflect (float) – Surface reflectivity

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• ubar1 (ndarray of float) – matrix of cosine of the observer angles

• cos_theta (float) – Cosine of the phase angle of the planet

• F0PI (array) – Downward incident solar radiation

• single_phase (str) – Single scattering phase function, default is the two-term henyey-greenstein phase function

• multi_phase (str) – Multiple scattering phase function, defulat is N=2 Legendre polynomial

• frac_a (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_b (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_c (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C), Default is : 1 - gcosb^2

• constant_back (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of back scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• constant_forward (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of forward scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• get_toa_intensity (int) – (Optional) Default=1 is to only return the TOA intensity you would need for a 1D spectrum (1) otherwise it will return zeros for TOA intensity

• get_lvl_flux (int) – (Optional) Default=0 is to only compute TOA intensity and NOT return the lvl fluxes so this needs to be flipped on for the climate calculations

• toon_coefficients (int) – (Optional) 0 for quadrature (default) 1 for eddington

Return type:

intensity at the top of the atmosphere for all the different ubar1 and ubar2

picaso.fluxes.get_reflected_1d_deprecate(nlevel, wno, nwno, numg, numt, dtau, tau, w0, cosb, gcos2, ftau_cld, ftau_ray, dtau_og, tau_og, w0_og, cosb_og, surf_reflect, ubar0, ubar1, cos_theta, F0PI, single_phase, multi_phase, frac_a, frac_b, frac_c, constant_back, constant_forward, toon_coefficients=0, b_top=0)

Computes toon fluxes given tau and everything is 1 dimensional. This is the exact same function as get_flux_geom_3d but is kept separately so we don’t have to do unecessary indexing for fast retrievals. :param nlevel: Number of levels in the model :type nlevel: int :param wno: Wave number grid in cm -1 :type wno: array of float :param nwno: Number of wave points :type nwno: int :param numg: Number of Gauss angles :type numg: int :param numt: Number of Chebyshev angles :type numt: int :param DTAU: This is the opacity contained within each individual layer (defined at midpoints of “levels”)

WITHOUT D-Eddington Correction Dimensions=# layer by # wave

Parameters:

• TAU (ndarray of float) – This is the cumulative summed opacity WITHOUT D-Eddington Correction Dimensions=# level by # wave

• W0 (ndarray of float) – This is the single scattering albedo, from scattering, clouds, raman, etc WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• COSB (ndarray of float) – This is the asymmetry factor WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• GCOS2 (ndarray of float) – Parameter that allows us to directly include Rayleigh scattering = 0.5*tau_rayleigh/(tau_rayleigh + tau_cloud)

• ftau_cld (ndarray of float) – Fraction of cloud extinction to total = tau_cloud/(tau_rayleigh + tau_cloud)

• ftau_ray (ndarray of float) – Fraction of rayleigh extinction to total = tau_rayleigh/(tau_rayleigh + tau_cloud)

• dtau_og (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT the delta eddington correction, if it was specified by user Dimensions=# layer by # wave

• tau_og (ndarray of float) – This is the cumulative summed opacity WITHOUT the delta eddington correction, if it was specified by user Dimensions=# level by # wave

• w0_og (ndarray of float) – Same as w0 but WITHOUT the delta eddington correction, if it was specified by user

• cosb_og (ndarray of float) – Same as cosbar buth WITHOUT the delta eddington correction, if it was specified by user

• surf_reflect (float) – Surface reflectivity

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• ubar1 (ndarray of float) – matrix of cosine of the observer angles

• cos_theta (float) – Cosine of the phase angle of the planet

• F0PI (array) – Downward incident solar radiation

• single_phase (str) – Single scattering phase function, default is the two-term henyey-greenstein phase function

• multi_phase (str) – Multiple scattering phase function, defulat is N=2 Legendre polynomial

• frac_a (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_b (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_c (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C), Default is : 1 - gcosb^2

• constant_back (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of back scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• constant_forward (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of forward scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• toon_coefficients (int) – 0 for quadrature (default) 1 for eddington

Returns:

• intensity at the top of the atmosphere for all the different ubar1 and ubar2

• To Do

• —–

• - F0PI Solar flux shouldn’t always be 1.. Follow up to make sure that this isn’t a bad – hardwiring to solar, despite “relative albedo”

picaso.fluxes.get_reflected_1d_gfluxv_deprecate(nlevel, wno, nwno, numg, numt, dtau, tau, w0, cosb, surf_reflect, b_top, b_surface, ubar0, F0PI, tridiagonal, delta_approx)

Computes upwelling and downwelling layer and level toon fluxes given tau and everything is 1 dimensional. This is the exact same function as `GFLUXV.f’. retrievals. :param nlevel: Number of levels in the model :type nlevel: int :param wno: Wave number grid in cm -1 :type wno: array of float :param nwno: Number of wave points :type nwno: int :param numg: Number of Gauss angles :type numg: int :param numt: Number of Chebyshev angles :type numt: int :param DTAU: This is the opacity contained within each individual layer (defined at midpoints of “levels”)

WITHOUT D-Eddington Correction Dimensions=# layer by # wave

Parameters:

• TAU (ndarray of float) – This is the cumulative summed opacity WITHOUT D-Eddington Correction Dimensions=# level by # wave

• W0 (ndarray of float) – This is the single scattering albedo, from scattering, clouds, raman, etc WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• COSB (ndarray of float) – This is the asymmetry factor WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• surf_reflect (float) – Surface reflectivity

• b_top (float) – Top Boundary Conditions

• b_surface (float) – Surface Boundary Conditions

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• F0PI (array) – Downward incident solar radiation

• delta_approx (int) – 0 for no Delta Approx, 1 for Delta Approx

Returns:

• intensity at the top of the atmosphere for all the different ubar1 and ubar2

• To Do

• —–

• - F0PI Solar flux shouldn’t always be 1.. Follow up to make sure that this isn’t a bad – hardwiring to solar, despite “relative albedo”

picaso.fluxes.get_reflected_3d(nlevel, wno, nwno, numg, numt, dtau_3d, tau_3d, w0_3d, cosb_3d, gcos2_3d, ftau_cld_3d, ftau_ray_3d, dtau_og_3d, tau_og_3d, w0_og_3d, cosb_og_3d, surf_reflect, ubar0, ubar1, cos_theta, F0PI, single_phase, multi_phase, frac_a, frac_b, frac_c, constant_back, constant_forward)

Computes toon fluxes given tau and everything is 3 dimensional. This is the exact same function as get_flux_geom_1d but is kept separately so we don’t have to do unecessary indexing for retrievals. :param nlevel: Number of levels in the model :type nlevel: int :param wno: Wave number grid in cm -1 :type wno: array of float :param nwno: Number of wave points :type nwno: int :param numg: Number of Gauss angles :type numg: int :param numt: Number of Chebyshev angles :type numt: int :param dtau_3d: This is the opacity contained within each individual layer (defined at midpoints of “levels”)

WITHOUT D-Eddington Correction Dimensions=# layer by # wave

Parameters:

• tau_3d (ndarray of float) – This is the cumulative summed opacity WITHOUT D-Eddington Correction Dimensions=# level by # wave

• w0_3d (ndarray of float) – This is the single scattering albedo, from scattering, clouds, raman, etc WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• cosb_3d (ndarray of float) – This is the asymmetry factor WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• gcos2_3d (ndarray of float) – Parameter that allows us to directly include Rayleigh scattering = 0.5*tau_rayleigh/(tau_rayleigh + tau_cloud)

• ftau_cld_3d (ndarray of float) – Fraction of cloud extinction to total = tau_cloud/(tau_rayleigh + tau_cloud)

• ftau_ray_3d (ndarray of float) – Fraction of rayleigh extinction to total = tau_rayleigh/(tau_rayleigh + tau_cloud)

• dtau_og_3d (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT the delta eddington correction, if it was specified by user Dimensions=# layer by # wave

• tau_og_3d (ndarray of float) – This is the cumulative summed opacity WITHOUT the delta eddington correction, if it was specified by user Dimensions=# level by # wave

• w0_og_3d (ndarray of float) – Same as w0 but WITHOUT the delta eddington correction, if it was specified by user

• cosb_og_3d (ndarray of float) – Same as cosbar buth WITHOUT the delta eddington correction, if it was specified by user

• surf_reflect (float) – Surface reflectivity

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• ubar1 (ndarray of float) – matrix of cosine of the observer angles

• cos_theta (float) – Cosine of the phase angle of the planet

• F0PI (array) – Downward incident solar radiation

• single_phase (str) – Single scattering phase function, default is the two-term henyey-greenstein phase function

• multi_phase (str) – Multiple scattering phase function, defulat is N=2 Legendre polynomial

• frac_a (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_b (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_c (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C), Default is : 1 - gcosb^2

• constant_back (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of back scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• constant_forward (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of forward scatterer. Remember, the output of A & M code does not separate back and forward scattering.

Returns:

• intensity at the top of the atmosphere for all the different ubar1 and ubar2

• To Do

• —–

• - F0PI Solar flux shouldn’t always be 1.. Follow up to make sure that this isn’t a bad – hardwiring to solar, despite “relative albedo”

• - take delta eddington option out of fluxes and move it all over to optics

picaso.fluxes.get_reflected_SH(nlevel, nwno, numg, numt, dtau, tau, w0, cosb, ftau_cld, ftau_ray, f_deltaM, dtau_og, tau_og, w0_og, cosb_og, surf_reflect, ubar0, ubar1, cos_theta, F0PI, w_single_form, w_multi_form, psingle_form, w_single_rayleigh, w_multi_rayleigh, psingle_rayleigh, frac_a, frac_b, frac_c, constant_back, constant_forward, stream, b_top=0, flx=0, single_form=0)

Computes rooney fluxes given tau and everything is 3 dimensional. This is the exact same function as get_flux_geom_1d but is kept separately so we don’t have to do unecessary indexing for retrievals.

Parameters:

• nlevel (int) – Number of levels in the model

• nwno (int) – Number of wave points

• numg (int) – Number of Gauss angles

• numt (int) – Number of Chebyshev angles

• dtau (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• tau (ndarray of float) – This is the cumulative summed opacity WITHOUT D-Eddington Correction Dimensions=# level by # wave

• w0 (ndarray of float) – This is the single scattering albedo, from scattering, clouds, raman, etc WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• cosb (ndarray of float) – This is the asymmetry factor WITHOUT D-Eddington Correction Dimensions=# layer by # wave

• ftau_cld (ndarray of float) – Fraction of cloud extinction to total = tau_cloud/(tau_rayleigh + tau_cloud)

• ftau_ray (ndarray of float) – Fraction of rayleigh extinction to total = tau_rayleigh/(tau_rayleigh + tau_cloud)

• f_deltaM (ndarray of float) – Fractional scattering coefficient for delta-M calculation f_deltaM = 0 if delta_eddington=False

• dtau_og (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT the delta eddington correction, if it was specified by user Dimensions=# layer by # wave

• tau_og (ndarray of float) – This is the cumulative summed opacity WITHOUT the delta eddington correction, if it was specified by user Dimensions=# level by # wave

• w0_og (ndarray of float) – Same as w0 but WITHOUT the delta eddington correction, if it was specified by user

• cosb_og (ndarray of float) – Same as cosbar buth WITHOUT the delta eddington correction, if it was specified by user

• surf_reflect (float) – Surface reflectivity

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• ubar1 (ndarray of float) – matrix of cosine of the observer angles

• cos_theta (float) – Cosine of the phase angle of the planet

• F0PI (array) – Downward incident solar radiation

• w_single_form (str) – Single scattering phase function approximation for SH

• w_multi_form (str) – Multiple scattering phase function approximation for SH

• psingle_form (str) – Scattering phase function approximation for psingle in SH

• w_single_rayleigh (str) – Toggle rayleigh scattering on/off for single scattering in SH

• w_multi_rayleigh (str) – Toggle rayleigh scattering on/off for multi scattering in SH

• psingle_rayleigh (str) – Toggle rayleigh scattering on/off for psingle in SH

• frac_a (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_b (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C)

• frac_c (float) – (Optional), If using the TTHG phase function. Must specify the functional form for fraction of forward to back scattering (A + B * gcosb^C), Default is : 1 - gcosb^2

• constant_back (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of back scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• constant_forward (float) – (Optional), If using the TTHG phase function. Must specify the assymetry of forward scatterer. Remember, the output of A & M code does not separate back and forward scattering.

• stream (int) – Order of expansion of Legendre polynomials (2 or 4)

• b_top (float) – Upper boundary condition for incoming intensity

• flx (int) – Toggle calculation of layerwise fluxes (0 = do not calculate, 1 = calculate)

• single_form (int) – Toggle which version of p_single to use (0 = explicit, 1 = legendre)

Returns:

• intensity at the top of the atmosphere for all the different ubar1 and ubar2

• To Do

• —–

• - F0PI Solar flux shouldn’t always be 1.. Follow up to make sure that this isn’t a bad – hardwiring to solar, despite “relative albedo”

• - take delta eddington option out of fluxes and move it all over to optics

picaso.fluxes.get_thermal_1d(nlevel, wno, nwno, numg, numt, tlevel, dtau, w0, cosb, plevel, ubar1, surf_reflect, hard_surface, dwno, calc_type)

This function uses the source function method, which is outlined here : https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JD094iD13p16287

The result of this routine is the top of the atmosphere thermal flux as a function of gauss and chebychev points accross the disk.

Everything here is in CGS units:

Fluxes - erg/s/cm^3 Temperature - K Wave grid - cm-1 Pressure ; dyne/cm2

Reminder: Flux = pi * Intensity, so if you are trying to compare the result of this with a black body you will need to compare with pi * BB !

Parameters:

• nlevel (int) – Number of levels which occur at the grid points (not to be confused with layers which are mid points)

• wno (numpy.ndarray) – Wavenumber grid in inverse cm

• nwno (int) – Number of wavenumber points

• numg (int) – Number of gauss points (think longitude points)

• numt (int) – Number of chebychev points (think latitude points)

• tlevel (numpy.ndarray) – Temperature as a function of level (not layer)

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• cosb (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the asymmetry of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• plevel (numpy.ndarray) – Pressure for each level (not layer, which is midpoints). CGS units (dyne/cm2)

• ubar1 (numpy.ndarray) – This is a matrix of ng by nt. This describes the outgoing incident angles and is generally computed in picaso.disco

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• hard_surface (int) – 0 for no hard surface (e.g. Jupiter/Neptune), 1 for hard surface (terrestrial)

• dwno (int) – delta wno needed for climate

• calc_type (int) – 0 for spectrum model, 1 for climate solver

Returns:

Thermal flux in CGS units (erg/cm3/s) in a matrix that is numg x numt x nwno

Return type:

numpy.ndarray

picaso.fluxes.get_thermal_1d_deprecate(nlevel, wno, nwno, numg, numt, tlevel, dtau, w0, cosb, plevel, ubar1, surf_reflect, hard_surface, tridiagonal)

This function uses the source function method, which is outlined here : https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JD094iD13p16287

The result of this routine is the top of the atmosphere thermal flux as a function of gauss and chebychev points accross the disk. Everything here is in CGS units: Fluxes - erg/s/cm^3 Temperature - K Wave grid - cm-1 Pressure ; dyne/cm2 Reminder: Flux = pi * Intensity, so if you are trying to compare the result of this with a black body you will need to compare with pi * BB ! :param nlevel: Number of levels which occur at the grid points (not to be confused with layers which are

mid points)

Parameters:

• wno (numpy.ndarray) – Wavenumber grid in inverse cm

• nwno (int) – Number of wavenumber points

• numg (int) – Number of gauss points (think longitude points)

• numt (int) – Number of chebychev points (think latitude points)

• tlevel (numpy.ndarray) – Temperature as a function of level (not layer)

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• cosb (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the asymmetry of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• plevel (numpy.ndarray) – Pressure for each level (not layer, which is midpoints). CGS units (dyne/cm2)

• ubar1 (numpy.ndarray) – This is a matrix of ng by nt. This describes the outgoing incident angles and is generally computed in picaso.disco

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• hard_surface (int) – 0 for no hard surface (e.g. Jupiter/Neptune), 1 for hard surface (terrestrial)

• tridiagonal (int) – 0 for tridiagonal, 1 for pentadiagonal

Returns:

Thermal flux in CGS units (erg/cm3/s) in a matrix that is numg x numt x nwno

Return type:

numpy.ndarray

picaso.fluxes.get_thermal_1d_gfluxi_deprecate(nlevel, wno, nwno, numg, numt, tlevel, dtau, w0, cosb, plevel, ubar1, surf_reflect, ugauss_angles, ugauss_weights, tridiagonal, calc_type, dwno)

This function uses the source function method, which is outlined here : https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JD094iD13p16287

The result of this routine is the top of the atmosphere thermal flux as a function of gauss and chebychev points accross the disk. Everything here is in CGS units: Fluxes - erg/s/cm^3 Temperature - K Wave grid - cm-1 Pressure ; dyne/cm2 Reminder: Flux = pi * Intensity, so if you are trying to compare the result of this with a black body you will need to compare with pi * BB ! :param nlevel: Number of levels which occur at the grid points (not to be confused with layers which are

mid points)

Parameters:

• wno (numpy.ndarray) – Wavenumber grid in inverse cm

• nwno (int) – Number of wavenumber points

• numg (int) – Number of gauss points (think longitude points)

• numt (int) – Number of chebychev points (think latitude points)

• tlevel (numpy.ndarray) – Temperature as a function of level (not layer)

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• cosb (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the asymmetry of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• plevel (numpy.ndarray) – Pressure for each level (not layer, which is midpoints). CGS units (dyne/cm2)

• ubar1 (numpy.ndarray) – This is a matrix of ng by nt. This describes the outgoing incident angles and is generally computed in picaso.disco

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• tridiagonal (int) – 0 for tridiagonal, 1 for pentadiagonal

• calc_type (int) – 0 for outgoing flux at top level output, 1 for upward and downward layer and level flux outputs

Returns:

• numpy.ndarray – Thermal flux in CGS units (erg/cm3/s) in a matrix that is numg x numt x nwno if calc_type=0

• numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray – Thermal flux in CGS units (erg/cm3/s) in four matrices if calc_type=0: that is level downward flux and level upward flux numg x numt x nlevel x nwno and then layer downward flux and layer upward flux numg x numt x nlayer x nwno

picaso.fluxes.get_thermal_3d(nlevel, wno, nwno, numg, numt, tlevel_3d, dtau_3d, w0_3d, cosb_3d, plevel_3d, ubar1, surf_reflect, hard_surface)

This function uses the source function method, which is outlined here : https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JD094iD13p16287

The result of this routine is the top of the atmosphere thermal flux as a function of gauss and chebychev points accross the disk. Everything here is in CGS units: Fluxes - erg/s/cm^3 Temperature - K Wave grid - cm-1 Pressure ; dyne/cm2 Reminder: Flux = pi * Intensity, so if you are trying to compare the result of this with a black body you will need to compare with pi * BB ! :param nlevel: Number of levels which occur at the grid points (not to be confused with layers which are

mid points)

Parameters:

• wno (numpy.ndarray) – Wavenumber grid in inverse cm

• nwno (int) – Number of wavenumber points

• numg (int) – Number of gauss points (think longitude points)

• numt (int) – Number of chebychev points (think latitude points)

• tlevel_3d (numpy.ndarray) – Temperature as a function of level (not layer). This 3d matrix has dimensions [nlevel,ngangle,ntangle].

• dtau_3d (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth. This 4d matrix has dimensions [nlevel, nwave,ngangle,ntangle].

• w0_3d (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations. This 4d matrix has dimensions [nlevel, nwave,ngangle,ntangle].

• cosb_3d (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the asymmetry of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations. This 4d matrix has dimensions [nlevel, nwave,ngangle,ntangle].

• plevel (numpy.ndarray) – Pressure for each level (not layer, which is midpoints). CGS units (dyne/cm2)

• ubar1 (numpy.ndarray) – This is a matrix of ng by nt. This describes the outgoing incident angles and is generally computed in picaso.disco

Returns:

Thermal flux in CGS units (erg/cm3/s) in a matrix that is numg x numt x nwno

Return type:

numpy.ndarray

picaso.fluxes.get_thermal_SH(nlevel, wno, nwno, numg, numt, tlevel, dtau, tau, w0, cosb, dtau_og, tau_og, w0_og, w0_no_raman, cosb_og, plevel, ubar1, surf_reflect, stream, hard_surface, flx=0)

The result of this routine is the top of the atmosphere thermal intensity as a function of gauss and chebychev points accross the disk.

Everything here is in CGS units:

Fluxes - erg/s/cm^3 Temperature - K Wave grid - cm-1 Pressure ; dyne/cm2

Parameters:

• nlevel (int) – Number of levels which occur at the grid points (not to be confused with layers which are mid points)

• wno (numpy.ndarray) – Wavenumber grid in inverse cm

• nwno (int) – Number of wavenumber points

• numg (int) – Number of gauss points (think longitude points)

• numt (int) – Number of chebychev points (think latitude points)

• tlevel (numpy.ndarray) – Temperature as a function of level (not layer)

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• tau (numpy.ndarray) – This is a matrix of nlevel by nwave. This describes the cumulative optical depth.

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• cosb (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the asymmetry of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• dtau_og (ndarray of float) – This is the opacity contained within each individual layer (defined at midpoints of “levels”) WITHOUT the delta eddington correction, if it was specified by user Dimensions=# layer by # wave

• tau_og (ndarray of float) – This is the cumulative summed opacity WITHOUT the delta eddington correction, if it was specified by user Dimensions=# level by # wave

• w0_og (ndarray of float) – Same as w0 but WITHOUT the delta eddington correction, if it was specified by user

• cosb_og (ndarray of float) – Same as cosbar buth WITHOUT the delta eddington correction, if it was specified by user

• plevel (numpy.ndarray) – Pressure for each level (not layer, which is midpoints). CGS units (dyne/cm2)

• ubar1 (numpy.ndarray) – This is a matrix of ng by nt. This describes the outgoing incident angles and is generally computed in picaso.disco

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• stream (int) – Order of expansion of Legendre polynomials (2 or 4)

• hard_surface (int) – 0 for no hard surface (e.g. Jupiter/Neptune), 1 for hard surface (terrestrial)

• flx (int) – Toggle calculation of layerwise fluxes (0 = do not calculate, 1 = calculate)

Returns:

Thermal flux in CGS units (erg/cm3/s) in a matrix that is numg x numt x nwno

Return type:

numpy.ndarray

picaso.fluxes.get_transit_1d(z, dz, nlevel, nwno, rstar, mmw, k_b, amu, player, tlayer, colden, DTAU)

Routine to get the transmission spectrum :param z: Altitude in decreasing order (cm) :type z: float, array :param dz: length of each atmospheric layer :type dz: float, array :param nlevel: Number of levels :type nlevel: int :param nwno: Number of wavelength points :type nwno: int :param rstar: Radius of star (cm) :type rstar: float :param mmw: Mean molecular weight :type mmw: float, array :param k_b: Boltzman constant cgs :type k_b: float :param amu: Atomic mass units cgs :type amu: float :param player: Pressure at layers (dyn/cm2) :type player: float, array :param tlayer: Temperature at layers (K) :type tlayer: float, array :param colden: Column density conputed in atmsetup.get_column_density() :type colden: float, array :param DTAU: Matrix of summed tau opacities from optics. This is

TAUGAS + TAURAY + TAUCLD

Returns:

Rp**2 /Rs**2 as a function of wavelength

Return type:

array

Notes

[1]

Brown, Timothy M. “Transmission spectra as diagnostics of extrasolar giant planet atmospheres.” The Astrophysical Journal 553.2 (2001): 1006.

picaso.fluxes.get_transit_3d(nlevel, nwno, radius, gravity, rstar, mass, mmw, k_b, G, amu, p_reference, plevel, tlevel, player, tlayer, colden, DTAU)

Routine to get the 3D transmission spectrum

picaso.fluxes.legP(mu)

Generate array of Legendre polynomials

picaso.fluxes.locate(array, value)

Parameters:

• array (array) – Array to be searched.

• value (float) – Value to be searched for.

Returns:

location of nearest point by bisection method

Return type:

int

picaso.fluxes.mat_dot(A, B)

Matrix-vector dot product

picaso.fluxes.numba_cumsum(mat)

Function to compute cumsum along axis=0 to bypass numba not allowing kwargs in cumsum

picaso.fluxes.pent_diag_solve(l, A, B, C, D, E, F)

Pentadiagonal Matrix solver :param A: :type A: array or list :param B: :type B: array or list :param C: :type C: array or list :param D: :type D: array or list :param E: :type E: array or list :param F: :type F: array or list

Returns:

Solution, x

Return type:

array

picaso.fluxes.planck_cgs_deprecate(wave, T, dwave)

picaso.fluxes.planck_rad_deprecate(iw, T, dT, tmin, tmax, bb, y2, tp)

picaso.fluxes.set_bb_deprecate(wno, delta_wno, nwno, ntmps, dt, tmin, tmax)

Function to compute a grid of black bodies before the code runs. This allows us to interpolate on a blackbody instead of computing the planck function repetitively. This was done because historically computing the planck function was a bottleneck in speed.

Parameters:

• wno (array, float) – Wavenumber array cm-1

• delta_wno (array, float) – Wavenumber bins cm-1

• nwno (int) – Number of wavenumbers (len(wno))

• ntmps (int) – Number of temperature points to compute. Default number is set in config.json

• dt (float) – Spacing in temperature to compute. Default number is set in config.json

• tmin (float) – Minimum temperature to compute the grid

• tmax (float) – Maximum temperature to compute the grid

Returns:

• array – black body grid (CGS), number of temperatures x number of wavenumbers

• array – spline values for interpolation, number of temperatures x number of wavenumbers

• array – temperature grid

picaso.fluxes.setup_2_stream_fluxes(nlayer, nwno, w0, b_top, b_surface, surf_reflect, ubar0, dtau, tau, a, b, B0=0.0, B1=0.0, fluxes=0, calculation=0)

Setup up matrices to solve flux problem for spherical harmonics method.

Parameters:

• nlayer (int) – Number of layers

• nwno (int) – Number of wavenumber points

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• b_top (array) – The diffuse radiation into the model at the top of the atmosphere

• b_surface (array) – The diffuse radiation into the model at the bottom. Includes emission, reflection of the unattenuated portion of the direct beam

• b_surface_SH4 (array) – Second bottom BC for SH4 method.

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• tau (numpy.ndarray) – This is a matrix of nlevel by nwave. This describes the cumulative optical depth.

• a (numpy.ndarray) – Coefficients of matrix capturing legendre expansion of phase function for multiple scattering

• b (numpy.ndarray) – Coefficients of source vector capturing legendre expansion of phase function for single scattering

• B0 (numpy.ndarray) – Matrix of blackbodies

• B1 (numpy.ndarray) – Eqn (26) Toon 89

• fluxes (int) – Toggle calculation of layerwise fluxes (0 = do not calculate, 1 = calculate)

• calculation (int) – Toggle calculation method (1 = linear, 2 = exponential)

Returns:

Matrices and vectors used to calculate fluxes and intensities at each level

Return type:

numpy.ndarrays

picaso.fluxes.setup_4_stream_fluxes(nlayer, nwno, w0, b_top, b_surface, b_surface_SH4, surf_reflect, ubar0, dtau, tau, a, b, B0=0.0, B1=0.0, fluxes=0, calculation=0)

Setup up matrices to solve flux problem for spherical harmonics method.

Parameters:

• nlayer (int) – Number of layers

• nwno (int) – Number of wavenumber points

• w0 (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the single scattering albedo of the atmosphere. Note this is free of any Raman scattering or any d-eddington correction that is sometimes included in reflected light calculations.

• b_top (array) – The diffuse radiation into the model at the top of the atmosphere

• b_surface (array) – The diffuse radiation into the model at the bottom. Includes emission, reflection of the unattenuated portion of the direct beam

• b_surface_SH4 (array) – Second bottom BC for SH4 method.

• surf_reflect (numpy.ndarray) – Surface reflectivity as a function of wavenumber.

• ubar0 (ndarray of float) – matrix of cosine of the incident angle from geometric.json

• dtau (numpy.ndarray) – This is a matrix of nlayer by nwave. This describes the per layer optical depth.

• tau (numpy.ndarray) – This is a matrix of nlevel by nwave. This describes the cumulative optical depth.

• a (numpy.ndarray) – Coefficients of matrix capturing legendre expansion of phase function for multiple scattering

• b (numpy.ndarray) – Coefficients of source vector capturing legendre expansion of phase function for single scattering

• B0 (numpy.ndarray) – Matrix of blackbodies

• B1 (numpy.ndarray) – Eqn (26) Toon 89

• fluxes (int) – Toggle calculation of layerwise fluxes (0 = do not calculate, 1 = calculate)

• calculation (int) – Toggle calculation method (0 = reflected, 1 = thermal)

Returns:

Matrices and vectors used to calculate fluxes and intensities at each level

Return type:

numpy.ndarrays

picaso.fluxes.setup_pent_diag(nlayer, nwno, c_plus_up, c_minus_up, c_plus_down, c_minus_down, b_top, b_surface, surf_reflect, gama, dtau, exptrm_positive, exptrm_minus, g1, g2, exptrm, lamda)

Parameters:

• nlayer (int) – number of layers in the model

• nwno (int) – number of wavelength points

• c_plus_up (array) – c-plus evaluated at the top of the atmosphere

• c_minus_up (array) – c_minus evaluated at the top of the atmosphere

• c_plus_down (array) – c_plus evaluated at the bottom of the atmosphere

• c_minus_down (array) – c_minus evaluated at the bottom of the atmosphere

• b_top (array) – The diffuse radiation into the model at the top of the atmosphere

• b_surface (array) – The diffuse radiation into the model at the bottom. Includes emission, reflection of the unattenuated portion of the direct beam

• surf_reflect (array) – Surface reflectivity

• g1 (array) – table 1 toon et al 1989

• g2 (array) – table 1 toon et al 1989

• g3 (array) – table 1 toon et al 1989

• lamba (array) – Eqn 21 toon et al 1989

• gama (array) – Eqn 22 toon et al 1989

• dtau (array) – Opacity per layer

• exptrm_positive (array) – Eqn 44, expoential terms needed for tridiagonal rotated layered, clipped at 35

• exptrm_minus (array) – Eqn 44, expoential terms needed for tridiagonal rotated layered, clipped at 35

Returns:

coefficient of the positive exponential term

Return type:

array

picaso.fluxes.setup_tri_diag(nlayer, nwno, c_plus_up, c_minus_up, c_plus_down, c_minus_down, b_top, b_surface, surf_reflect, gama, dtau, exptrm_positive, exptrm_minus)

Before we can solve the tridiagonal matrix (See Toon+1989) section “SOLUTION OF THE TwO-STREAM EQUATIONS FOR MULTIPLE LAYERS”, we need to set up the coefficients. :param nlayer: number of layers in the model :type nlayer: int :param nwno: number of wavelength points :type nwno: int :param c_plus_up: c-plus evaluated at the top of the atmosphere :type c_plus_up: array :param c_minus_up: c_minus evaluated at the top of the atmosphere :type c_minus_up: array :param c_plus_down: c_plus evaluated at the bottom of the atmosphere :type c_plus_down: array :param c_minus_down: c_minus evaluated at the bottom of the atmosphere :type c_minus_down: array :param b_top: The diffuse radiation into the model at the top of the atmosphere :type b_top: array :param b_surface: The diffuse radiation into the model at the bottom. Includes emission, reflection

of the unattenuated portion of the direct beam

Parameters:

• surf_reflect (array) – Surface reflectivity

• g1 (array) – table 1 toon et al 1989

• g2 (array) – table 1 toon et al 1989

• g3 (array) – table 1 toon et al 1989

• lamba (array) – Eqn 21 toon et al 1989

• gama (array) – Eqn 22 toon et al 1989

• dtau (array) – Opacity per layer

• exptrm_positive (array) – Eqn 44, expoential terms needed for tridiagonal rotated layered, clipped at 35

• exptrm_minus (array) – Eqn 44, expoential terms needed for tridiagonal rotated layered, clipped at 35

Returns:

coefficient of the positive exponential term

Return type:

array

picaso.fluxes.slice_eq(array, lim, value)

Funciton to replace values with upper or lower limit

picaso.fluxes.slice_gt(array, lim)

Funciton to replace values with upper or lower limit

picaso.fluxes.slice_lt(array, lim)

Funciton to replace values with upper or lower limit

picaso.fluxes.slice_lt_cond(array, cond_array, cond, newval)

Funciton to replace values with upper or lower limit

picaso.fluxes.slice_lt_cond_arr(array, cond_array, cond, newarray)

Funciton to replace values with upper or lower limit

picaso.fluxes.slice_rav(array, lim)

Funciton to replace values with upper or lower limit

picaso.fluxes.solve_4_stream_banded(M, B, stream)

Solve the Spherical Harmonics Problem

Returns:

• intgrl_new (numpy.ndarray) – Integrated source function for source function technique

• flux (numpy.ndarray) – Upwards lux at bottom of atmosphere

• X (numpy.ndarray) – Coefficients of flux/intensity matrix problem

picaso.fluxes.spline_deprecate(x, y, n, yp0, ypn)

picaso.fluxes.tidal_flux(T_e, wave_in, nlevel, pressure, pm, hratio, col_den)

Computes Tidal Fluxes in all levels. Py of TIDALWAVE subroutine.

Parameters

T_efloat

Temperature (internal?)

wave_infloat

what is this?

nlevelint

# of levels

pressurearray

pressure array

pmfloat

Some pressure (?)

hratiofloat

Ratio of Scale Height over Chapman Scale Height

col_denarray

Column density array

Returns

Tidal Fluxes and DE/DM in ergs/g sec

picaso.fluxes.tri_diag_solve(l, a, b, c, d)

Tridiagonal Matrix Algorithm solver, a b c d can be NumPy array type or Python list type. refer to this wiki and to this explanation.

A, B, C and D refer to: .. math:: A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = D(I) This solver returns X. :param A: :type A: array or list :param B: :type B: array or list :param C: :type C: array or list :param C: :type C: array or list

Returns:

Solution, x

Return type:

array

picaso.fluxes.vec_dot(A, B)

Vector-vector dot product

picaso.io_utils module

picaso.io_utils.read_hdf(filename, requires, raise_except=False, **kwargs)

read in H5 format file and pull out specific table based on specific header information in the requires field. return None if file is not there

Parameters:

• filename (str) – name of H5 file to read

• requires (dict) – dictionary with keys and values to get correct table e.g. {‘wmin’:0.3, ‘wmax’:1, ‘R’:5000} for opacity grid

Returns:

d – data from H5 file decoded into python object

Return type:

python object

picaso.io_utils.read_json(filename, **kwargs)

read in a JSON format file. return None if the file is not there.

Parameters:

• filename (string) – name of the JSON file to read

• except (bool) – if true, raise exception if file is missing. if false, return empty dict.

Returns:

d – data from JSON file decoded into python object

Return type:

python object

picaso.justdoit module

picaso.justdoit.HJ_cld()

Function to get rough Jupiter Cloud model with fsed=3

picaso.justdoit.HJ_pt()

Function to get Jupiter’s PT profile

picaso.justdoit.HJ_pt_3d(as_xarray=False, add_kz=False, input_file='/data/reference_data/picaso/reference/base_cases/HJ_3d.pt')

Function to get Jupiter’s PT profile

Parameters:

• as_xarray (bool) – Returns as xarray, instead of dictionary

• add_kz (bool) – Returns kzz along with PT info

• input_file (str) – point to input file in the same format as mitgcm example file in base_cases/HJ_3d.pt

picaso.justdoit.OH_conc(temp, press, x_h2o, x_h2)

picaso.justdoit.SH_calculate_fluxes_options(printout=True)

Retrieve options for calculating layerwise fluxes

picaso.justdoit.SH_psingle_form_options(printout=True)

Retrieve options for direct scattering form approximation

picaso.justdoit.SH_rayleigh_options(printout=True)

Retrieve options for rayleigh scattering

picaso.justdoit.SH_scattering_options(printout=True)

Retrieve all the options for scattering radiation in SH

picaso.justdoit.all_planets()

Load all planets from https://exoplanetarchive.ipac.caltech.edu

picaso.justdoit.brown_dwarf_cld()

Function to get rough Brown Dwarf cloud model with Teff=1270 K M/H=1xSolar, C/O=1xSolar, fsed=1

picaso.justdoit.brown_dwarf_pt()

Function to get rough Brown Dwarf climate model with Teff=1270 K M/H=1xSolar, C/O=1xSolar, fsed=1

picaso.justdoit.check_units(unit)

picaso.justdoit.evolution_track(mass=1, age='all')

Plot or grab an effective temperature for a certain age and mass planet. :param mass: (Optional) Mass of planet, in Jupiter Masses. Only valid options = 1, 2, 4, 6, 8, 10,’all’

If another value is entered, it will find the nearest option.

Parameters:

age (float or str, optional) – (Optional) Age of planet, in years or ‘all’ to return full model

Returns:

• if age=None (data = {‘cold’:all_data_cold_start, ‘hot’:all_data_hot_start})

• if age=float (data = {‘cold’:data_at_age, ‘hot’:data_at_age})

• if plot=False (returns data)

• else (returns data, plot)

picaso.justdoit.find_press(at_tau, a, b, c)

picaso.justdoit.find_strat(mieff_dir, pressure, temp, dtdp, FOPI, nofczns, nstr, x_max_mult, t_table, p_table, grad, cp, opacityclass, grav, rfaci, rfacv, nlevel, tidal, tmin, tmax, dwni, bb, y2, tp, cloudy, cld_species, mh, fsed, flag_hack, save_profile, all_profiles, opd_cld_climate, g0_cld_climate, w0_cld_climate, flux_net_ir_layer, flux_plus_ir_attop, verbose=1)

Function iterating on the TP profile by calling tstart and changing opacities as well :param mieff_dir: path to directory with mieff files for virga :type mieff_dir: str :param it_max: Maximum iterations allowed in the inner no opa change loop :type it_max: int :param itmx: Maximum iterations allowed in the outer opa change loop :type itmx: int :param conv: :type conv: float :param convt: Convergence criteria , if max avg change in temp is less than this then outer loop converges :type convt: float :param nofczns: # of conv zones :type nofczns: int :param nstr: dimension of 20

NSTR vector describes state of the atmosphere: 0 is top layer 1 is top layer of top convective region 2 is bottom layer of top convective region 3 is top layer of lower radiative region 4 is top layer of lower convective region 5 is bottom layer of lower convective region

Parameters:

• xmaxmult –

• temp (array) – Guess temperatures to start with

• pressure (array) – Atmospheric pressure

• t_table (array) – Visible flux addition fraction

• nlevel (int) – # of levels

• temp – Guess Temperature array, dimension is nlevel

• pressure – Pressure array

• t_table – Tabulated Temperature array for convection calculations

• p_table (array) – Tabulated pressure array for convection calculations

• grad (array) – Tabulated grad array for convection calculations

• cp (array) – Tabulated cp array for convection calculations

• opacityclass (class) – Opacity class created with jdi.oppanection

• grav (float) – Gravity of planet in SI

• rfaci (float) – IR flux addition fraction

• rfacv (float) – Visible flux addition fraction

• nlevel – # of levels, not layers

• tidal (array) – Tidal Fluxes dimension = nlevel

• tmin (float) – Minimum allwed Temp in the profile

• tmax (float) – Maximum allowed Temp in the profile

• dwni (array) – Spectral interval corrections (dimension= nwvno)

• verbose (int) – If 0 nothing gets printed If 1 this prints everything out

Returns:

Temperature array and lapse ratio array if converged else Temperature array twice

Return type:

array

picaso.justdoit.find_strat_deq(mieff_dir, pressure, temp, dtdp, FOPI, nofczns, nstr, x_max_mult, t_table, p_table, grad, cp, opacityclass, grav, rfaci, rfacv, nlevel, tidal, tmin, tmax, dwni, bb, y2, tp, cloudy, cld_species, mh, fsed, flag_hack, quench_levels, kz, mmw, save_profile, all_profiles, self_consistent_kzz, save_kzz, all_kzz, opd_cld_climate, g0_cld_climate, w0_cld_climate, flux_net_ir_layer, flux_plus_ir_attop, photo_inputs_dict=None, on_fly=False, gases_fly=None, verbose=1)

Function iterating on the TP profile by calling tstart and changing opacities as well :param mieff_dir: path to directory with mieff files for virga :type mieff_dir: str :param it_max: Maximum iterations allowed in the inner no opa change loop :type it_max: int :param itmx: Maximum iterations allowed in the outer opa change loop :type itmx: int :param conv: :type conv: float :param convt: Convergence criteria , if max avg change in temp is less than this then outer loop converges :type convt: float :param nofczns: # of conv zones :type nofczns: int :param nstr: dimension of 20

NSTR vector describes state of the atmosphere: 0 is top layer 1 is top layer of top convective region 2 is bottom layer of top convective region 3 is top layer of lower radiative region 4 is top layer of lower convective region 5 is bottom layer of lower convective region

Parameters:

• xmaxmult –

• temp (array) – Guess temperatures to start with

• pressure (array) – Atmospheric pressure

• t_table (array) – Visible flux addition fraction

• nlevel (int) – # of levels

• temp – Guess Temperature array, dimension is nlevel

• pressure – Pressure array

• t_table – Tabulated Temperature array for convection calculations

• p_table (array) – Tabulated pressure array for convection calculations

• grad (array) – Tabulated grad array for convection calculations

• cp (array) – Tabulated cp array for convection calculations

• opacityclass (class) – Opacity class created with jdi.oppanection

• grav (float) – Gravity of planet in SI

• rfaci (float) – IR flux addition fraction

• rfacv (float) – Visible flux addition fraction

• nlevel – # of levels, not layers

• tidal (array) – Tidal Fluxes dimension = nlevel

• tmin (float) – Minimum allwed Temp in the profile

• tmax (float) – Maximum allowed Temp in the profile

• dwni (array) – Spectral interval corrections (dimension= nwvno)

• verbose (int) – If 1= prints out all output If 0= does not print anything out

Returns:

Temperature array and lapse ratio array if converged else Temperature array twice

Return type:

array

picaso.justdoit.get_contribution(bundle, opacityclass, at_tau=1, dimension='1d')

Currently top level program to run albedo code

Parameters:

• bundle (dict) – This input dict is built by loading the input = justdoit.load_inputs()

• opacityclass (class) – Opacity class from justdoit.opannection

• at_tau (float) – (Optional) Default = 1, This is to compute the pressure level at which cumulative opacity reaches at_tau. Usually users want to know when the cumulative opacity reaches a tau of 1.

• dimension (str) – (Optional) Default = ‘1d’. Currently only 1d is supported.

Returns:

• taus_by_species (dict) – Each dictionary entry is a nlayer x nwave that represents the per layer optical depth for that molecule. If you do not see a molecule that you have added as input, check to make sure it is propertly formatted (e.g. Sodium must be Na not NA, Titanium Oxide must be TiO not TIO)

• cumsum_taus (dict) – Each dictionary entry is a nlevel x nwave that represents the cumulative summed opacity for that molecule. If you do not see a molecule that you have added as input, check to make sure it is propertly formatted (e.g. Sodium must be Na not NA, Titanium Oxide must be TiO not TIO)

• at_pressure_array (dict) – Each dictionary entry is a nwave array that represents the pressure level where the cumulative opacity reaches the value specieid by the user through at_tau. If you do not see a molecule that you have added as input, check to make sure it is a molecule that you have added as input, check to make sure it is propertly formatted (e.g. Sodium must be Na not NA, Titanium Oxide must be TiO not TIO)

picaso.justdoit.get_targets()

Function to grab available targets using exoplanet archive data.

Return type:

Dataframe from Exoplanet Archive

picaso.justdoit.input_xarray(xr_usr, opacity, p_reference=10, calculation='planet')

This takes an input based on these standards and runs: -gravity -phase_angle -star -approx (p_reference=10) -atmosphere -clouds (if there are any)

Parameters:

• xr_usr (xarray) – xarray based on ERS formatting requirements

• opacity (justdoit.opannection) – opacity connection

• p_reference (float) – Default is to take from xarray reference pressure in bars

• calculation (str) – ‘planet’ or ‘browndwarf’

Example

case = jdi.input_xarray(xr_user) case.spectrum(opacity,calculation=’transit_depth’)

class picaso.justdoit.inputs(calculation='planet', climate=False)

Bases: object

Class to setup planet to run

Parameters:

• calculation (str) – (Optional) Controls planet or brown dwarf calculation. Default = ‘planet’. Other option is “browndwarf”.

• climate (bool) – (Optional) If true, this do a thorough iterative calculation to compute a temperature pressure profile.

phase_angle

Type:

set phase angle

gravity

Type:

set gravity

star

Type:

set stellar spec

atmosphere

Type:

set atmosphere composition and PT

clouds

Type:

set cloud profile

approx

Type:

set approximation

spectrum

Type:

create spectrum

TP_line_earth(P, Tsfc=294.0, Psfc=1.0, gam_trop=0.18, Ptrop=0.199, gam_strat=-0.045, Pstrat=0.001, nlevel=150)

Author: Mike R. Line Estimates Earth’s pressure-temperature profile. All default values have been tuned to semi reproduce Earth’s temperature pressure profile. :param P: Pressure array usually np.logspace(-6,2,nlevel) :type P: array :param Tsfc: Surface Temperature (K). Earth is 294 K :type Tsfc: float,optional :param Psfc: Surface Pressure (bar). Earth is 1 bar. :type Psfc: float ,optional :param gam_trop: Tropospheric dry lapse rate. Earth is ~0.18 :type gam_trop: float ,optional :param Ptrop: Tropospheric pressure. Earth 0.199 bar. :type Ptrop: float ,optional :param gam_strat: Stratospheric lapse rate. Earth is -0.045 :type gam_strat: float ,optional :param Pstrat: Stratospheric pressure (bars). Earth is 0.001 bar. / :type Pstrat: float ,optional :param nlevel: Number of grid levels :type nlevel: int ,optional

Returns:

Temperature array (K)Psfc

Return type:

array

T_eff(Teff=None)

Get Teff for climate run

Parameters:

T_eff (float) – (Optional) Effective temperature of Planet

add_pt(T, P)

Adds temperature pressure profile to atmosphere :param T: Temperature Array :type T: array :param P: Pressure Array :type P: array :param nlevel: # of atmospheric levels :type nlevel: int

Returns:

• T (numpy.array) – Temperature grid

• P (numpy.array) – Pressure grid

approx(single_phase='TTHG_ray', multi_phase='N=2', delta_eddington=True, raman='pollack', tthg_frac=[1, -1, 2], tthg_back=-0.5, tthg_forward=1, p_reference=1, rt_method='toon', stream=2, toon_coefficients='quadrature', single_form='explicit', calculate_fluxes='off', query='nearest_neighbor', w_single_form='TTHG', w_multi_form='TTHG', psingle_form='TTHG', w_single_rayleigh='on', w_multi_rayleigh='on', psingle_rayleigh='on', get_lvl_flux=False)

This function REsets all the default approximations in the code from what is in config file. This means that it will rewrite what is specified via config file defaults. It transforms the string specificatons into a number so that they can be used in numba nopython routines.

To see the str cases such as TTHG_ray users see all the options by using the function justdoit.single_phase_options or justdoit.multi_phase_options, etc.

single_phasestr

Single scattering phase function approximation

multi_phasestr

Multiple scattering phase function approximation

delta_eddingtonbool

Turns delta-eddington on and off

ramanstr

Uses various versions of raman scattering default is to use the pollack approximation

tthg_fraclist

Functional of forward to back scattering with the form of polynomial : tthg_frac[0] + tthg_frac[1]*g_b^tthg_frac[2] See eqn. 6 in picaso paper

tthg_backfloat

Back scattering asymmetry factor gf = g_bar*tthg_back

tthg_forwardfloat

Forward scattering asymmetry factor gb = g_bar * tthg_forward

p_referencefloat

Reference pressure (bars) This is an arbitrary pressure that corresponds do the user’s input of radius. Usually something “at depth” around 1-10 bars.

methodstr

Toon (‘toon’) or spherical harmonics (‘SH’).

streamint

Two stream or four stream (options are 2 or 4). For 4 stream need to set method=’SH’

toon_coefficients: str

Decide whether to use Quadrature (“quadrature”) or Eddington (“eddington”) schemes to define Toon coefficients in two-stream approximation (see Table 1 in Toon et al 1989)

single_formstr

form of the phase function can either be written as an ‘explicit’ henyey greinstein or it can be written as a ‘legendre’ expansion. Default is ‘explicit’

querystr

method to grab opacities. either “nearest_neighbor” or “interp” which interpolates based on 4 nearest neighbors. Default is nearest_neighbor which is significantly faster.

w_single_formstr

Single scattering phase function approximation for SH

w_multi_formstr

Multiple scattering phase function approximation for SH

psingle_formstr

Scattering phase function approximation for psingle in SH

w_single_rayleighstr

Toggle rayleigh scattering on/off for single scattering in SH

w_multi_rayleighstr

Toggle rayleigh scattering on/off for multi scattering in SH

psingle_rayleighstr

Toggle rayleigh scattering on/off for psingle in SH

get_lvl_fluxbool

This parameter returns the level by level and layer by layer fluxes in the full output Default is False

atmosphere(df=None, filename=None, exclude_mol=None, verbose=True, **pd_kwargs)

Builds a dataframe and makes sure that minimum necessary parameters have been suplied. Sets number of layers in model.

Parameters:

• df (pandas.DataFrame or dict) – (Optional) Dataframe with volume mixing ratios and pressure, temperature profile. Must contain pressure (bars) at least one molecule

• filename (str) – (Optional) Filename with pressure, temperature and volume mixing ratios. Must contain pressure at least one molecule

• exclude_mol (list of str) – (Optional) List of molecules to ignore from opacity. It will NOT change other aspects of the calculation like mean molecular weight. This should be used as exploratory ONLY. if you actually want to remove the contribution of a molecule entirely from your profile you should remove it from your input data frame.

• verbose (bool) – (Optional) prints out warnings. Default set to True

• pd_kwargs (kwargs) – Key word arguments for pd.read_csv to read in supplied atmosphere file

atmosphere_3d(ds, regrid=True, plot=True, iz_plot=0, verbose=True)

Checks your xarray input to make sure the necessary elements are included. If requested, it will regrid your output according to what you have specified in phase_angle() routine. If you have not requested a regrid, it will check to make sure that the latitude/longitude grid that you have specified in your xarray is the same one that you have set in the phase_angle() routine.

Parameters:

• ds (xarray.DataArray) – xarray input grid (see GCM 3D input tutorials)

• regrid (bool) – If True, this will auto regrid your data, based on the input to the phase_angle function you have supllied If False, it will skip regridding. However, this assumes that you have already regridded your data to the necessary gangles and tangles. PICASO will double check for you by comparing latitude/longitudes of what is in your xarray to what was computed in the phase_angle function.

• plot (bool) – If True, this will auto output a regridded plot

• iz_plot (int) – Altitude index to plot if it is requested

• verbose (bool) – If True, this will plot out messages, letting you know if your input data is being transformed

atmosphere_4d(ds=None, shift=None, plot=True, iz_plot=0, verbose=True, zero_point='night_transit')

Regrids xarray

Parameters:

• ds (xarray.DataArray) – xarray input grid (see GCM 3D input tutorials) Only optional if you have already defined your dataframe to self.inputs[‘atmosphere’][‘profile’]

• shift (array) – Degrees, for each orbital phase, picaso will rotate the longitude grid phase_i`+`shift_i. For example, for tidally locked planets, shift`=0 at all phase angles. Therefore, `shift must be input as an array of length n_phase, set by phase_angle() routine. Use plot=True to understand how your grid is being shifted.

• plot (bool) – If True, this will auto output a regridded plot

• iz_plot (bool) – pressure index to plot

• verbose (bool) – If True, this will plot out messages, letting you know if your input data is being transformed

• zero_point (str) – Is your zero point “night_transit”, or “secondary_eclipse” Default, “night_transit”

call_photochem(temp, pressure, logMH, cto, pressure_surf, mass, radius, kzz, tstop=10000000.0, filename=None, stfilename=None, network=None, network_ct=None, first=True, pc=None, user_psurf=True, user_psurf_add=3)

channon_grid_high(filename=None)

channon_grid_low(filename=None)

Interpolate from visscher grid

chem_interp(chem_grid)

Interpolates chemistry based on dataframe input of either 1460 or 1060 grid This particular function needs to have all molecules as columns as well as pressure and temperature

chemeq_3d(c_o=1.0, log_mh=0.0, n_cpu=1)

You must have already ran atmosphere_3d or pre-defined an xarray gcm before running this function.

This function will post-process sonora chemical equillibrium chemistry onto your 3D grid.

CURRENT options log m/h: 0.0, 0.5, 1.0, 1.5, 1.7, 2.0 C/O: 0.5X, 1.0X, 1.5X, 2.0X, 2.5X

Parameters:

• c_o (float,optional) – default = 1 (solar), options= 0.5X, 1.0X, 1.5X, 2.0X, 2.5X

• log_mh (float, optional) – default = 0 (solar), options = 0.0, 0.5, 1.0, 1.5, 1.7, 2.0

• n_cpu (int) – Number of cpu to use for parallelization of chemistry

chemeq_deprecate(CtoO, Met)

This interpolates from a precomputed grid of CEA runs (run by M.R. Line) :param CtoO: C to O ratio (solar = 0.55) :type CtoO: float :param Met: Metallicity relative to solar (solar = 1) :type Met: float

chemeq_visscher(c_o, log_mh)

Author of Data: Channon Visscher

Find nearest neighbor from visscher grid JUNE 2015 MODELS BASED ON 1060-POINT MARLEY GRID GRAPHITE ACTIVITY ADDED IN TEXT FILES (AFTER OCS) “ABUNDANCE” INDICATES CONDENSATION CONDITION (O OR 1) CURRENT GRID

FE/H: 0.0, 0.5, 1.0, 1.5, 1.7, 2.0

C/O: 0.5X, 1.0X, 1.5X, 2.0X, 2.5X

The solar carbon-to-oxygen ratio is calculated from Lodders (2010):

CARBON = 7.19E6 ATOMS OXYGEN = 1.57E7 ATOMS

This gives a “solar” C/O ratio of 0.458

The C/O ratio adjusted by keeping C + O = constant and adjusting the carbon-to-oxygen ratio by a factor relative to the solar value (i.e., a factor of “1” means 1x the solar value, i.e. a C/O ratio of 0.458).

This approach keeps the heavy-element-to-hydrogen ratio (Z/X) constant for a given [Fe/H]

These abundances are then multiplied by the metallicity factor (10**[Fe/H]) along with every other element in the model.

Parameters:

• co (int) – carbon to oxygen ratio relative to solar. Solar = 1

• log_mh (int) – metallicity (relative to solar) Will find the nearest value to 0.0, 0.5, 1.0, 1.5, 1.7, 2.0 Solar = 0

climate(opacityclass, save_all_profiles=False, as_dict=True, with_spec=False, save_all_kzz=False, diseq_chem=False, self_consistent_kzz=False, kz=None, on_fly=False, gases_fly=None, chemeq_first=True, verbose=True)

Top Function to run the Climate Model

Parameters:

• opacityclass (class) – Opacity class from justdoit.opannection

• save_all_profiles (bool) – If you want to save and return all iterations in the T(P) profile,True/False

• with_spec (bool) – Runs picaso spectrum at the end to get the full converged outputs, Default=False

• save_all_kzz (bool) – If you want to save and return all iterations in the kzz profile,True/False

• diseq_chem (bool) – If you want to run `on-the-fly’ mixing (takes longer),True/False

• self_consistent_kzz (bool) – If you want to run MLT in convective zones and Moses in the radiative zones

• kz (array) – Kzz input array if user wants constant or whatever input profile (cgs)

• verbose (bool) – If True, triggers prints throughout code

clouds(filename=None, g0=None, w0=None, opd=None, p=None, dp=None, df=None, **pd_kwargs)

Cloud specification for the model. Clouds are parameterized by a single scattering albedo (w0), an assymetry parameter (g0), and a total extinction per layer (opd). g0,w0, and opd are both wavelength and pressure dependent. Our cloud models come from eddysed. Their output look something like this where pressure is in bars and wavenumber is inverse cm. We will sort pressure and wavenumber before we reshape so the exact order doesn’t matter pressure wavenumber opd w0 g0 1. 1. … . . 1. 2. … . . 1. 3. … . . . . … . . . . … . . 1. M. … . . 2. 1. … . . . . … . . N. . … . . If you are creating your own file you have to make sure that you have a pressure (bars) and **wavenumber**(inverse cm) column. We will use this to make sure that your cloud and atmospheric profiles are on the same grid. **If there is no pressure or wavelength parameter we will assume that you are on the same grid as your atmospheric input, and on the eddysed wavelength grid! ** Users can also input their own fixed cloud parameters, by specifying a single value for g0,w0,opd and defining the thickness and location of the cloud. :param filename: (Optional) Filename with info on the wavelength and pressure-dependent single scattering

albedo, asymmetry factor, and total extinction per layer. Input associated pd_kwargs so that the resultant output has columns named : g0, w0 and opd. If you are not using the eddysed output, you will also need a wavenumber and pressure column in units of inverse cm, and bars.

Parameters:

• g0 (float, list of float) – (Optional) Asymmetry factor. Can be a single float for a single cloud. Or a list of floats for two different cloud layers

• w0 (list of float) – (Optional) Single Scattering Albedo. Can be a single float for a single cloud. Or a list of floats for two different cloud layers

• opd (list of float) – (Optional) Total Extinction in dp. Can be a single float for a single cloud. Or a list of floats for two different cloud layers

• p (list of float) – (Optional) Bottom location of cloud deck (LOG10 bars). Can be a single float for a single cloud. Or a list of floats for two different cloud layers

• dp (list of float) – (Optional) Total thickness cloud deck above p (LOG10 bars). Can be a single float for a single cloud or a list of floats for two different cloud layers Cloud will span 10**(np.log10(p-dp))

• df (pd.DataFrame, dict) – (Optional) Same as what would be included in the file, but in DataFrame or dict form

clouds_3d(ds, regrid=True, plot=True, iz_plot=0, iw_plot=0, verbose=True)

Checks your cloud xarray input to make sure the necessary elements are included. If requested, it will regrid your output according to what you have specified in phase_angle() routine. If you have not requested a regrid, it will check to make sure that the latitude/longitude grid that you have specified in your xarray is the same one that you have set in the phase_angle() routine.

Parameters:

• ds (xarray.DataArray) – xarray input grid (see cloud GCM 3D input tutorials)

• regrid (bool) – If True, this will auto regrid your data, based on the input to the phase_angle function you have supllied If False, it will skip regridding. However, this assumes that you have already regridded your data to the necessary gangles and tangles. PICASO will double check for you by comparing latitude/longitudes of what is in your xarray to what was computed in the phase_angle function.

• plot (bool) – If True, this will auto output a regridded plot

• iz_plot (int) – Altitude index to plot if a plot is requested

• iw_plot (int) – Wavelength index to plot if plot is requested

• verbose (bool) – If True, this will plot out messages, letting you know if your input data is being transformed

clouds_4d(ds=None, plot=True, iz_plot=0, iw_plot=0, verbose=True)

Regrids xarray

Parameters:

• ds (xarray.DataArray) – xarray input grid (see GCM 3D input tutorials) Only optional if you have already defined your dataframe to self.inputs[‘clouds’][‘profile’]

• plot (bool) – If True, this will auto output a regridded plot

• iz_plot (bool) – pressure index to plot

• iw_plot (bool) – wavelength index to plot

• verbose (bool) – If True, this will plot out messages, letting you know if your input data is being transformed

clouds_reset()

Reset cloud dict to zeros

effective_temp(teff=None)

Same as T_eff with different notation

Parameters:

teff (float) – (Optional) Effective temperature of Planet

gravity(gravity=None, gravity_unit=None, radius=None, radius_unit=None, mass=None, mass_unit=None)

Get gravity based on mass and radius, or gravity inputs

Parameters:

• gravity (float) – (Optional) Gravity of planet

• gravity_unit (astropy.unit) – (Optional) Unit of Gravity

• radius (float) – (Optional) radius of planet MUST be specified for thermal emission!

• radius_unit (astropy.unit) – (Optional) Unit of radius

• mass (float) – (Optional) mass of planet

• mass_unit (astropy.unit) – (Optional) Unit of mass

guillot_pt(Teq, T_int=100, logg1=-1, logKir=-1.5, alpha=0.5, nlevel=61, p_bottom=1.5, p_top=-6)

Creates temperature pressure profile given parameterization in Guillot 2010 TP profile called in fx() :param Teq: equilibrium temperature :type Teq: float :param T_int: Internal temperature, if low (100) currently set to 100 for everything :type T_int: float :param kv1: see parameterization Guillot 2010 (10.**(logg1+logKir)) :type kv1: float :param kv2: see parameterization Guillot 2010 (10.**(logg1+logKir)) :type kv2: float :param kth: see parameterization Guillot 2010 (10.**logKir) :type kth: float :param alpha: set to 0.5 :type alpha: float , optional :param nlevel: Number of atmospheric layers :type nlevel: int, optional :param p_bottom: Log pressure (bars) of the lower bound pressure :type p_bottom: float, optional :param p_top: Log pressure (bars) of the TOA :type p_top: float , optional

Returns:

• T (numpy.array) – Temperature grid

• P (numpy.array) – Pressure grid

inputs_climate(temp_guess=None, pressure=None, rfaci=1, nofczns=1, nstr=None, rfacv=None, m_planet=None, r_planet=None, cloudy=False, mh=None, CtoO=None, species=None, fsed=None, mieff_dir=None, photochem=False, photochem_file=None, photochem_stfile=None, photonetwork_file=None, photonetworkct_file=None, tstop=10000000.0, psurf=10)

Get Inputs for Climate run

Parameters:

• temp_guess (array) – Guess T(P) profile to begin with

• pressure (array) – Pressure Grid for climate code (this wont change on the fly)

• nstr (array) – NSTR vector describes state of the atmosphere: 0 is top layer [0] 1 is top layer of top convective region 2 is bottom layer of top convective region 3 is top layer of lower radiative region 4 is top layer of lower convective region 5 is bottom layer of lower convective region [nlayer-1]

• nofczns (integer) – Number of guessed Convective Zones. 1 or 2

• rfacv (float) – Fractional contribution of reflected light in net flux. =0 for no stellar irradition, =0.5 for full day-night heat redistribution =1 for dayside

• rfaci (float) – Default=1, Fractional contribution of thermal light in net flux Usually this is kept at one and then the redistribution is controlled via rfacv

• cloudy (bool) – Include Clouds or not (True or False)

• mh (string) – Metallicity string for 1060 grid, ‘+0.5’,’0.0’,’-0.5’.

• CtoO (string) – C/O ratio string for 1060 grid

• species (string) – Cloud species to be included if cloudy

• fsed (float) – Sedimentation Efficiency (f_sed) if cloudy

• mieff_dir (str) – path to directory with mieff files for virga

• photochem (bool) – Turns off (False) and on (True) Photochem

phase_angle(phase=0, num_gangle=10, num_tangle=1, symmetry=False, phase_grid=None, calculation=None)

Define phase angle and number of gauss and tchebychev angles to compute. Controls all geometry of the calculation. Computes latitude and longitudes. Computes cos theta and incoming and outgoing angles. Adds everything to class.

Please see geometry notebook for a deep dive into how these angles are incorporated into the calculation.

Typical numbers: - 1D Thermal: num_gangle = 10, num_tangle=1 - 1D Reflected light with zero phase : num_gangle = 10, num_tangle=1 - 1D Reflected light with non zero phase : num_gangle = 8, num_tangle=8 - 3D Thermal or Reflected : angles will depend on 3D features interested in resolving.

Parameters:

• phase (float,int) – Phase angle in radians

• num_gangle (int) – Number of Gauss angles to integrate over facets (Default is 10). This is defined as angles over the full sphere. So 10 as the default compute the RT with 5 points on the west hemisphere and 5 points on the west. Higher numbers will slow down code.

• num_tangle (int) – Number of Tchebyshev angles to integrate over facets (Default is 1, which automatically assumes symmetry). Must be even if considering symmetry.

• symmetry (bool) – Default is False. Note, if num_tangle=1 you are automatically considering symmetry. This will only compute the unique incoming angles so that the calculation isn’t doing repetetive models. E.g. if num_tangles=10, num_gangles=10. Instead of running a 10x10=100 FT calculations, it will only run 5x5 = 25 (or a quarter of the sphere).

• phase_grid (array) – This is ONLY for computing phase curves and ONLY can be used

• calculation (str) – Used for phase curve calculation. Needs to be specified in order to determine integrable angles correctly.

phase_curve(opacityclass, full_output=False, plot_opacity=False, n_cpu=1)

Run phase curve :param opacityclass: Opacity class from justdoit.opannection :type opacityclass: class :param full_output: (Optional) Default = False. Returns atmosphere class, which enables several

plotting capabilities.

Parameters:

n_cpu (int) – (Optional) Default = 1 (no parallelization). Number of cpu to parallelize calculation.

phase_curve_geometry(calculation, phase_grid, num_gangle=10, num_tangle=10)

Geometry setup for phase curve calculation. This computes all the necessary

Parameters:

• calculation (str) – ‘reflected’ or ‘thermal’

• phase_grid (float) – phase angle grid to compute phase curve in radians

• num_gangle (int) – number of gauss angles, equivalent to longitude

• num_tangle (int) – number of tchebychev angles, equivalent to latitude

premix_3d(opa, n_cpu=1)

You must have already ran atmosphere_3d or pre-defined an xarray gcm before running this function.

This function will post-process sonora chemical equillibrium chemistry onto your 3D grid.

CURRENT options log m/h: 0.0, 0.5, 1.0, 1.5, 1.7, 2.0 C/O: 0.5X, 1.0X, 1.5X, 2.0X, 2.5X

Parameters:

• c_o (float,optional) – default = 1 (solar), options= 0.5X, 1.0X, 1.5X, 2.0X, 2.5X

• log_mh (float, optional) – default = 0 (solar), options = 0.0, 0.5, 1.0, 1.5, 1.7, 2.0

• n_cpu (int) – Number of cpu to use for parallelization of chemistry

premix_atmosphere(opa, df=None, filename=None, **pd_kwargs)

Builds a dataframe and makes sure that minimum necessary parameters have been suplied. Sets number of layers in model. :param opa: Opacity class from opannection : RetrieveCks() :type opa: class :param df: (Optional) Dataframe with volume mixing ratios and pressure, temperature profile.

Must contain pressure (bars) at least one molecule

Parameters:

• filename (str) – (Optional) Filename with pressure, temperature and volume mixing ratios. Must contain pressure at least one molecule

• exclude_mol (list of str) – (Optional) List of molecules to ignore from file

• pd_kwargs (kwargs) – Key word arguments for pd.read_csv to read in supplied atmosphere file

premix_atmosphere_diseq(opa, quench_levels, t_mix=None, df=None, filename=None, **pd_kwargs)

Builds a dataframe and makes sure that minimum necessary parameters have been suplied. Sets number of layers in model. :param opa: Opacity class from opannection : RetrieveCks() :type opa: class :param df: (Optional) Dataframe with volume mixing ratios and pressure, temperature profile.

Must contain pressure (bars) at least one molecule

Parameters:

• filename (str) – (Optional) Filename with pressure, temperature and volume mixing ratios. Must contain pressure at least one molecule

• exclude_mol (list of str) – (Optional) List of molecules to ignore from file

• pd_kwargs (kwargs) – Key word arguments for pd.read_csv to read in supplied atmosphere file

setup_climate()

Turns off planet specific things, so program can run as usual

setup_nostar()

Turns off planet specific things, so program can run as usual

sonora(sonora_path, teff, chem='low')

This queries Sonora temperature profile that can be downloaded from profiles.tar on Zenodo:

• Bobcat Models: [profile.tar file](https://zenodo.org/record/1309035#.Xo5GbZNKjGJ)

• Elf OWL Models: [L Type Models](https://zenodo.org/records/10385987), [T Type Models](https://zenodo.org/records/10385821), [Y Type Models](https://zenodo.org/records/10381250)

Note gravity is not an input because it grabs gravity from self.

Parameters:

• sonora_path (str) – Path to the untarred profile.tar file from sonora grid

• teff (float) – teff to query. does not have to be exact (it will query one the nearest neighbor)

• chem (str) – Default = ‘low’. There are two sonora chemistry grids. the default is use low. There is sublety to this that is going to be explained at length in the sonora grid paper. Until then, ONLY use chem=’low’ unless you have reached out to one of the developers.

spectrum(opacityclass, calculation='reflected', dimension='1d', full_output=False, plot_opacity=False, as_dict=True)

Run Spectrum :param opacityclass: Opacity class from justdoit.opannection :type opacityclass: class :param calculation: Either ‘reflected’ or ‘thermal’ for reflected light or thermal emission.

If running a brown dwarf, this will automatically default to thermal

Parameters:

• dimension (str) – (Optional) Dimensions of the calculation. Default = ‘1d’. But ‘3d’ is also accepted. In order to run ‘3d’ calculations, user must build 3d input (see tutorials)

• full_output (bool) – (Optional) Default = False. Returns atmosphere class, which enables several plotting capabilities.

• plot_opacity (bool) – (Optional) Default = False, Creates pop up of the weighted opacity

• as_dict (bool) – (Optional) Default = True. If true, returns a condensed dictionary to the user. If false, returns the atmosphere class, which can be used for debugging. The class is clunky to navigate so if you are consiering navigating through this, ping one of the developers.

star(opannection, temp=None, metal=None, logg=None, radius=None, radius_unit=None, semi_major=None, semi_major_unit=None, database='ck04models', filename=None, w_unit=None, f_unit=None)

Get the stellar spectrum using pysynphot and interpolate onto a much finer grid than the planet grid.

Parameters:

• opannection (class picaso.RetrieveOpacities) – This is the opacity class and it’s needed to get the correct wave info and raman scattering cross sections

• temp (float) – (Optional) Teff of the stellar model if using the stellar database feature. Not needed for filename option.

• metal (float) – (Optional) Metallicity of the stellar model if using the stellar database feature. Not needed for filename option.

• logg (float) – (Optional) Logg cgs of the stellar model if using the stellar database feature. Not needed for filename option.

• radius (float) – (Optional) Radius of the star. Only needed as input if you want relative flux units (Fp/Fs)

• radius_unit (astropy.unit) – (Optional) Any astropy unit (e.g. radius_unit=astropy.unit.Unit(“R_sun”))

• semi_major (float) – (Optional) Semi major axis of the planet. Only needed to compute fp/fs for albedo calculations.

• semi_major_unit (astropy.unit) – (Optional) Any astropy unit (e.g. radius_unit=astropy.unit.Unit(“au”))

• database (str) – (Optional)The database to pull stellar spectrum from. See documentation for pysynphot. Most popular are ‘ck04models’, phoenix’ and

• filename (str) – (Optional) Upload your own stellar spectrum. File format = two column white space (wave, flux). Must specify w_unit and f_unit

• w_unit (str) – (Optional) Used for stellar file wave units. Needed for filename input. Pick: ‘um’, ‘nm’, ‘cm’, ‘hz’, or ‘Angs’

• f_unit (str) – (Optional) Used for stellar file flux units. Needed for filename input. Pick: ‘FLAM’ or ‘Jy’ or ‘erg/cm2/s/Hz’

surface_reflect(albedo, wavenumber, old_wavenumber=None)

Set atmospheric surface reflectivity. This preps the code to run a terrestrial planet. This will automatically change the run to “hardsurface”, which alters the lower boundary condition of the thermal_1d flux calculation. :param albedo: Set constant albedo for surface reflectivity :type albedo: float

virga(condensates, directory, fsed=1, b=1, eps=0.01, param='const', mh=1, mmw=2.2, kz_min=100000.0, sig=2, Teff=None, alpha_pressure=None, supsat=0, gas_mmr=None, do_virtual=False, verbose=True)

Runs virga cloud code based on the PT and Kzz profiles that have been added to inptus class. :param condensates: Condensates to run in cloud model :type condensates: str :param fsed: Sedimentation efficiency coefficient :type fsed: float :param b: Denominator of exponential in sedimentation efficiency (if param is ‘exp’) :type b: float :param eps: Minimum value of fsed function (if param=exp) :type eps: float :param param: fsed parameterisation

‘const’ (constant), ‘exp’ (exponential density derivation), ‘pow’ (power-law)

Parameters:

• mh (float) – Metallicity

• mmw (float) – Atmospheric mean molecular weight

• gas_mmr (dict) – Gas MMR as a dictionary for individual gases. This allows users to override virga’s chemistry. E.g. {‘SiO2’:1e-6}

• kz_min (float) – Minimum kzz value

• sig (float) – Width of the log normal distribution for the particle sizes

• Teff (float, optional) – Effective temperature. If None, Teff set to temperature at 1 bar

• alpha_pressure (float, optional) – Pressure at which we want fsed=alpha for variable fsed calculation. If None, pressure set to the top of the atmosphere

• do_virtual (bool) – Turn on and off the “virtual” cloud which is a cloud that forms below the pressure grid defined by the user.

• verbose (bool) – Turn off warnings

virga_3d(condensates, directory, fsed=1, mh=1, mmw=2.2, kz_min=100000.0, sig=2, n_cpu=1, verbose=True, smooth_kz=False, full_output=False)

Runs virga cloud code based on the PT and Kzz profiles that have been added to inptus class.

Parameters:

• condensates (str) – Condensates to run in cloud model

• fsed (float) – Sedimentation efficiency

• mh (float) – Metallicity

• mmw (float) – Atmospheric mean molecular weight

• n_cpu (int) – number cpu to parallelize

• verbose (bool) – Print statements to help user

• smooth_kz (bool) – If true, it uses the min_kz value and does a UnivariateSpline accross the kz values to smooth out the profile

• full_output (bool) – Returns full output of virga model run

picaso.justdoit.jupiter_cld()

Function to get rough Jupiter Cloud model with fsed=3

picaso.justdoit.jupiter_pt()

Function to get Jupiter’s PT profile

picaso.justdoit.load_planet(df, opacity, phase_angle=0, stellar_db='ck04models', verbose=False, **planet_kwargs)

Wrapper to simplify PICASO run. This really turns picaso into a black box. This was created specifically Sagan School tutorial. It grabs planet parameters from the user supplied df, then adds a parametric PT profile (Guillot et al. 2010). The user still has to add chemistry.

Parameters:

• df (pd.DataFrame) – This is single row from all_planets()

• opacity (np.array) – Opacity loaded from opannection

• phase_angle (float) – Observing phase angle (radians)

• stellar_db (str) – Stellar database to pull from. Default is ck04models but you can also use phoenix if you have those downloaded.

• verbose (bool , options) – Print out warnings

• planet_kwargs (dict) – List of parameters to supply NexSci database is values don’t exist

picaso.justdoit.multi_phase_options(printout=True)

Retrieve all the options for multiple scattering radiation

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

picaso.justdoit.output_xarray(df, picaso_class, add_output={}, savefile=None)

This function converts all picaso output to xarray which is easy to save It is recommended for reuse and reproducibility of models. the returned xarray can be used with input_xarray to rerun models. See Preservation notebook.

Parameters:

• df (dict) – This is the output of your spectrum and must include “full_output=True” For example, df = case.spectrum(opa, full_output=True) It can also be the output of your climate run. For example, out=case.climate(…, withSpec=True) If running climate you must run withSpec=True OR add output of case.spectrum to your output dictionary by doing (for example): >> out = case.climate(…, withSpec=False) >> #add any post processing steps you desire (e.g. case.atmosphere or case.clouds) >> df = case.spectrum(opa, calculation=’thermal’) >> out[‘spectrum_output’] = df

• picaso_class (picaso.inputs) – This is the original class that you made to do your run. For example, case=jdi.inputs();followed by case.atmosphere(), case.clouds(), etc

• add_output (dict) – These are any additional outputs you want included in your xarray metadata This dictionary has a very specific struture if you want it to work correctly. To see the structure you can run justdoit.standard_metadata() which will return an empty dictionary that you can fill.

• savefile (str) – Optional, if string it will save a xarray file for you

Returns:

this xarray dataset can be easily passed to justdoit.input_xarray to run a spectrum

Return type:

xarray.Dataset

Todo

• figure out why pandas index are being returned for pressure and wavelenth

• fix clouds wavenumber_layer which doesnt seem like it would work

• add clima inputs : teff (or tint), number of convective zones

• cvs_locs array should be more clear

picaso.justdoit.picaso(bundle, opacityclass, dimension='1d', calculation='reflected', full_output=False, plot_opacity=False, as_dict=True)

Currently top level program to run albedo code :param bundle: This input dict is built by loading the input = justdoit.load_inputs() :type bundle: dict :param opacityclass: Opacity class from justdoit.opannection :type opacityclass: class :param dimension: (Optional) Dimensions of the calculation. Default = ‘1d’. But ‘3d’ is also accepted.

In order to run ‘3d’ calculations, user must build 3d input (see tutorials)

Parameters:

• full_output (bool) – (Optional) Default = False. Returns atmosphere class, which enables several plotting capabilities.

• plot_opacity (bool) – (Optional) Default = False, Creates pop up of the weighted opacity

• as_dict (bool) – (Optional) Default = True. If true, returns a condensed dictionary to the user. If false, returns the atmosphere class, which can be used for debugging. The class is clunky to navigate so if you are consiering navigating through this, ping one of the developers.

Return type:

dictionary with albedos or fluxes or both (depending on what calculation type)

picaso.justdoit.profile(mieff_dir, it_max, itmx, conv, convt, nofczns, nstr, x_max_mult, temp, pressure, FOPI, t_table, p_table, grad, cp, opacityclass, grav, rfaci, rfacv, nlevel, tidal, tmin, tmax, dwni, bb, y2, tp, final, cloudy, cld_species, mh, fsed, flag_hack, save_profile, all_profiles, opd_cld_climate, g0_cld_climate, w0_cld_climate, flux_net_ir_layer=None, flux_plus_ir_attop=None, first_call_ever=False, verbose=True)

Function iterating on the TP profile by calling tstart and changing opacities as well :param mieff_dir: path to directory with mieff files for virga :type mieff_dir: str :param it_max: Maximum iterations allowed in the inner no opa change loop :type it_max: int :param itmx: Maximum iterations allowed in the outer opa change loop :type itmx: int :param conv: :type conv: float :param convt: Convergence criteria , if max avg change in temp is less than this then outer loop converges :type convt: float :param nofczns: # of conv zones :type nofczns: int :param nstr: dimension of 20

NSTR vector describes state of the atmosphere: 0 is top layer 1 is top layer of top convective region 2 is bottom layer of top convective region 3 is top layer of lower radiative region 4 is top layer of lower convective region 5 is bottom layer of lower convective region

Parameters:

• xmaxmult –

• temp (array) – Guess temperatures to start with

• pressure (array) – Atmospheric pressure

• t_table (array) – Visible flux addition fraction

• nlevel (int) – # of levels

• temp – Guess Temperature array, dimension is nlevel

• pressure – Pressure array

• t_table – Tabulated Temperature array for convection calculations

• p_table (array) – Tabulated pressure array for convection calculations

• grad (array) – Tabulated grad array for convection calculations

• cp (array) – Tabulated cp array for convection calculations

• opacityclass (class) – Opacity class created with jdi.oppanection

• grav (float) – Gravity of planet in SI

• rfaci (float) – IR flux addition fraction

• rfacv (float) – Visible flux addition fraction

• nlevel – # of levels, not layers

• tidal (array) – Tidal Fluxes dimension = nlevel

• tmin (float) – Minimum allwed Temp in the profile

• tmax (float) – Maximum allowed Temp in the profile

• dwni (array) – Spectral interval corrections (dimension= nwvno)

Returns:

Temperature array and lapse ratio array if converged else Temperature array twice

Return type:

array

picaso.justdoit.profile_deq(mieff_dir, it_max, itmx, conv, convt, nofczns, nstr, x_max_mult, temp, pressure, FOPI, t_table, p_table, grad, cp, opacityclass, grav, rfaci, rfacv, nlevel, tidal, tmin, tmax, dwni, bb, y2, tp, final, cloudy, cld_species, mh, fsed, flag_hack, quench_levels, kz, mmw, save_profile, all_profiles, self_consistent_kzz, save_kzz, all_kzz, opd_cld_climate, g0_cld_climate, w0_cld_climate, flux_net_ir_layer, flux_plus_ir_attop, photo_inputs_dict=None, on_fly=False, gases_fly=None, verbose=True)

Function iterating on the TP profile by calling tstart and changing opacities as well :param mieff_dir: path to directory with mieff files for virga :type mieff_dir: str :param it_max: Maximum iterations allowed in the inner no opa change loop :type it_max: int :param itmx: Maximum iterations allowed in the outer opa change loop :type itmx: int :param conv: :type conv: float :param convt: Convergence criteria , if max avg change in temp is less than this then outer loop converges :type convt: float :param nofczns: # of conv zones :type nofczns: int :param nstr: dimension of 20

NSTR vector describes state of the atmosphere: 0 is top layer 1 is top layer of top convective region 2 is bottom layer of top convective region 3 is top layer of lower radiative region 4 is top layer of lower convective region 5 is bottom layer of lower convective region

Parameters:

• xmaxmult –

• temp (array) – Guess temperatures to start with

• pressure (array) – Atmospheric pressure

• t_table (array) – Visible flux addition fraction

• nlevel (int) – # of levels

• temp – Guess Temperature array, dimension is nlevel

• pressure – Pressure array

• t_table – Tabulated Temperature array for convection calculations

• p_table (array) – Tabulated pressure array for convection calculations

• grad (array) – Tabulated grad array for convection calculations

• cp (array) – Tabulated cp array for convection calculations

• opacityclass (class) – Opacity class created with jdi.oppanection

• grav (float) – Gravity of planet in SI

• rfaci (float) – IR flux addition fraction

• rfacv (float) – Visible flux addition fraction

• nlevel – # of levels, not layers

• tidal (array) – Tidal Fluxes dimension = nlevel

• tmin (float) – Minimum allwed Temp in the profile

• tmax (float) – Maximum allowed Temp in the profile

• dwni (array) – Spectral interval corrections (dimension= nwvno)

• verbose (bool) – If True, prints out messages

Returns:

Temperature array and lapse ratio array if converged else Temperature array twice

Return type:

array

picaso.justdoit.query_options()

Retrieve options for querying opacities

picaso.justdoit.raman_options()

Retrieve options for raman scattering approximtions

picaso.justdoit.rt_methodology_options(printout=True)

Retrieve all the options for methodology

picaso.justdoit.single_phase_options(printout=True)

Retrieve all the options for direct radation

picaso.justdoit.standard_metadata()

picaso.justdoit.stream_options(printout=True)

Retrieve all the options for stream

picaso.justdoit.toon_phase_coefficients(printout=True)

Retrieve options for coefficients used in Toon calculation

picaso.justdoit.young_planets()

Load planets from ZJ’s paper

picaso.justplotit module

picaso.justplotit.all_optics_1d(full_output, wave_range, return_output=False, legend=None, ng=None, nt=None, colors=('#0072B2', '#E69F00', '#F0E442', '#009E73', '#56B4E9', '#D55E00', '#CC79A7', '#000000'), **kwargs)

Plots 1d profiles of optical depth per layer, single scattering, and asymmetry averaged over the user input wave_range.

Parameters:

• full_output (list,dict) – Dictonary of outputs or list of dicts

• wave_range (list) – min and max wavelength in microns

• return_output (bool) – Default is just to return a figure but you can also return all the 1d profiles

• legend (bool) – Default is none. Legend for each component of out

• ng (int) – gauss index

• nt (int) – chebychev intex

• **kwargs (keyword arguments) – Key word arguments will be supplied to each bokeh figure function

picaso.justplotit.animate_convergence(clima_out, picaso_bundle, opacity, calculation='thermal', wave_range=[0.3, 6], molecules=['H2O', 'CH4', 'CO', 'NH3'])

Function to animate climate convergence given all profiles that were computed throughout the climate run Animates from the first guess through to the final converged state. Runs spectra for each of those cases.

Parameters:

• clima_out (dict) – Output from, for example: clima_out = cl_run.climate(opacity_ck, save_all_profiles=True)

• picaso_bundle (<picaso.justdoit.inputs>) – This would be the equivalent of cl_run that you used to run climate_model(). This ensures that your spectra can be rerun

• opacity (<picaso.optics.RetrieveCKs>) – This is the jdi.opannection that you used as input in the run_climate_model

• wave_range (list) – sets the wavelngth range of the spectra

• molecules (list) – list of strings for what molecules to animate

Return type:

matplotlib.animation

picaso.justplotit.bin_errors(newx, oldx, dy)

Bin errors properly to account for reduction in noise

Parameters:

• newx (array) – New x axis (either micron or wavenumber)

• oldx (array) – Old x axis (either micron or wavenumber)

• dy (array) – Error bars

Returns:

new dy

Return type:

array

picaso.justplotit.brightness_temperature(out_dict, plot=True, R=None, with_guide=True)

Plots and returns brightness temperature

$T_{

m bright}=dfrac{a}{{lambda}logleft(dfrac{{b}}{F(lambda){lambda}^5}+1 ight)}$

where a = 1.43877735$ imes$10$^{-2}$ m.K and b = 11.91042952$ imes$10$^{-17}$ m$^4$kg/s$^3$

out_dictdict

output of bundle.spectrum(opa,full_output=True) or a dictionary with {“wavenumber”: np.array, “thermal”: np.array of standard output in erg/cm2/s/cm}

plotbool

If true creates and returns a plot

Rfloat

If not None, rebins the brightness temperature

with_guidebool

Plots points from the minimum and maximum pt values

picaso.justplotit.cloud(full_output)

Plotting the cloud input from picaso.

The plot itselfs creates maps of the wavelength dependent single scattering albedo and cloud opacity as a function of altitude.

Parameters:

full_output –

Return type:

A row of two bokeh plots with the single scattering and optical depth map

picaso.justplotit.create_heat_map(data, rayleigh=True, extend=False, plot_height=300, plot_width=300, font_size='12px')

picaso.justplotit.disco(full_output, wavelength, calculation='reflected')

Plot disco ball with facets. Bokeh is not good with 3D things. So this is in matplotlib

Parameters:

• full_output (class) – Full output from picaso

• wavelength (list) – Where to plot 3d facets. Can input as many wavelengths as wanted. Must be a list, must be in microns.

• calculation (str, optional) – Default is to plot ‘reflected’ light but can also switch to ‘thermal’ if it has been computed

picaso.justplotit.explore(df, key)

Function to explore a dictionary that is THREE levels deep and return the data sitting at the end of key.

Parameters:

df (dict) – Dictionary that you want to search through.

Examples

Consider a dictionary that has df[‘layer’][‘cloud’][‘w0’] = [0,0,0]

>>>explore(df, ‘w0’) [0,0,0]

picaso.justplotit.find_nearest_1d(array, value)

picaso.justplotit.find_nearest_2d(array, value, axis=1)

picaso.justplotit.find_nearest_old(array, value)

picaso.justplotit.flux_at_top(full_output, plot_bb=True, R=None, pressures=[0.1, 0.01, 0.001], ng=None, nt=None, **kwargs)

Routine to plot the OLR with overlaying black bodies.

Flux units are CGS = erg/s/cm^3

Parameters:

• full_output (dict) – full dictionary output with {‘wavenumber’,’thermal’,’full_output’}

• plot_bb (bool , optional) – Default is to plot black bodies for three pressures specified by pressures

• R (float) – New constant R to bin to

• pressures (list, optional) – Default is a list of three pressures (in bars) = [1e-1,1e-2,1e-3]

• ng (int) – Used for 3D calculations to select point on the sphere (equivalent to longitude point)

• nt (int) – Used for 3D calculations to select point on the sphere (equivalent to latitude point)

• **kwargs (dict) – Any key word argument for bokeh.figure()

picaso.justplotit.heatmap_taus(out, R=0)

Plots a heatmap of the tau fields (taugas, taucld, tauray)

Parameters:

• out (dict) – full_ouput dictionary

• R (int) – Resolution to bin to (if zero, no binning)

picaso.justplotit.lon_lat_to_cartesian(lon_r, lat_r, R=1)

calculates lon, lat coordinates of a point on a sphere with radius R

picaso.justplotit.map(full_output, pressure=[0.1], plot='temperature', wavelength=None, igauss=0)

Plot disco ball with facets. Bokeh is not good with 3D things. So this is in matplotlib

Parameters:

• full_output (class) – Full output from picaso

• pressure (list) – What pressure (in bars) to make the map on. Note: this assumes that the pressure grid is the same for each ng,nt point accross the grid.

• plot (str, optional) – Default is to plot ‘temperature’ map but can also switch to any 3d [nlayer, nlong, nlat] or 4d [nlayer, nwave, nlong, nlat] output in full_output. You can check what is available to plot by printing: print(full_output[‘layer’].keys(). If you are plotting something that is ALSO wavelength dependent you have to also supply a single wavelength.

• wavelength – This allows users to plot maps of things that are wavelength dependent, like taugas and taucld.

• float – This allows users to plot maps of things that are wavelength dependent, like taugas and taucld.

• optional – This allows users to plot maps of things that are wavelength dependent, like taugas and taucld.

• igauss (int) – Gauss point to plot if using ktables this can be greater than 0 up to ngauss-1. Otherwise, This must be zero for monochromatic opacities.

picaso.justplotit.mean_regrid(x, y, newx=None, R=None)

Rebin the spectrum at a minimum R or on a fixed grid!

Parameters:

• x (array) – Wavenumbers

• y (array) – Anything (e.g. albedo, flux)

• newx (array) – new array to regrid on.

• R (float) – create grid with constant R

• binwidths (array) – bin widths centered at x

Return type:

final x, and final y

picaso.justplotit.mixing_ratio(full_output, limit=50, ng=None, nt=None, plot_type='bokeh', molecules=None, **kwargs)

Returns plot of mixing ratios

Parameters:

• full_output (class) – picaso.atmsetup.ATMSETUP

• limit (int) – Limits the number of curves to 20. Will plot the top 20 molecules with highest max(abundance). Limit must be >=3.

• molecules (list) – Directly specify which molecule to plot

• plot_type (str) – Options are “bokeh” or “matplotlib”. Default is bokeh

• **kwargs (dict) – Any key word argument for bokeh.figure()

picaso.justplotit.molecule_contribution(contribution_out, opa, min_pressure=4.5, R=100, **kwargs)

Function to plot & graph the Tau~1 Pressure (bars) of various elements

Parameters:

• contribution_out (dict) – contribution_out from jdi.get_contribution. This function will grab contribution_out[‘tau_p_surface’] Pressure vs. wavelength optical depth surface at tau specified by user in get_contribution function (user input for at_tau)

• opa (picaso.opannection) – Picaso opacity connection to get wavelength

• min_pressure (float, int) – Minimum pressure contribution in bars for molecules you want to plot. Ignores all molecules that are optically thick higher than min_pressure (bars)

• R (int) – Resolution defined as lambda/dlambda

• Outputs –

• ------- –

• figure (bokeh.plotting.figure.Figure) – Shows a default graph of Tau 1 Surface of various molecules and a graph based on user input based on their parameters

picaso.justplotit.numba_cumsum(mat)

Function to compute cumsum along axis=0 to bypass numba not allowing kwargs in cumsum

picaso.justplotit.phase_curve(allout, to_plot, collapse=None, R=100, palette=('#5e4fa2', '#3288bd', '#66c2a5', '#abdda4', '#e6f598', '#ffffbf', '#fee08b', '#fdae61', '#f46d43', '#d53e4f', '#9e0142'), verbose=True, **kwargs)

Plots phase curves

Parameters:

• allouts (dict) – picaso allouts element that comes from jdi.phase_curve

• to_plot (str) – either fpfs_reflected, fpfs_thermal, or thermal, or albedo

• collapse (str or float or list of float) – Allowable options to collapse wavelength axis: - ‘np.mean’ or np.sum - float or list of float: wavelength(s) in microns (will find the nearest value to this wavelength). Must be in wavenumber range.

• R (float) – Resolution to regrid before finding nearest wavelength element

• palette (list) – list of hex from bokeh or other palette

• verbose (bool) – Print out low level warnings

• kwargs (dict) – Bokeh plotting kwargs for bokeh.Figure

picaso.justplotit.phase_snaps(allout, x='longitude', y='pressure', z='temperature', palette='RdBu_r', y_log=True, x_log=False, z_log=False, col_wrap=3, collapse='np.mean', igauss=0)

Parameters:

• x (str) – What to plot on the x axis options = (‘longitude’ or ‘latitude’ or ‘pressure’)

• y (str) – What to plot on the y axis (‘longitude’ or ‘latitude’ or ‘pressure’)

• z (str) – What to plot in the heatmap (‘temperature’,’taugas’,’taucld’,’tauray’,’w0’,’g0’,’opd’)

• y_log (bool) – Makes y axis log

• x_log (bool) – Makes x axis log

• z_log (bool) – Makes z axis log (colorbar)

• palette (str) – Color pallete

• col_wrap (int) – Column wrap, determines number of columns to split runs into

• collapse (str or int) – Collapse lets us know how to collapse the axis, not used. For instance, if plotting x=longitude, and y=pressure, with collapse=mean, it will take an average along the latitude axis. If collapse=0, it will take the 0th latitude point. Allowed collapse functions = np.mean, np.median, np.min, np.max

• igauss (int) – If using k-coeff gauss points, this can be changed to get different gauss quadrature points.

picaso.justplotit.photon_attenuation(full_output, at_tau=0.5, return_output=False, igauss=0, **kwargs)

Plot breakdown of gas opacity, cloud opacity, Rayleigh scattering opacity at a specified pressure level.

Parameters:

• full_output (class) – picaso.atmsetup.ATMSETUP

• at_tau (float) – Opacity at which to plot the cumulative opacity. Default = 0.5.

• return_output (bool) – Return photon attenuation plot values

• igauss (int) – Gauss angle to plot if using correlated-k method. If not, should always be 0.

• **kwargs (dict) – Any key word argument for bokeh.plotting.figure()

Returns:

• if return_output=False (bokeh plot)

• else (bokeh plot,wave,at_pressures_gas,at_pressures_cld,at_pressures_ray)

picaso.justplotit.plot_cld_input(nwno, nlayer, filename=None, df=None, pressure=None, wavelength=None, **pd_kwargs)

This function was created to investigate CLD input file for PICASO.

The plot itselfs creates maps of the wavelength dependent single scattering albedo and cloud opacity and assymetry parameter as a function of altitude.

Parameters:

• nwno (int) – Number of wavenumber points. For runs from Ackerman & Marley, this will always be 196.

• nlayer (int) – Should be one less than the number of levels in your pressure temperature grid. Cloud opacity is assigned for slabs.

• file (str , optional) – (Optional)Path to cloud input file

• df (str) – (Optional)Dataframe of cloud input file

• wavelength (array , optional) – (Optional) this allows you to reset the tick marks to wavelengths instead of indicies

• pressure (array, optional) – (Optional) this allows you to reset the tick marks to pressure instead of indicies

• pd_kwargs (kwargs) – Pandas key word arguments for pandas.read_csv

Return type:

Three bokeh plots with the single scattering, optical depth, and assymetry maps

picaso.justplotit.plot_errorbar(x, y, e, plot=None, point_kwargs={}, error_kwargs={}, plot_type='bokeh', plot_kwargs={})

Plot only symmetric y error bars in bokeh plot

Parameters:

• x (array) – x data

• y (array) – y data

• e (array) – +- error for y which will be distributed as y+e, y-e on data point

• plot (bokeh.figure, optional) – Bokeh figure to add error bars to

• plot_type (str, optional) – type of plot (either bokeh or matplotlib)

• point_kwargs (dict) – formatting for circles

• error_kwargs (dict) – formatting for error bar lines

• plot_kwargs (dict) – plot attriutes for bokeh figure

picaso.justplotit.plot_evolution(evo, y='Teff', **kwargs)

Plot evolution of tracks. Requires input from justdoit:

evo = justdoit.evolution_track(mass=’all’,age=’all’)

Parameters:

• evo (dict) – Output from the function justdoit.evolution_track(mass=’all’,age=’all’)

• y (str) – What to plot on the y axis. Can be anything in the pandas that as the mass attached to the end. E.g. “Teff” is an option because there exists “Teff1Mj”. But, age_years is not an option as it is not a function of mass. Current options : [logL, Teff, grav_cgs]

picaso.justplotit.plot_format(df)

Function to reformat plots

picaso.justplotit.plot_multierror(x, y, plot, dx_up=0, dx_low=0, dy_up=0, dy_low=0, point_kwargs={}, error_kwargs={})

Plot non-symmetric x and y error bars in bokeh plot

Parameters:

• x (array) – x data

• y (array) – y data

• dx_up (array or int or float) – upper error bar to be distributed as x + dx_up

• dx_low (array or int or float) – lower error bar to be distributed as x + dx_low

• dy_up (array or int or float) – upper error bar to be distributed as y + dy_up

• dy_low (array or int or float) – lower error bar to be distributed as y + dy_low

• plot (bokeh.figure) – Bokeh figure to add error bars to

• point_kwargs (dict) – formatting for circles

• error_kwargs (dict) – formatting for error bar lines

picaso.justplotit.pt(full_output, ng=None, nt=None, **kwargs)

Returns plot of pressure temperature profile

Parameters:

• full_output (class) – picaso.atmsetup.ATMSETUP

• **kwargs (dict) – Any key word argument for bokeh.figure()

picaso.justplotit.pt_adiabat(clima_out, input_class, plot=True)

Plot the PT profile with the adiabat

Parameters:

• clima_out (dict) – Full output from picaso

• input_class (justdoit.input) – PICASO input class (e.g. the result of jdi.inputs())

Return type:

adiabat, dTdP, pressure

picaso.justplotit.rt_heatmap(data, figure_kwargs={}, cmap_kwargs={})

This creates a heat map of the dataframe created from picaso.test.thermal_sh_test() or picaso.test.dlugach_test() The tutorial 10c_AnalyzingApproximationsThermal.ipynb and 10b_AnalyzingApproximationsReflectedLightSH.ipynb shows it in use.

Parameters:

• data (pandas.DataFrame) – Technically, this can be any dataframe, but this function was specifically built to replicate Figure 9 in Batalha et al. 2019 and Figure 6 from Rooney et al. 2023. Part II Thermal. It is assumed that the index of the dataframe is asymmetry and the columns are single scattering albedo.

• figure_kwargs (dict) – keywords for bokeh.plotting.figure function

• cmap_kwargs (dict) – keywords for bokeh.models.LinearColorMapper funtion. Most common parameters to input are palette (list or tuple of colors), high (upper lim of color map), and low (lower lim of color map)

picaso.justplotit.spectrum(xarray, yarray, legend=None, wno_to_micron=True, palette=('#0072B2', '#E69F00', '#F0E442', '#009E73', '#56B4E9', '#D55E00', '#CC79A7', '#000000'), muted_alpha=0.2, **kwargs)

Plot formated albedo spectrum

Parameters:

• xarray (float array, list of arrays) – wavenumber or micron

• yarray (float array, list of arrays) – albedo or fluxes

• legend (list of str , optional) – legends for plotting

• wno_to_micron (bool , optional) – Converts wavenumber to micron

• palette (list,optional) – List of colors for lines. Default only has 8 colors so if you input more lines, you must give a different pallete

• muted_alpha (float) – number 0-1 to indicate how muted you want the click functionaity

• **kwargs (dict) – Any key word argument for bokeh.plotting.figure()

Return type:

bokeh plot

picaso.justplotit.spectrum_hires(wno, alb, legend=None, **kwargs)

Plot formated albedo spectrum

Parameters:

• wno (float array, list of arrays) – wavenumber

• alb (float array, list of arrays) – albedo

• legend (list of str) – legends for plotting

• **kwargs (dict) – Any key word argument for hv.opts

Return type:

bokeh plot

picaso.justplotit.taumap(full_output, at_tau=1, wavelength=1, igauss=0)

Plot breakdown of gas opacity, cloud opacity, Rayleigh scattering opacity at a specified pressure level.

Parameters:

• full_output (class) – full_output from dictionary picaso output

• at_tau (float) – Opacity at which to plot the cumulative opacity. Default = 0.5.

• igauss (int) – Gauss point to plot if using ktables this can be greater than 0 up to ngauss-1. Otherwise, This must be zero for monochromatic opacities.

• **kwargs (dict) – Any key word argument for bokeh.plotting.figure()

Return type:

bokeh plot

picaso.justplotit.thermal_contribution(full_output, tau_max=1.0, R=100, **kwargs)

Computer the contribution function from https://doi.org/10.3847/1538-4357/aadd9e equation 4 Note the equation in the paper is missing the - sign in the exponent :param full_output: full dictionary output with {‘wavenumber’,’thermal’,’full_output’} :type full_output: dict :param tau_max: Maximum tau to consider as opaque, largely here to help prevent NaNs and the colorbar from being weird. :type tau_max: float, optional :param **kwargs: Any key word argument for pcolormesh :type **kwargs: dict

picaso.justplotit.transmission_contribution(full_output, R=None, **kwargs)

Compute the transmission contribution function from

https://petitradtrans.readthedocs.io/en/latest/content/notebooks/analysis.html#Transmission-contribution-functions

Parameters:

• full_output (dict) – full dictionary output with {‘wavenumber’,’thermal’,’full_output’}

• R (float, optional) – Resolution for rebinning

• **kwargs (dict) – Any key word argument for pcolormesh

Returns:

• fig – matplotlib figure

• ax – matplotlib ax

• um – micron array

• CF_bin – contribution function

picaso.opacity_factory module

picaso.opacity_factory.adapt_array(arr)

picaso.opacity_factory.build_skeleton(db_f)

This functionb builds a skeleton sqlite3 database with three tables: 1) header 2) molecular 3) continuum

picaso.opacity_factory.compute_ck_molecular(molecule, og_directory, wv_file_name=None, order=4, gfrac=0.95, dir_kark_ch4=None, alkali_dir=None, min_wavelength=None, max_wavelength=None, R=None, verbose=True)

Function to generate correlated-K tables for each individual gas

Parameters:

• molecule (str) – Name of molecule

• og_directory (str) – Directory of all the cross sections that include folders e.g. “H2O”, “CH4”

• alkali_dir (str) – Alakalis directory

• wv_file_name (str) – (optional) file name with wavelength. First column wavenumber, second column delta wavenumber Must supply this OR combo of min, max wavelength and R

• order (int) – (Optional) Gauss Legendre order which by default is set to 4, with the double gauss method

• gfrac (int) – (Optional) Double-gauss method of splitting up the gauss points, by default we use 0.95

• dir_kark_ch4 (str) – (Optional) directory of additional karkochka methane

• alkali_dir – (Optional) Directory of alkalis or “individual_file” which is the method to use for Roxana Lupus data.

• verbose (bool) – (Optional) prints out status of which p,t, point the code is at

picaso.opacity_factory.continuum_avail(db_file)

picaso.opacity_factory.convert_array(text)

picaso.opacity_factory.create_grid(min_wavelength, max_wavelength, constant_R)

Simple function to create a wavelength grid defined with a constant R.

Parameters:

• min_wavelength (float) – Minimum wavelength in microns

• max_wavelength (float) – Maximum wavelength in microns

• constant_R (float) – Constant R spacing

Return type:

wavenumber grid defined at constant Resolution

picaso.opacity_factory.create_grid_minR(min_wavelength, max_wavelength, minimum_R)

Simple function to create a wavelength grid defined with a minimum R. This does not create a “constant” resolution grid. Rather, it defines the R at the minimum R so that the wavelength grid satisfies all_Rs>R

picaso.opacity_factory.delete_molecule(mol, db_filename)

Delete a record from the database

Parameters:

• mol (str) – molecule to delete

• db_filename (str) – database to delete record from

picaso.opacity_factory.find_nearest(array, value)

picaso.opacity_factory.find_nearest_1d(array, value)

picaso.opacity_factory.fit_linsky(t, wno, va=3)

Adds H2-H2 opacity from Linsky (1969) and Lenzuni et al. (1991) to values that were set to -33 as place holders.

Parameters:

• t (float) – Temperature

• wno (numpy.array) – wave number

• va (int or float) – (Optional) 1,2 or 3 (depending on what overtone to compute

Return type:

H2-H2 absorption in cm-1 amagat-2

picaso.opacity_factory.func_read_gas(molecule, og_directory, wv_file_name=None, order=4, gfrac=0.95, dir_kark_ch4=None, alkali_dir=None, min_wavelength=None, max_wavelength=None, R=None, verbose=True)

Function to generate correlated-K tables for each individual gas

Parameters:

• molecule (str) – Name of molecule

• og_directory (str) – Directory of all the cross sections that include folders e.g. “H2O”, “CH4”

• alkali_dir (str) – Alakalis directory

• wv_file_name (str) – (optional) file name with wavelength. First column wavenumber, second column delta wavenumber Must supply this OR combo of min, max wavelength and R

• order (int) – (Optional) Gauss Legendre order which by default is set to 4, with the double gauss method

• gfrac (int) – (Optional) Double-gauss method of splitting up the gauss points, by default we use 0.95

• dir_kark_ch4 (str) – (Optional) directory of additional karkochka methane

• alkali_dir – (Optional) Directory of alkalis or “individual_file” which is the method to use for Roxana Lupus data.

• verbose (bool) – (Optional) prints out status of which p,t, point the code is at

picaso.opacity_factory.g_w_2gauss(order=4, gfrac=0.95)

Gets gauss and weight poitns for correlated-k function Specifically for the 2-point gauss method

Parameters:

• order (int) – Gauss order

• gfrac (float) – Adjustable parameter which sets the division in the values of G which will be sampled by the first and second set of Gauss points.

Return type:

gauss, weights

picaso.opacity_factory.get_continuum(db_file, species, temperature)

Grab continuum opacity from sqlite database

db_filestr

sqlite3 database filename

specieslist of str

Single species name. you can run continuum_avail to see what species are present.

temperaturelist of float

Temperature to grab. can grab available temperatures from continuum_avail

picaso.opacity_factory.get_h2minus(t, new_wno)

This results gives the H2 minus opacity, needed for temperatures greater than 600 K. K L Bell 1980 J. Phys. B: At. Mol. Phys. 13 1859, Table 1 theta=5040/T(K)

The result is given in cm4/dyn, which is why will will multiply by nh2*ne*k*T where: nh2: number of h2 molecules/cm3 ne: number of electrons/cm3 This will happen when we go to sum opacities. For now returns will be in cm4/dyn

Parameters:

• t (numpy.ndarray) – temperature in K

• wno (numpy.ndarray) – wave number cm-1

Return type:

H2 minus opacity in units of cm4/dyn

picaso.opacity_factory.get_hminusbf(wno)

H- bound free opacity, which is only dependent on wavelength. From John 1988 http://adsabs.harvard.edu/abs/1988A%26A…193..189J

Parameters:

wno (numpy.ndarray) – Wavenumber in cm-1

Returns:

Absorption coefficient in units of cm2

Return type:

array of floats

picaso.opacity_factory.get_hminusff(t, wno)

H- free free opacity, which is both wavelength and temperature dependent. From Bell & Berrington (1987) Also includes factor for simulated emission

Parameters:

• t (float) – Temperature in K

• wno (numpy.ndarray) – Wavenumber in cm-1

Returns:

Gives cross section in cm^5

Return type:

array of float

picaso.opacity_factory.get_kark_CH4(kark_file, new_wno, T, dset)

picaso.opacity_factory.get_kark_CH4_noTdependence(kark_dir, new_wave, temperature)

Files from kark+2010 paper

Return type:

opacity in cm2/species

picaso.opacity_factory.get_molecular(db_file, species, temperature, pressure)

Grab molecular opacity from sqlite database

db_filestr

sqlite3 database filename

specieslist of str

Single species name. you can run molecular_avail to see what species are present.

temperaturelist of flaot,int

Temperature to grab. This will grab the nearest pair neighbor so you dont have to be exact. This should be used in conjuction with pressure. For example, if you want one temperature at two pressures, then this should be input as t = [100,100] (and pressure would be p=[0.1,1] for instance)

pressurelist of flaot,int

Pressure to grab. This will grab the nearest pair neighbor so you dont have to be exact. This should be used in conjuction with pressure. For example, if you want one temperature at two pressures, then this should be input as t = [100,100] (and pressure would be p=[0.1,1] for instance)

picaso.opacity_factory.get_optical_o3(file_o3, new_wvno_grid)

This hacked ozone is from here: http://satellite.mpic.de/spectral_atlas/cross_sections/Ozone/O3.spc

picaso.opacity_factory.get_original_data(original_file, colnames, new_db, overwrite=False)

The continuum factory takes the CIA opacity file and adds in extra sources of opacity from other references to fill in empty bands. It assumes that the original file is structured as following,with the first column wavenumber and the rest molecules: 1000 198 75. 0.0 -33.0000 -33.0000 -33.0000 -33.0000 -33.0000 20.0 -7.4572 -7.4518 -6.8038 -6.0928 -5.9806 40.0 -6.9547 -6.9765 -6.6322 -5.7934 -5.4823 … … … etc

Where 1000 corresponds to the number of wavelengths, 198 corresponds to the number of temperatures and 75 is the first temperature. If this structure changes, we will need to restructure this top level __init__ function that reads it in the original opacity file.

Parameters:

• original_file (str) – Filepath that points to original opacity file (see description above)

• colnames (list) – defines the sources of opacity in the original file. For the example file above, colnames would be [‘wno’,’H2H2’,’H2He’,’H2H’,’H2CH4’,’H2N2’]

• new_db (str) – New database name

• overwrite (bool) – Default is set to False as to not overwrite any existing files. This parameter controls overwriting cia database

picaso.opacity_factory.get_wvno_grid(filename, min_wavelength=None, max_wavelength=None, R=None)

picaso.opacity_factory.h2h2_overtone(t, wno)

Add in special CIA h2h2 band at 0.8 microns

Parameters:

• t (float) – Temperature

• wno (numpy.array) – wave number

Return type:

H2-H2 absorption in cm-1 amagat-2

picaso.opacity_factory.insert(cur, conn, mol, T, opacity)

Insert into

picaso.opacity_factory.insert_hitran_cia(original_file, molname, new_db, new_wno)

This continuum function takes a HTIRAN CIA file and adds it as an extra source of opacity

Parameters:

• original_file (str) – Filepath that points to HITRAN file (e.g. H2-CH4_eq_2011.cia)

• new_db (str) – New database name

• new_wno (numpy.ndarray, list) – wavenumber grid to interpolate onto (units of inverse cm)

• overwrite (bool) – Default is set to False as to not overwrite any existing files. This parameter controls overwriting cia database

picaso.opacity_factory.insert_molecular_1060(molecule, min_wavelength, max_wavelength, new_R, og_directory, new_db, dir_kark_ch4=None, dir_optical_o3=None)

Function to resample 1060 grid data onto lower resolution grid. The general procedure in this function is to interpolate original 1060 data onto very high resolution grid (R=1e6). Then, determine number of bins to take given input ‘new_R’. The final opacity grid will be : original_opacity[::BINS]

NOTE: From several tests “new_R” should be at least 100x higher than the ultimate planet spectrum you want to bin down to.

Parameters:

• molecule (str) – Name of molecule (should match a directory with 1060 files)

• min_wavelength (float) – Minimum wavelength in database in units of micron

• max_wavelength (float) – Maximum wavelength in database in units of micron

• new_R (float) – New R to regrid to. If new_R=None, it will retain the original 1e6 million resolution of the lbl grid.

picaso.opacity_factory.insert_molecular_1460(molecule, min_wavelength, max_wavelength, og_directory, new_db, new_R=None, new_dwno=None, old_R=1000000.0, old_dwno=0.0035, alkali_dir='alkalis', dir_kark_ch4=None, dir_optical_o3=None, insert_direct=False)

DEVELOPER USE ONLY. Function to resample 1060 grid data onto lower resolution grid. The general procedure in this function is to interpolate original 1060 data onto very high resolution grid (R=1e6). Then, determine number of bins to take given input ‘new_R’. The final opacity grid will be : original_opacity[::BINS]

NOTE: From several tests “new_R” should be at least 100x higher than the ultimate planet spectrum you want to bin down to.

Parameters:

• molecule (str) – Name of molecule (should match a directory with 1060 files)

• min_wavelength (float) – Minimum wavelength in database in units of micron

• max_wavelength (float) – Maximum wavelength in database in units of micron

• new_R (int,float , optional) – Optional, new R to regrid to. This will create a new wavelength solution that is constant in R. The other option (new_dwno is to do constant wavenumber bin) This must be smaller than old_R, which is default to 1e6.

• old_R (int,float, optional) – This is set to match the approximate cross sections that are computed LBL. If you do not want any resampling, then you need to set new_R=old_R

• new_dwno (float, optional) – Optional, new constant wavenumber bin to create wavelength solution. This must be bigger than old_dwno.

• old_dwno (float, optional) – This is set to 0.0035 to match the smallest wavenumber than we compute in our grid. If you do not want any resampling, then you need to set new_dwno=old_dwno

• og_directory (str) – Directory of all the cross sections that include folders e.g. “H2O”, “CH4”

• alkali_dir (str) – Alakalis directory

• new_db (str) – New database name

• dir_kark_ch4 (str) – Karkoschka methane to hack in

• dir_optical_o3 (str) – optical ozone to hack in

• insert_direct (bool) – this ignores all other inputs except min and max wavelength and inserts opacities into the database diretly as is

picaso.opacity_factory.insert_molecular_1460_old(molecule, min_wavelength, max_wavelength, new_R, og_directory, new_db, dir_kark_ch4=None, dir_optical_o3=None)

Function to resample Ehsan’s 1460 grid data onto lower resolution grid, 1060 grid. The general procedure in this function is to interpolate original 1060 data onto very high resolution grid (R=1e6). Then, determine number of bins to take given input ‘new_R’. The final opacity grid will be : original_opacity[::BINS]

NOTE: From several tests “new_R” should be at least 100x higher than the ultimate planet spectrum you want to bin down to.

Parameters:

• molecule (str) – Name of molecule (should match a directory with 1060 files)

• min_wavelength (float) – Minimum wavelength in database in units of micron

• max_wavelength (float) – Maximum wavelength in database in units of micron

• new_R (float) – New R to regrid to. If new_R=None, it will retain the original 1e6 million resolution of the lbl grid.

picaso.opacity_factory.insert_wno_grid(wno_grid, cur, con)

Inserts basics into the header file. This puts all the units to the database.

picaso.opacity_factory.listdir(path)

picaso.opacity_factory.molecular_avail(db_file)

picaso.opacity_factory.open_local(db_f)

Code needed to open up local database, interpret arrays from bytes and return cursor

picaso.opacity_factory.regrid(x, y, newx=None, R=None, statistic='mean')

Rebin the spectrum at a minimum R or on a fixed grid

Parameters:

• x (array) – Wavenumbers

• y (array) – Anything (e.g. albedo, flux)

• newx (array) – new array to regrid on.

• R (float) – create grid with constant R

Return type:

final x, and final y

picaso.opacity_factory.restruct_continuum(original_file, colnames, new_wno, new_db, overwrite)

The continuum factory takes the CIA opacity file and adds in extra sources of opacity from other references to fill in empty bands. It assumes that the original file is structured as following,with the first column wavenumber and the rest molecules: 1000 198 75. 0.0 -33.0000 -33.0000 -33.0000 -33.0000 -33.0000 20.0 -7.4572 -7.4518 -6.8038 -6.0928 -5.9806 40.0 -6.9547 -6.9765 -6.6322 -5.7934 -5.4823 … … … etc

Where 1000 corresponds to the number of wavelengths, 198 corresponds to the number of temperatures and 75 is the first temperature. If this structure changes, we will need to restructure this top level __init__ function that reads it in the original opacity file.

Parameters:

• original_file (str) – Filepath that points to original opacity file (see description above)

• colnames (list) – defines the sources of opacity in the original file. For the example file above, colnames would be [‘wno’,’H2H2’,’H2He’,’H2H’,’H2CH4’,’H2N2’]

• new_wno (numpy.ndarray, list) – wavenumber grid to interpolate onto (units of inverse cm)

• new_db (str) – New database name

• overwrite (bool) – Default is set to False as to not overwrite any existing files. This parameter controls overwriting cia database

picaso.opacity_factory.restructure_opacity(new_db, ntemp, temperatures, molecules, og_opacity, old_wno, new_wno)

picaso.opacity_factory.vectorize_rebin_mean(bins, Fp)

picaso.opacity_factory.vectorize_rebin_median(bins, Fp)

picaso.opacity_factory.vresample_and_insert_molecular(molecule, min_wavelength, max_wavelength, new_R, og_directory, new_db)

This function is identical to resample_and_insert_molecular but it was written to determine if taking the median was better than taking every BIN’th data point. It uses a special vectorize function to speed up the median. Results are about the same but this is much slower to run. So dont bother using this unless there is something specific to try.

picaso.optics module

class picaso.optics.RetrieveCKs(ck_dir, cont_dir, wave_range=None, deq=False, on_fly=False, gases_fly=None)

Bases: object

This will be the class to retrieve correlated-k tables from the database. Right now this is in beta mode and is retrieving the direct heritage 196 grid files. :param ck_dir: Directory of the pre-mixed correlated k table data. Check that you are pointing to a

directory with the files “ascii_data” and “full_abunds”.

Parameters:

• cont_dir (str) – This should be in the references directory on Github. This has all the continuum opacity that you need (e.g. H2H2). It also has some extra cross sections that you may want to add in like Karkoschka methane and the optical ozone band.

• wave_range (list) – NOT FUNCTIONAL YET. Wavelength range to compuate in the format [min micron, max micron]

adapt_array()

needed to interpret bytes to array

convert_array(text)

needed to interpret bytes to array

get_available_continuum()

Get the pressures and temperatures that are available for the continuum and molecules

get_available_rayleigh()

get_continuum(atmosphere)

get_gauss_pts_661_1460_deprecate()

NOT needed anymore, should be replaced with 1460 function Function to read the legacy data of the 1060 grid computed by Roxana Lupu. .. note:

This function is **highly** sensitive to the file format. You cannot edit the ascii file and then
run this function. Each specific line is accounted for.

get_legacy_data_1060(wave_range, deq=False)

Function to read the legacy data of the 1060 grid computed by Roxana Lupu.

Note

This function is highly sensitive to the file format. You cannot edit the ascii file and then run this function. Each specific line is accounted for.

get_legacy_data_1460(wave_range)

Function to read the legacy data of the 1060 grid computed by Roxana Lupu.

Note

This function is highly sensitive to the file format. You cannot edit the ascii file and then run this function. Each specific line is accounted for.

get_mixing_indices(atmosphere)

Takes in atmosphere profile and returns an array which is nlayer by ngauss by nwno

get_new_wvno_grid_661()

get_opacities(atmosphere, exclude_mol=1)

atmosphereclass

picaso atmosphere class

exclude_molint

Not yet functional for CK option since they are premixed. For individual CK molecules, this will ignore the optical contribution from one molecule.

get_opacities_deq(bundle, atmosphere)

get_opacities_deq_onfly(bundle, atmosphere, gases_fly=None)

get_pre_mix_ck(atmosphere)

Takes in atmosphere profile and returns an array which is nlayer by ngauss by nwno

get_pre_mix_ck_nearest(atmosphere)

Takes in atmosphere profile and returns an array which is nlayer by ngauss by nwno

get_pre_mix_ck_sm(atmosphere)

Takes in atmosphere profile and returns an array which is nlayer by ngauss by nwno

load_kcoeff_arrays(path)

This loads and returns the kappa tables from Theodora’s bin_mol files will have this very hardcoded. use this func only once before run begins

load_kcoeff_arrays_first(path, gases_fly)

This loads and returns the kappa tables from opacity_factory.compute_ck_molecular() Currently we have a single file for each molecule

Parameters:

• path (str) – Path to the individual ck files

• gases_fly (bool) – Specifies what gasses to mix on the fly

mix_my_opacities(bundle, atmosphere)

Top Function to perform “on-the-fly” mixing and then interpolating of 5 opacity sources from Amundsen et al. (2017)

mix_my_opacities_gasesfly(bundle, atmosphere, gases_fly)

Top Function to perform “on-the-fly” mixing and then interpolating of 5 opacity sources from Amundsen et al. (2017)

open_local()

Code needed to open up local database, interpret arrays from bytes and return cursor

run_cia_spline()

run_cia_spline_661()

class picaso.optics.RetrieveOpacities(db_filename, raman_data, wave_range=None, location='local', resample=1)

Bases: object

This will be the class that will retrieve the opacities from the sqlite3 database. Right now we are just employing nearest neighbors to grab the respective opacities.

Parameters:

• db_filename (str) – Opacity file name

• raman_data (str) – file name of the raman cross section data

• wave_range (list of float) – Wavelength range (in microns) to run. Default None grabs all the full wavelength range available

• location (str) – Default = local. This is the only option currently available since we dont have AWS or other services enabled.

• resample (int) – Default =1 (no resampling)

raman_db

Database of Raman cross sections

Type:

pandas.DataFrame

wno

Wavenumber grid from opacity database (CGS cm-1)

Type:

array

wave

Wavelength grid from opacity database (micron)

Type:

array

nwno

Number of wavenumber points in the grid

Type:

int

cia_temps

CIA temps available from database

Type:

array

molecules

Molecules availalbe from database

Type:

array

pt_pairs

pressure-temperature pairs available from grid (our grid is usually computed on a fixed 1060 PT grid)

Type:

array

molecular_opa

Molecular Opacity dictionary that contains all the queried information necessary to do calculation

Type:

dict

continuum_opa

Continuum opacity dictionary that contains all the queried information necessary for the calculation

Type:

dict

rayleigh_opa

Rayleigh opacity dictionary for all the sources we have. This is computed as cm^2/species when the users runs opannection. We technically could just store this in the DB as well. However, because its so fast, it doesn’t take a lot of time, and because its non-temperature dependent, we only have to do it once.

Type:

dict

db_connect()

Inherents functionality of open_local. Eventually it will inherent functionality of “open_non_locals” (e.g. cloud or other sources)

open_local()

Opens sqlite3 connection and enables array interpretation from bytes

get_available_data()

Gets opacity database attributes such as available T-Ps , availalbe molecules. etc.

get_available_rayleigh()

Gets rayleigh cross sections when opannection is opened.

get_opacities()

Gets opacities after user specifies atmospheric profile (e.g. full PT and Composition)

adapt_array()

needed to interpret bytes to array

compute_stellar_shits(wno_star, flux_star)

Takes the hires stellar spectrum and computes Raman shifts

Parameters:

• wno_star (array) – hires wavelength grid of the star in wave number

• flux_star (array) – Flux of the star in whatever units. Doesn’t matter for this since it’s returning a ratio

Returns:

array of the shifted stellar spec divided by the unshifted stellar spec on the model wave number grid

Return type:

matrix

convert_array(text)

needed to interpret bytes to array

find_needed_pts(tlayer, player)

Function to identify the PT pairs needed given a set of temperature and pressures. When the opacity connection is established the PT pairs are read in and added to the class.

get_available_data(wave_range, resample)

Get the pressures and temperatures that are available for the continuum and molecules

get_available_rayleigh()

get_opacities(atmosphere, exclude_mol=1)

Get’s opacities using the atmosphere class using interpolation for molecular, but not continuum. Continuum opacity is grabbed via nearest neighbor methodology.

get_opacities_nearest(atmosphere, exclude_mol=1)

Get’s opacities using the atmosphere class

open_local()

Code needed to open up local database, interpret arrays from bytes and return cursor

preload_opacities(molecules, p_range, t_range)

Function that pre-loads opacities for certain molecules and p-t range

Parameters:

• molecules (list) – list of molecule names

• p_range (list) – list of floats in bars e.g [1e-2,1]

• t_range (list) – List of floats in kelvin e.g [300,1000]

picaso.optics.adapt_array(arr)

needed to interpret bytes to array

picaso.optics.bin_star(wno_new, wno_old, Fp)

Takes average of group of points using uniform tophat

Parameters:

• wno_new (numpy.array) – inverse cm grid to bin to (cm-1)

• wno_old (numpy.array) – inverse cm grid (old grid)

• Fp (numpy.array) – transmission spectra, which is on wno grid

picaso.optics.compute_opacity(atmosphere, opacityclass, ngauss=1, stream=2, delta_eddington=True, test_mode=False, raman=0, plot_opacity=False, full_output=False, return_mode=False)

Returns total optical depth per slab layer including molecular opacity, continuum opacity. It should automatically select the molecules needed

Parameters:

• atmosphere (class ATMSETUP) – This inherets the class from atmsetup.py

• opacityclass (class opacity) – This inherets the class from optics.py. It is done this way so that the opacity db doesnt have to be reloaded in a retrieval

• ngauss (int) – Number of gauss angles if using correlated-k. If using monochromatic opacites, ngauss should one.

• stream (int) – Number of streams (only affects detal eddington approximation)

• delta_eddington (bool) – (Optional) Default=True, With Delta-eddington on, it incorporates the forward peak contribution by adjusting optical properties such that the fraction of scattered energy in the forward direction is removed from the scattering parameters

• raman (int) – (Optional) Default =0 which corresponds to oklopcic+2018 raman scattering cross sections. Other options include 1 for original pollack approximation for a 6000K blackbody. And 2 for nothing.

• test_mode (bool) – (Optional) Default = False to run as normal. This will overwrite the opacities and fix the delta tau at 0.5.

• full_output (bool) – (Optional) Default = False. If true, This will add taugas, taucld, tauray to the atmosphere class. This is done so that the users can debug, plot, etc.

• plot_opacity (bool) – (Optional) Default = False. If true, Will create a pop up plot of the weighted of each absorber at the middle layer

• return_mode (bool) – (Optional) Default = False, If true, will only return matrices for all the weighted opacity contributions

Returns:

• DTAU (ndarray) – This is a matrix with # layer by # wavelength. It is the opacity contained within a layer including the continuum, scattering, cloud (if specified), and molecular opacity If requested, this is corrected for with Delta-Eddington.

• TAU (ndarray) – This is a matrix with # level by # wavelength. It is the cumsum of opacity contained including the continuum, scattering, cloud (if specified), and molecular opacity If requested, this is corrected for with Delta-Eddington.

• WBAR (ndarray) – This is the single scattering albedo that includes rayleigh, raman and user input scattering sources. It has dimensions: # layer by # wavelength If requested, this is corrected for with Delta-Eddington.

• COSB (ndarray) – This is the asymettry factor which accounts for rayleigh and user specified values It has dimensions: # layer by # wavelength If requested, this is corrected for with Delta-Eddington.

• ftau_cld (ndarray) – This is the fraction of cloud opacity relative to the total TAUCLD/(TAUCLD + TAURAY)

• ftau_ray (ndarray) – This is the fraction of rayleigh opacity relative to the total TAURAY/(TAUCLD + TAURAY)

• GCOS2 (ndarray) – This is used for Cahoy+2010 methodology for accounting for rayleigh scattering. It replaces the use of the actual rayleigh phase function by just multiplying ftau_ray by 2

• DTAU (ndarray) – This is a matrix with # layer by # wavelength. It is the opacity contained within a layer including the continuum, scattering, cloud (if specified), and molecular opacity If requested, this is corrected for with Delta-Eddington.

• TAU (ndarray) – This is a matrix with # level by # wavelength. It is the cumsum of opacity contained including the continuum, scattering, cloud (if specified), and molecular opacity Original, never corrected for with Delta-Eddington.

• WBAR (ndarray) – This is the single scattering albedo that includes rayleigh, raman and user input scattering sources. It has dimensions: # layer by # wavelength Original, never corrected for with Delta-Eddington.

• COSB (ndarray) – This is the asymettry factor which accounts for rayleigh and user specified values It has dimensions: # layer by # wavelength Original, never corrected for with Delta-Eddington.

Notes

This was baselined against jupiter with the old fortran code. It matches 100% for all cases except for hotter cases where Na & K are present. This differences is not a product of the code but a product of the different opacities (1060 grid versus old 736 grid)

Todo

Add a better approximation than delta-scale (e.g. M.Marley suggests a paper by Cuzzi that has a better methodology)

picaso.optics.compute_raman(nwno, nlayer, wno, stellar_shifts, tlayer, cross_sections, j_initial, deltanu)

The Ramam scattering will alter the rayleigh scattering. The returned value is modified single scattering albedo. Cross sectiosn from: http://iopscience.iop.org/0004-637X/832/1/30/suppdata/apjaa3ec7t2_mrt.txt This method is described in Pollack+1986. Albeit not the best method. Sromovsky+2005 pointed out the inconsistencies in this method. You can see from his comparisons that the Pollack approximations don’t accurately capture the depths of the line centers. Since then, OKLOPCIC+2016 recomputed cross sections for J<=9. We are using those cross sections here with a different star. Huge improvement over what we had before. Number of J levels is hard coded to 10 ! Will be added to the rayleigh scattering as : TAURAY*RAMAN :param nwno: Number of wave points :type nwno: int :param nlayer: Number of layers :type nlayer: int :param wno: Array of output grid of wavenumbers :type wno: array :param stellar_shifts: Array of floats that has dimensions (n wave pts x n J levels: 0-9) :type stellar_shifts: ndarray :param tlayer: Array of floats that has dimensions nlayer :type tlayer: array :param cross_sections: The row of “C’s” from Antonija’s table. :type cross_sections: ndarray :param j_initial: The row of initial rotational energy states from Antonija’s table :type j_initial: ndarray :param deltanu: The row of delta nu’s from Antonia’s table :type deltanu: ndarray

picaso.optics.convert_array(clas, text)

needed to interpret bytes to array

picaso.optics.find_nearest(array, value)

picaso.optics.interp_matrix(a1, a2, a3, a4, t_interp, p_interp)

picaso.optics.j_fraction(j, T)

This computes the fraction of molecules at each J level. :param j: The initial rotational levels ranging from J=0 to 9 for hydrogen only :type j: int

Returns:

f_J in eqn. 3 https

Return type:

//arxiv.org/pdf/1605.07185.pdf

picaso.optics.numba_cumsum(mat)

Function to compute cumsum along axis=0 to bypass numba not allowing kwargs in cumsum

picaso.optics.partition_function(j, T)

Given j and T computes partition function for ro-vibrational levels of H2. This is the exponential and the statistical weight g_J in Eqn 3 in https://arxiv.org/pdf/1605.07185.pdf It is also used to compute the partition sum Z. :param j: Energy level :type j: int :param T: Temperature at each atmospheric layer in the atmosphere :type T: array float

Return type:

Returns partition function

picaso.optics.partition_sum(T)

This is the total partition sum. I am truncating it at 20 since it seems to approach 1 around then. This is also pretty fast to compute so 20 is fine for now. :param T: Array of temperatures for each layer :type T: array

Return type:

Z, the partition sum

picaso.optics.raman_pollack(nlayer, wave)

Mystery raman scattering. Couldn’t figure out where it came from.. so discontinuing. Currently function doesnt’ totally work. In half fortran-half python. Legacy from fortran albedo code. This method is described in Pollack+1986. Albeit not the best method. Sromovsky+2005 pointed out the inconsistencies in this method. You can see from his comparisons that the Pollack approximations don’t accurately capture the depths of the line centers. Since then, OKLOPCIC+2016 did a much better model of raman scattring (with ghost lines). Might be worth it to consider a more sophisticated version of Raman scattering. Will be added to the rayleigh scattering as : TAURAY*RAMAN

#OLD FORTRAN CODE #constants h = 6.6252e-27 c = 2.9978e10 bohrd = 5.2917e-9 hmass = 1.6734e-24 rmu = .5 * hmass #set wavelength shift of the ramam scatterer shift_v0 = 4161.0 facip = h * c / ( 1.e-4 * 27.2 * 1.602e-12 ) facray = 1.e16 * bohrd ** 3 * 128. * np.pi ** 5 * bohrd ** 3 / 9. facv = 2.08 / 2.38 * facray / bohrd ** 2 / ( 8. * np.pi * np.pi * rmu * c * shift_v0 ) * h #cross section of the unshifted rayleigh and the vibrationally shifted rayleigh gli = np.zeros(5) wli = np.zeros(5) gri = np.zeros(5) wri = np.zeros(5) gli[:] = [1.296, .247, .297, .157, .003] wli[:] = [.507, .628, .733, 1.175, 2.526] gri[:] = [.913, .239, .440, .344, .064] wri[:] = [.537, .639, .789, 1.304, 3.263] alp = np.zeros(7) arp = np.zeros(7) alp[:] = [6.84, 6.96, 7.33, 8.02, 9.18, 11.1, 14.5 ] arp[:] = [3.66, 3.71, 3.88, 4.19, 4.70, 5.52, 6.88 ] omega = facip / wavelength #first compute extinction cross section for unshifted component #e.g. rayleigh alpha_l=0 alpha_r=0 for i in range(5):

alpha_l += gli[i] / ( wli[i] ** 2 - omega ** 2 ) alpha_r += gri[i] / ( wri[i] ** 2 - omega ** 2 )

alpha2 = (( 2. * alpha_r + alpha_l ) / 3. ) ** 2 gamma2 = ( alpha_l - alpha_r ) ** 2 qray = facray * ( 3. * alpha2 + 2./3. * gamma2 ) / wavelength ** 4 #next, compute the extinction cross section for vibrationally #shifted component ip = np.zeros(2) ip[:] = [int(omega/0.05), 5.0] ip = np.min(ip) + 1 f = omega / 0.05 - float(ip-1) alpha_pl = ( 1. - f ) * alp[ip] + f * alp[ip+1] alpha_pr = ( 1. - f ) * arp[ip] + f * arp[ip+1] alpha_p2 = (( 2. * alpha_pr + alpha_pl ) / 3. ) ** 2 gamma_p2 = ( alpha_pl - alpha_pr ) ** 2 qv = facv / SHIFT( WAVEL, -SHIFTV0 ) ** 4 * ( 3. * alpha_p2 + 2./3. * gamma_p2 )

picaso.optics.rayleigh_old(colden, gasmixing, wave, xmu, amu)

DISCONTINUED Rayleigh function taken from old albedo code. Keeping this modular, as we may want to swap out different methods to calculate rayleigh opacity :param colden: Column Density in CGS units :type colden: array of float :param mixingratios: Mixing ratio as function of altitude for [H2 He CH4 ] :type mixingratios: array of float :param wave: wavelength (microns) of grid. This should be in DECENDING order to comply with

everything else being in wave numberes

Parameters:

• xmu (array of float) – mean molecular weight of atmosphere in amu

• amu (float) – amu constant in grams

picaso.optics.spline(x, y, n, yp0, ypn)

picaso.test module

picaso.test.dlugach_test(single_phase='OTHG', multi_phase='N=1', output_dir=None, rayleigh=True, phase=True, method='Toon', stream=2, opd=0.2, toon_coefficients='quadrature', delta_eddington=False)

Test the flux against against Dlugach & Yanovitskij https://www.sciencedirect.com/science/article/pii/0019103574901675?via%3Dihub

There are two tests: rayleigh and constant_tau. This will run through the entire table of values (XXI). Users can run both or specify one. Constant_tau picks single value for the asymmetry factors g = [0.5,0.75,0.8,0.85,0.9] and tests single scattering abled. Rayleigh only varies the ssa (with rayleigh phase fun)

This routine contains the functionality to reproduce table XXI and determine the % error in various realms

Parameters:

• single_phase (str) – Single phase function to test with the constant_tau test. Can either be: “cahoy”,”TTHG_ray”,”TTHG”, and “OTHG”

• output_dir (str) – Output directory for results of test. Default is to just return the dataframe (None).

• rayleigh (bool) – Default is True. Turns on and off the rayleigh phase function test

• phase (float) – Default=True. Turns on and off the constant g phase function test

• Retuns –

• ------ –

• XXI (DataFrame of % deviation from Dlugach & Yanovitskij Table) –

picaso.test.madhu_test(rayleigh=True, isotropic=True, asymmetric=True, single_phase='TTHG_ray', output_dir=None)

Test the flux against against Madhu and Burrows https://arxiv.org/pdf/1112.4476.pdf

There are three tests: rayleigh, isotropic, and lambert. Reproduces Figure 2 of Madhu and Burrows paper.

Parameters:

• rayleigh (bool) – Tests rayleigh phase function

• isotropic (bool) – Tests isotropic phase function

• single_phase (str) – Single phase function to test with the constant_tau test. Can either be: “cahoy”,”TTHG_ray”,”TTHG”, and “OTHG”

• output_dir (str) – Output directory for results of test. Default is to just return the dataframe (None).

Return type:

DataFrame of % deviation from Dlugach & Yanovitskij Table XXI

picaso.test.test_it_all()

picaso.test.thermal_sh_test(single_phase='OTHG', output_dir=None, phase=True, method='toon', stream=2, toon_coefficients='quadrature', delta_eddington=True, disort_data=False, phangle=0, tau=0.2)

Generate data to compare against DISORT. Must also run picaso_compare.py in pyDISORT :param single_phase: Single phase function to test with the constant_tau test. Can either be:

“cahoy”,”TTHG_ray”,”TTHG”, and “OTHG”

Parameters:

• output_dir (str) – Output directory for results of test. Default is to just return the dataframe (None).

• rayleigh (bool) – Default is True. Turns on and off the rayleigh phase function test

• phase (float) – Default=True. Turns on and off the constant g phase function test

• tau (float) – constant per layer opacity, 0.2 used as default in rooney+2023 tests

• Retuns –

• ------ –

• DataFrame –

picaso.wavelength module

picaso.wavelength.get_cld_input_grid(filename_or_grid='wave_EGP.dat', grid661=False)

The albedo code relies on the cloud code input, which is traditionally on a 196 wavelength grid. This method is to retrieve that grid. This file should be kept in the package reference data. Alternatively, it might be useful to add it directly to the input.cld file.

Parameters:

filename_or_grid (str or ndarray) – filename of the input grid OR a numpy array of wavenumber corresponding to cloud input

Returns:

array of wave numbers in increasing order

Return type:

array

picaso.wavelength.regrid(matrix, old_wno, new_wno)

This takes in a matrix that is (number of something versus number of wavelength points) and regrids to a new wave grid.

Parameters:

• matrix (ndarray) – matrix that is (number of something versus number of wavelength points)

• old_wno (array) – array that represents the old wavelength grid of the matrix

• new_wno (array) – array that represents the desired wavelength grid

Returns:

matrix that has been reinterpolated to new wave number grid

Return type:

matrix

Module contents