Including clouds

Powerlaw cloud models

Let’s initialize an atmosphere in the usual way, see “Getting Started”:

import numpy as np
from petitRADTRANS import Radtrans

atmosphere = Radtrans(line_species = ['H2O_HITEMP',
                                      'CO_all_iso_HITEMP',
                                      'CH4',
                                      'CO2',
                                      'Na_allard',
                                      'K_allard'],
                      rayleigh_species = ['H2', 'He'],
                      continuum_opacities = ['H2-H2', 'H2-He'],
                      wlen_bords_micron = [0.3, 15])

pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

/Users/molliere/Documents/Project_Docs/petitRADTRANS/petitRADTRANS/radtrans.py:99: FutureWarning: pRT_input_data_path was set by an environment variable. In a future update, the path to the petitRADTRANS input_data will be set within a .ini file that will be automatically generated into the user home directory (OS agnostic), inside a .petitradtrans directory
  warnings.warn(f"pRT_input_data_path was set by an environment variable. In a future update, the path to "

  Read line opacities of H2O_HITEMP...
 Done.
  Read line opacities of CO_all_iso_HITEMP...
 Done.
  Read line opacities of CH4...
 Done.
  Read line opacities of CO2...
 Done.
  Read line opacities of Na_allard...
 Done.
  Read line opacities of K_allard...
 Done.

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
Done.

Units in petitRADTRANS: remember that all units in petitRADTRANS are in cgs, except for pressure, which is in bars, and the mean molecular weight (MMW), which is in units of atomic mass units.

Next we set up the atmospheric parameters, with a parameter selection identical to the “Getting Started” case.

from petitRADTRANS import nat_cst as nc
from petitRADTRANS.physics import guillot_global

R_pl = 1.838*nc.r_jup_mean
gravity = 1e1**2.45
P0 = 0.01

kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperature = guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ)

mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperature)
mass_fractions['He'] = 0.24 * np.ones_like(temperature)
mass_fractions['H2O_HITEMP'] = 0.001 * np.ones_like(temperature)
mass_fractions['CO_all_iso_HITEMP'] = 0.01 * np.ones_like(temperature)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperature)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Na_allard'] = 0.00001 * np.ones_like(temperature)
mass_fractions['K_allard'] = 0.000001 * np.ones_like(temperature)

MMW = 2.33 * np.ones_like(temperature)

Abundances in petitRADTRANS: remember that abundances in pRT are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). You can convert between mass fractions and VMRs by using \begin{equation} X_i = \frac{\mu_i}{\mu}n_i, \end{equation} where \(X_i\) is the mass fraction of species \(i\), \(\mu_i\) the mass of a single molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molecular weight, and \(n_i\) is the VMR of species \(i\).

Now, let’s calculate cloudy spectra! The first available mode is the so-called power law cloud. In this case the following opacity is added to the scattering cross-section: \begin{equation} \kappa = \kappa_0\left(\frac{\lambda}{\lambda_0}\right)^\gamma \end{equation} Where \(\kappa_0\) is the opacity in units of cm\(^2\)/g, at \(\lambda_0=0.35 \ {\rm \mu m}\). The power law index \(\gamma\) fixes the wavelength dependence, Rayleigh-like scattering would be obtained for \(\gamma=-4\). Hence \(\kappa_0\) and \(\gamma\) are free parameters.

We now calculate a clear atmosphere, and then four power law cloudy models, with \(\kappa_0 = 0.01\) cm\(^2\)/g, and four different \(\gamma\) values.

Clear model:

import pylab as plt
plt.rcParams['figure.figsize'] = (10, 6)

# Clear
atmosphere.calc_transm(temperature, mass_fractions, \
                       gravity, MMW, R_pl=R_pl, P0_bar=P0)
clear = atmosphere.transm_rad/nc.r_jup_mean

First, \(\gamma = -4\) (Rayleigh-like):

kappa_zero = 0.01
gamma_scat = -4.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m4 = atmosphere.transm_rad/nc.r_jup_mean

Second, \(\gamma = -2\) (weaker scattering power law):

kappa_zero = 0.01
gamma_scat = -2.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m2 = atmosphere.transm_rad/nc.r_jup_mean

Third, \(\gamma = 0\) (flat opacity):

kappa_zero = 0.01
gamma_scat = 0.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m0 = atmosphere.transm_rad/nc.r_jup_mean

Fourth, the exoctic case of, \(\gamma = 1\) (positive opacity slope):

kappa_zero = 0.01
gamma_scat = 1.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

p1 = atmosphere.transm_rad/nc.r_jup_mean

# Make plot

plt.plot(nc.c/atmosphere.freq/1e-4, clear, label = 'Clear')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m4, \
         label = r'Powerlaw cloud, $\gamma = -4$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m2, \
         label = r'Powerlaw cloud, $\gamma = -2$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m0, \
         label = r'Powerlaw cloud, $\gamma = 0$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         p1,\
         label = r'Powerlaw cloud, $\gamma = 1$')
plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

Gray cloud deck and/or scaled Rayhleigh scattering

In addition to power law clouds, petitRADTRANS can also simply add a gray cloud deck, or scale the Rayleigh scattering opacities of the gas by a factor specified by the user. For these examples we will reuse the Radtrans object from the demonstration of the power law cloud above, there’s no need to intialize it again.

Let’s calculate three spectra: - A clear model - A cloudy model with a gray cloud deck at 0.01 bar - A haze-like model, mimicked by scaling the Rayleigh scattering of the gas by a factor 10 - A combination of the two above cloud models: gray cloud deck plus haze

import pylab as plt
plt.rcParams['figure.figsize'] = (10, 6)

# Clear
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Clear')

# Gray cloud deck at 0.01 bar
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      Pcloud = 0.01)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Gray cloud deck at 0.01 bar')

# Haze (10 x gas Rayleigh scattering)
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      haze_factor = 10)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Rayleigh haze')

# Haze + cloud deck
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      haze_factor = 10, Pcloud = 0.01)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Rayleigh haze + cloud deck')

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

Defining arbitrary opacity functions (can be used to parameterize clouds)

What happens if you have a favorite cloud parameterization, but its implementation is not offered in pRT? Fear not, there is an option for you! When running calc_transm() (or calc_flux()) you can hand it any arbitrary function that returns a scattering or absorption opacity as a function of wavelength and pressure to pRT. This is then going to be used in that particular calc_transm() or calc_flux() run. An example can be found below, where we implement the following cloud model: \begin{equation} \kappa_{\rm cloud}(\lambda, P)= \frac{\kappa_0}{1+\left(\lambda/\lambda_0\right)^p} \left(\frac{P}{P_{\rm base}}\right)^{f_{\rm sed}} {\rm \ if \ } P and \(0\) if \(P\geq P_{\rm base}\). Here, \(\lambda\) and \(P\) are the wavelength in micron and pressure in bar, respectively, and \(\kappa_0\), \(\lambda_0\), \(p\), \(P_{\rm base}\) and \(f_{\rm sed}\) are free parameters, thought to parameterize the opacity (in \(\rm cm^2/g\)) at the reference wavelength \(\lambda_0\), the power law dependence of the opacity with wavelengh at \(\lambda \gg \lambda_0\), the cloud base pressure in bar, and the cloud scale height descrease factor, respectively.

Here an example implementation of this cloud model:

def cloud_opas(kappa_0,
               lambda_0,
               p,
               P_base,
               f_sed):
                     # in micron, in bar
    def give_opacity(wlen, press):

        retVal = np.zeros((len(wlen), len(press)))

        for i_p in range(len(press)):
            if press[i_p] < P_base:
                retVal[:, i_p] = kappa_0 / (1. + (wlen/lambda_0)**p) * (press[i_p]/P_base)**f_sed
            else:
                retVal[:, i_p] = 0.
        return retVal

    return give_opacity

We can now hand this function to calc_transm() or (calc_flux()). If we want to use it as an absorption opacity, we have to use the give_absorption_opacity keyword. If we want to use it as a scattering opacity, we have to use the give_scattering_opacity keyword. Note that in the transmission spectrum case the results will not differ between the two options, because both scattering and absorption cross-sections simply contribute to the total extinction opacity.

Below some examples:

# Clear
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0)
plt.plot(nc.c/atmosphere.freq/1e-4,
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Clear')

# External cloud absorption function, test 1
atmosphere.calc_transm(temperature,
                       mass_fractions,
                       gravity,
                       MMW,
                       R_pl=R_pl,
                       P0_bar=P0,
                       Pcloud = 0.01,
                       give_absorption_opacity=cloud_opas(1e-2, 1., 1., 0.01, 2))
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = r'User specified cloud function ' \
                   r'$(\kappa_0=0.01 \ {\rm cm^2/g}$, $\lambda_0=1 \ {\rm \mu m}$, $p=1$, ' \
                   r'$P_{\rm base} = 0.01  \ {\rm bar}$ and $f_{\rm sed}=2$)')

# External cloud absorption function, test 2
atmosphere.calc_transm(temperature,
                       mass_fractions,
                       gravity,
                       MMW,
                       R_pl=R_pl,
                       P0_bar=P0,
                       Pcloud = 0.01,
                       give_absorption_opacity=cloud_opas(1e1, 0.5, 4, 0.001, 1))
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = r'User specified cloud function ' \
                   r'$(\kappa_0=10 \ {\rm cm^2/g}$, $\lambda_0=0.5 \ {\rm \mu m}$, $p=4$, ' \
                   r'$P_{\rm base} = 0.001  \ {\rm bar}$ and $f_{\rm sed}=1$)')


plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

If you use this option in correlated-k (``c-k``) mode, please only parameterize opacities that vary slowly enough with wavelength. Molecular and atomic line opacities cannot be added in this way, for example, for these the correlated-k opacity treatment is required. Clouds (even if they have an absorption feature) and other pseudo-continuum opacities such as collision induced absorption (CIA) can be handled well with the treatment above, however.

Condensate clouds from real optical constants

Let’s calculate some spectra using opacities derived from optical constants of various materials. In this example we will start of with forsterite, that is \(\rm Mg_2SiO_4\). For the list of available cloud species see “Available opacity species”.

Choose your favorite particle setup: petitRADTRANS offers multiple opacity versions for every given condensate, assuming either spherical (Mie scattering) or irreagularly shaped (DHS method = distribution of hollow spheres) particles. Moreover, opacities assuming a crystalline and/or amorphous internal structure can be used, where available.

We set up the atmosphere like before, this time loading the opacity of crystalline (c), irregularly shaped particles (DHS method, abbreviated d). The mode identifier for the \(\rm Mg_2SiO_4\) opacity therefore is Mg2SiO4(c)_cd. (c) stands for condensed. These are all possible modes, note that for some species only crystalline or amorphous cross-sections are available.

Opacity code	Description
Species(c)_cm	crystalline Mie opacity of Species
Species(c)_cd	crystalline DHS opacity of Species
Species(c)_am	amorphous Mie opacity of Species
Species(c)_ad	amorphous DHS opacity of Species

But now, let’s start! We simply give one additional cloud_species list to the Radtrans object, containing the name of the opacity species we want to use.

import numpy as np
from petitRADTRANS import Radtrans

atmosphere = Radtrans(line_species = ['H2O_HITEMP',
                                      'CO_all_iso_HITEMP',
                                      'CH4',
                                      'CO2',
                                      'Na_allard',
                                      'K_allard'], \
      cloud_species = ['Mg2SiO4(c)_cd'], \
      rayleigh_species = ['H2', 'He'], \
      continuum_opacities = ['H2-H2', 'H2-He'], \
      wlen_bords_micron = [0.3, 15])

pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

/Users/molliere/Documents/Project_Docs/petitRADTRANS/petitRADTRANS/radtrans.py:99: FutureWarning: pRT_input_data_path was set by an environment variable. In a future update, the path to the petitRADTRANS input_data will be set within a .ini file that will be automatically generated into the user home directory (OS agnostic), inside a .petitradtrans directory
  warnings.warn(f"pRT_input_data_path was set by an environment variable. In a future update, the path to "

  Read line opacities of H2O_HITEMP...
 Done.
  Read line opacities of CO_all_iso_HITEMP...
 Done.
  Read line opacities of CH4...
 Done.
  Read line opacities of CO2...
 Done.
  Read line opacities of Na_allard...
 Done.
  Read line opacities of K_allard...
 Done.

  Read in opacity of cloud species Mg2SiO4 ...
 Done.

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
Done.

Units in petitRADTRANS: remember that all units in petitRADTRANS are in cgs, except for pressure, which is in bars, and the mean molecular weight (MMW), which is in units of atomic mass units.

Next we set up the atmospheric parameters, with a parameter selection almost identical to the “Getting Started” case. We now also add the mass fraction of the cloud species we are interested in. Note that the mass fractions only contains the name Mg2SiO4(c) now, the _cd suffix has to be omitted.

from petitRADTRANS import nat_cst as nc
R_pl = 1.838*nc.r_jup_mean
gravity = 1e1**2.45
P0 = 0.01

kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperature = guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ)

mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperature)
mass_fractions['He'] = 0.24 * np.ones_like(temperature)
mass_fractions['H2O_HITEMP'] = 0.001 * np.ones_like(temperature)
mass_fractions['CO_all_iso_HITEMP'] = 0.01 * np.ones_like(temperature)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperature)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Na_allard'] = 0.00001 * np.ones_like(temperature)
mass_fractions['K_allard'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

MMW = 2.33 * np.ones_like(temperature)

Abundances in petitRADTRANS: remember that abundances in petitCODE are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). You can convert between mass fractions and VMRs by using \begin{equation} X_i = \frac{\mu_i}{\mu}n_i, \end{equation} where \(X_i\) is the mass fraction of species \(i\), \(\mu_i\) the mass of a single molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molecular weight, and \(n_i\) is the VMR of species \(i\).

Setting the cloud particle size

We have to define a few additional parameters for the clouds, to start the calculation. We can give the mean particle size and width of the log-normal particle distribution:

radius = {}
radius['Mg2SiO4(c)'] = 0.00005*np.ones_like(temperature) # I.e. a 0.5-micron particle size (0.00005 cm)

sigma_lnorm = 1.05

Now, let’s calculate a clear and cloudy spectrum to compare:

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       radius = radius, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'cloudy', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       radius = radius, sigma_lnorm = sigma_lnorm)
plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'clear', zorder = 1)

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc='best')
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

Here one sees that there is a lot of additional absorption in the optical and near-IR. Also note the silicate (Si-O stretching mode) becoming visible at 10 micron.

Calculating the particle size

Alternatively, one can specify an eddy diffusion parameter (\(K_{zz}\), with units of cm\(^2\)/s), which expresses how strongly atmospheric mixing is. The cloud particle size is then controlled by specifying \(K_{zz}\), and the unitless settling parameter \(f_{\rm sed}\), which expresses the particles’ mass averaged settling velocity, when compared to the local atmospheric mixing speed. Also, a width for the log-normal particle size distribution must be specified:

Kzz = np.ones_like(temperature)*1e1**7.5
fsed = 2.
sigma_lnorm = 1.05

Again, let’s calculate a clear and cloudy spectrum to compare:

mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'cloudy', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'clear', zorder = 1)

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc='best')
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

The resulting mean particle sizes can be accessed like this:

plt.yscale('log')
plt.xscale('log')

plt.ylim([1e2,1e-6])

plt.ylabel('P (bar)')
plt.xlabel('Average particle size (microns)')

plt.plot(atmosphere.r_g[:,atmosphere.cloud_species.index('Mg2SiO4(c)')]/1e-4, pressures)
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

Hansen Distribution Clouds

While the Ackermann-Marley model used here typically uses a log-normal distribution for the cloud particle radii, it is also possible to use the Hansen+ (1971) distribution. To do this, you’ll need to set the dist parameter in calc_flux to hansen, and set the effective width of the distribution with the b_hans parameter. Alternately, you can specify both the distribution width and center using the b_hans and a_hans arguments. Note that the value of b_hans will be different from that of sigma_lnorm to give a similar distribution, as the effective width is weighted by the particle area. See the Hansen+ (1971) paper for more details.

Here, we’ll calculate the mean particle size as a function of altitude

mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)
Kzz = np.ones_like(temperature)*1e1**7.5
fsed = 2.
b_hans = 0.01
atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW,
                       Kzz = Kzz, fsed=fsed, b_hans=b_hans,dist='hansen')

plt.yscale('log')
plt.xscale('log')

plt.ylim([1e2,1e-6])

plt.ylabel('P (bar)')
plt.xlabel('Average particle size (microns)')

plt.plot(atmosphere.r_g[:,atmosphere.cloud_species.index('Mg2SiO4(c)')]/1e-4, pressures)
plt.show()
plt.clf()

<Figure size 1000x600 with 0 Axes>

Minimum particle size: the opacity of particles smaller than 1 nm is set to zero.

Emission spectra

We can use the Radtrans object created above to calculate emission spectra as well. Here we will also generate an instance of pRT with scattering turned on with do_scat_emis = True, (see “Scattering for Emission Spectra” for an example on how to do this in detail).

mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, \
         color = 'black', label = 'cloudy, no scattering', zorder = 1)

# Load scattering version of pRT
atmosphere = Radtrans(line_species = ['H2O_HITEMP',
                                      'CO_all_iso_HITEMP',
                                      'CH4',
                                      'CO2',
                                      'Na_allard',
                                      'K_allard'], \
      cloud_species = ['Mg2SiO4(c)_cd'], \
      rayleigh_species = ['H2', 'He'], \
      continuum_opacities = ['H2-H2', 'H2-He'], \
      wlen_bords_micron = [0.3, 15], \
      do_scat_emis = True)

pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm, \
                       add_cloud_scat_as_abs = True)
plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, \
         label = 'cloudy, including scattering', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, '-', \
         color = 'red', label = 'clear', zorder = 0)

plt.legend(loc='best')
plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Planet flux $F_\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$)')
plt.show()
plt.clf()

Emission scattering is enabled: enforcing test_ck_shuffle_comp = True
  Read line opacities of H2O_HITEMP...
 Done.
  Read line opacities of CO_all_iso_HITEMP...
 Done.
  Read line opacities of CH4...
 Done.
  Read line opacities of CO2...
 Done.
  Read line opacities of Na_allard...
 Done.
  Read line opacities of K_allard...
 Done.

  Read in opacity of cloud species Mg2SiO4 ...
 Done.

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
Done.

<Figure size 1000x600 with 0 Axes>

Here we plotted the clear spectrum, neglecting the cloud opacity, the cloudy spectrum only considering the absorption of the cloud particles, and the cloudy spectrum with scattering, using the scattering version of pRT.

Scattering and petitRADTRANS: remember that scattering is included for emission spectra in petitRADTRANS only if requested specifically when generating the Radtrans object, as it increases the runtime (see “Scattering for Emission Spectra” for an example how to do this).