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
Have the following column names (any order) opd g0 w0
Be white space delimeted
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:
- Version:
2015.01.29
Requirements¶
References
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
- 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, calculation='reflected')¶
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, verbose=True)¶
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, reorder_output=False, **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
reorder_output (bool) – Returns sorted outputs, for better phase curve plotting. re-orders phases such that brightest value is at phase=0
- 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
- 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