Given the refractive indices at a certain wavelength this module
calculates the Mie scattering by a stratified sphere.The basic code used
was that described in the report: ” Subroutines for computing the parameters of
the electromagnetic radiation scattered by a sphere ” J.V. Dave,
I B M Scientific Center, Palo Alto , California.
Report NO. 320 - 3236 .. MAY 1968 .
Parameters:
RO (float) – Outer Shell Radius (cm)
RFR (float) – Real refractive index of shell layer (in the form n= RFR-i*RFI)
RFI (float) – Imaginary refractive index of shell layer (in the form n= RFR-i*RFI)
THET (ndarray) – Angle in degrees between the directions of the incident and the scattered radiation.
JX (integer) – Total number of THET for which calculations are required
R (float) – Radius of core (cm)`
RE2 (float) – Real refractive index of core (in the form n= RE2-i*TMAG2)
TMAG2 (float) – Imaginary refractive index of core (in the form n= RE2-i*TMAG2)
WVNO (float) – Wave-number corresponding to the wavelength. (cm^-1)
Returns:
QEXT (float) – Efficiency factor for extinction,VAN DE HULST,P.14 ‘ 127
QSCAT (float) – Efficiency factor for scattering,VAN DE HULST,P.14 ‘ 127
CTBQRS (float) – Average(cos(theta))*QSCAT,VAN DE HULST,P.14 ‘ 127
ISTATUS (integer) – Convergence indicator, 0 if converged, -1 if otherwise.
Reads reference files with wavelength, and refractory indecies.
This function relies on input files being structured as a 4 column file with
columns: index, wavelength (micron), nn, kk
Parameters:
igas (str) – Gas name
directory (str) – Directory were reference files are located.
mw_atmos (float) – Mean molecular weight of the atmosphere amu
gas_mmr (float , optional) – Gas mass mixing ratio
None points to the default value of :
1xSolar = 2.55E-07
10xSolar = 2.1829E-06
50xSolar = 8.1164E-06
mh (float , optional) – Metallicity, Default is 1=1xSolar
Defines all the atmospheric parameters needed for the calculation
Note: Some of this is repeated in the layer() function.
This is done on purpose because layer is used to get a “virtual”
layer off the user input’s grid. These parameters are also
needed, though, to get the initial mixing parameters from the user.
Therefore, we put them in both places.
Computes kz profile and mixing length given user input. In brief the options are:
Input Kz
Input constant kz
3) Input convective heat flux (supply chf in df)
3a) Input convective heat flux, correct for latent heat (supply chf in df and set latent_heat=True)
and/or 3b) Input convective heat flux, correct for convective overshoot (supply chf, convective_overshoot=1/3)
4) Set kz_min to prevent kz from going too low (any of the above and set kz_min~1e5)
Parameters:
df (dataframe or dict) – Dataframe from input with “pressure”(bars),”temperature”(K). MUST have at least two
columns with names “pressure” and “temperature”.
Optional columns include the eddy diffusion “kz” in cm^2/s CGS units, and
the convective heat flux ‘chf’ also in cgs (e.g. sigma_csg T^4)
constant_kz (float) – Constant value for kz, if kz is supplied in df or filename,
it will inheret that value and not use this constant_value
Default = None
latent_heat (bool) – optional, Default = False. The latent heat factors into the mixing length.
When False, the mixing length goes as the scale height
When True, the mixing length is scaled by the latent heat
convective_overshoot (float) – Optional, Default is None. But the default value used in
Ackerman & Marley 2001 is 1./3. If you are unsure of what to pick, start
there. This is only used when the
This is ONLY used when a chf (convective heat flux) is supplied
kz_min (float) – Minimum Kz value. This will reset everything below kz_min to kz_min.
Default = 1e5 cm2/s
Define Kz in CGS. Should be on same grid as pressure. This overwrites whatever was
defined in get_pt ! Users can define kz by:
Defining a DataFrame with keys ‘pressure’ (in bars), and ‘kz’
Defining constant kz
Supplying a convective heat flux and prescription for latent_heat
Parameters:
df (pandas.DataFrame, dict) – Dataframe or dictionary with ‘kz’ as one of the fields.
constant_kz (float) – Constant value for kz in units of cm^2/s
chf (ndarray) – Convective heat flux in cgs units (e.g. sigma T^4). This will be used to compute
the kzz using the methodology of Gierasch and Conrath (1985)
latent_heat (bool) – optional, Default = False. The latent heat factors into the mixing length.
When False, the mixing length goes as the scale height
When True, the mixing length is scaled by the latent heat
This is ONLY used when a chf (convective heat flux) is supplied
df (dataframe or dict) – Dataframe with “pressure”(bars),”temperature”(K). MUST have at least two
columns with names “pressure” and “temperature”.
Optional columns include the eddy diffusion “kz” in cm^2/s CGS units, and
the convective heat flux ‘chf’ also in cgs (e.g. sigma_csg T^4)
filename (str) – Filename read in. Will be read in with pd.read_csv and should
result in two named headers “pressure”(bars),”temperature”(K).
Optional columns include the eddy diffusion “kz” in cm^2/s CGS units, and
the convective heat flux ‘chf’ also in cgs (e.g. sigma_csg T^4)
Use pd_kwargs to ensure file is read in properly.
kz_min (float, optional) – Minimum Kz value. This will reset everything below kz_min to kz_min.
Default = 1e5 cm2/s
constant_kz (float, optional) – Constant value for kz, if kz is supplied in df or filename,
it will inheret that value and not use this constant_value
Default = None
latent_heat (bool) – optional, Default = False. The latent heat factors into the mixing length.
When False, the mixing length goes as the scale height
When True, the mixing length is scaled by the latent heat
convective_overshoot (float) – Optional, Default is None. But the default value used in
Ackerman & Marley 2001 is 1./3. If you are unsure of what to pick, start
there. This is only used when the
This is ONLY used when a chf (convective heat flux) is supplied
Teff (float, optional) – Effective temperature. If None (default), Teff set to temperature at 1 bar
alpha_pressure (float) – Pressure at which we want fsed=alpha for variable fsed calculation
pd_kwargs (kwargs) – Pandas key words for file read in.
If reading old style eddysed files, you would need:
skiprows=3, sep=’s+’, header=None, names=[“ind”,”pressure”,”temperature”,”kz”]
Function that calculations new Mie database using PyMieScatt.
:param gas_name: List of names of gasses. Or a single gas name.
See pyeddy.available() to see which ones are currently available.
Parameters:
dir_refrind (str) – Directory where you store optical refractive index files that will be created.
dir_out (str) – Directory where you want to store Mie parameter files. Will be stored as gas_name.Mieff.
BEWARE FILE OVERWRITES.
rmin (float , optional) – (Default=1e-8) Units of cm. The minimum radius to compute Mie parameters for.
Usually 0.1 microns is small enough. However, if you notice your mean particle radius
is on the low end, you may compute your grid to even lower particle sizes.
nradii (int, optional) – (Default=60) number of radii points to compute grid on. 40 grid points for exoplanets/BDs
is generally sufficient.
Calculate spectrally-resolved profiles of optical depth, single-scattering
albedo, and asymmetry parameter for a user-input particle radius distribution
Parameters:
wave_in (ndarray) – your wavelength grid in microns
ndz (float) –
Column density of total particle concentration (#/cm^2)
Note: set to whatever, it’s your free knob
—- this does not directly translate to something physical because it’s for all particles in your slab
May have to use values of 1e8 or so
radius (ndarray) – Radius bin values - the range of particle sizes of interest. Maybe measured in the lab,
Ensure radius_unit is specified
radius_unit (astropy.unit.Units) – Astropy compatible unit
Top level program to run eddysed. Requires running Atmosphere class
before running this.
Parameters:
atmo (class) – Atmosphere class
directory (str, optional) – Directory string that describes where refrind files are
as_dict (bool, optional) – Default = False. Option to view full output as dictionary
og_solver (bool, optional) –
Default=True. BETA. Contact developers before changing to False.
Option to change mmr solver (True = original eddysed, False = new direct solver)
direct_tol (float , optional) – Only used if og_solver =False. Default = True.
Tolerance for direct solver
refine_TP (bool, optional) – Only used if og_solver =False.
Option to refine temperature-pressure profile for direct solver
og_vfall (bool, optional) – Option to use original A&M or new Khan-Richardson method for finding vfall
analytical_rg (bool, optional) – Only used if og_solver =False.
Option to use analytical expression for rg, or alternatively deduce rg from calculation
Calculation option will be most useful for future
inclusions of alternative particle size distributions
do_virtual (bool) – If the user adds an upper bound pressure that is too low. There are cases where a cloud wants to
form off the grid towards higher pressures. This enables that.
Returns:
opd, w0, g0 – Extinction per layer, single scattering abledo, asymmetry parameter,
All are ndarrays that are nlayer by nwave
dict – Dictionary output that contains full output. See tutorials for explanation of all output.
mh (float) – Atmospheric metallicity in NON log units (e.g. 1 for 1x solar)
sig (float) – Width of the log normal particle distribution
d_molecule (float) – diameter of atmospheric molecule (cm) (Rosner, 2000)
(3.711e-8 for air, 3.798e-8 for N2, 2.827e-8 for H2)
Set in Atmosphere constants
eps_k (float) – Depth of the Lennard-Jones potential well for the atmosphere
Used in the viscocity calculation (units are K) (Rosner, 2000)
c_p_factor (float) – specific heat of atmosphere (erg/K/g) . Usually 7/2 for ideal gas
diatomic molecules (e.g. H2, N2). Technically does slowly rise with
increasing temperature
og_vfall (bool , optional) – optional, default = True. True does the original fall velocity calculation.
False does the updated one which runs a tad slower but is more consistent.
The main effect of turning on False is particle sizes in the upper atmosphere
that are slightly bigger.
do_virtual (bool,optional) – optional, Default = True which adds a virtual layer if the
species condenses below the model domain.
Reads reference files with wavelength, and refractory indecies.
This function relies on input files being structured as a 4 column file with
columns: index, wavelength (micron), nn, kk
Parameters:
igas (str) – Gas name
directory (str) – Directory were reference files are located.
Calculate layer condensate properties by iterating on optical depth
in one model layer (convering on optical depth over sublayers)
gas_namestr
Name of condenstante
rho_pfloat
density of condensed vapor (g/cm^3)
t_layerfloat
Temperature of layer mid-pt (K)
p_layerfloat
Pressure of layer mid-pt (dyne/cm^2)
t_topfloat
Temperature at top of layer (K)
t_botfloat
Temperature at botton of layer (K)
p_topfloat
Pressure at top of layer (dyne/cm2)
p_botfloat
Pressure at botton of layer
kzfloat
eddy diffusion coefficient (cm^2/s)
mixlfloat
Mixing length (cm)
gravityfloat
Gravity of planet cgs
mw_atmosfloat
Molecular weight of the atmosphere
gas_mwfloat
Gas molecular weight
q_belowfloat
total mixing ratio (vapor+condensate) below layer (g/g)
supsatfloat
Super saturation factor
fsedfloat
Sedimentation efficiency coefficient (unitless)
bfloat
Denominator of exponential in sedimentation efficiency (if param is ‘exp’)
eps: float
Minimum value of fsed function (if param=exp)
z_topfloat
Altitude at top of layer
z_botfloat
Altitude at bottom of layer
z_alphafloat
Altitude at which fsed=alpha for variable fsed calculation
paramstr
fsed parameterisation
‘const’ (constant), ‘exp’ (exponential density derivation)
sigfloat
Width of the log normal particle distribution
mhfloat
Metallicity NON log soar (1=1xSolar)
rminfloat
Minium radius on grid (cm)
nradint
Number of radii on Mie grid
d_moleculefloat
diameter of atmospheric molecule (cm) (Rosner, 2000)
(3.711e-8 for air, 3.798e-8 for N2, 2.827e-8 for H2)
Set in Atmosphere constants
eps_kfloat
Depth of the Lennard-Jones potential well for the atmosphere
Used in the viscocity calculation (units are K) (Rosner, 2000)
c_p_factorfloat
specific heat of atmosphere (erg/K/g) . Usually 7/2 for ideal gas
diatomic molecules (e.g. H2, N2). Technically does slowly rise with
increasing temperature
og_vfallbool
Use original or new vfall calculation
Returns:
qc_layer (ndarray) – condenstate mixing ratio (g/g)
qt_layer (ndarray) – gas + condensate mixing ratio (g/g)
rg_layer (ndarray) – geometric mean radius of condensate cm
reff_layer (ndarray) – droplet effective radius (second moment of size distrib, cm)
ndz_layer (ndarray) – number column density of condensate (cm^-3)
q_below (ndarray) – total mixing ratio (vapor+condensate) below layer (g/g)
Sets up a PICASO-readable dataframe that inserts a wavelength dependent aerosol layer at the user’s
given pressure bounds, i.e., a wavelength-dependent slab of clouds or haze.
Parameters:
p_bottom (float) – the cloud/haze base pressure
the upper bound of pressure (i.e., lower altitude bound) to set the aerosol layer. (Bars)
opd (ndarray) – wavelength-dependent optical depth of the aerosol
w0 (ndarray) – wavelength-dependent single scattering albedo of the aerosol
wavenumber_grid (ndarray) – wavenumber grid in (cm^-1)
pressure_grid (ndarray) – bars, user-defined pressure grid for the model atmosphere
p_top (float) – bars, the cloud/haze-top pressure
This cuts off the upper cloud region as a step function.
You must specify either p_top or p_decay.
p_decay (ndarray) – noramlized to 1, unitless
array the same size as pressure_grid which specifies a
height dependent optical depth. The usual format of p_decay is
a fsed like exponential decay ~np.exp(-fsed*z/H)
Return type:
Dataframe of aerosol layer with pressure (in levels - non-physical units!), wavenumber, opd, w0, and g0 to be read by PICASO
Plot PT profiles with condensation curves.
Thick lines in this pt plot signify those gases that
were turned on in the run, NOT those that were recommended.
Parameters:
out (dict) – Dictionary output from pyeddy run
with_condensation (bool) – Plots condensation curves of gases. Also plots those that were turned on in the
calculation with line_width=5. All others are set to 1.
return_condensation (bool) – If true, it returns list of condenation temps for each gas, pressure grid,
and a list of the gas names
**kwargs (kwargs) – Kwargs for bokeh.figure()
Returns:
if return_condensation (fig, cond temperature curves, ,cond pressure grid , name of all available gases)
”
Define force balance for spherical particles falling in atmosphere, namely equate
gravitational and viscous drag forces.
Viscous drag assumed to be quadratic (Benchaita et al. 1983) and is a function of
the Reynolds number dependent drag coefficient.
Drag coefficient taken from Khan-Richardson model (Richardson et al. 2002) and
is valid for 1e-2 < Re < 1e5.
Calculates the pressure of saturation mixing ratio for a gas
that is extrapolated below the model domain.
This is specifically used if the gas has saturated
below the model grid
Calculate fallspeed for a particle at one layer in an
atmosphere, assuming low Reynolds number and in the molecular regime (Epstein drag),
i.e., where Knudsen number (mfp/r) is high (>>1).
User chooses whether particle is spherical or an aggregate.
If aggregrate, the monomer size is given by r and the outer effective radius of the aggregrate is Ragg.
Calculates fallspeed for a fractal aggegrate particle as performed
by Ohno et al., 2020, with an outer
effective radius of R_agg made up of monomers of size r, at one layer in an
atmosphere. ***Requires characteristic radius ragg to have D > or equal 2, i.e., not very fluffy aggregrates.
Essentially a more explicit version of the regular “vfall_aggregates” function,
and is not dependent on being in free molecular (Esptein) regime
all units are cgs
Parameters:
r (float) – monomer particle radius (cm)
grav (float) – acceleration of gravity (cm/s^2)
rhop (float) – density of aggregrate particle (g/cm^3)
This is used to find F(X)-y=0 where F(X) is the fall speed for
a spherical particle and y is the convective velocity scale.
When the two are balanced we arrive at our correct particle radius.
Therefore, it is the same as the vfall function but with the
subtraction of the convective velocity scale (cm/s) w_convect.