# Including clouds¶

## Powerlaw cloud models¶

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

```
[1]:
```

```
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.

```
[2]:
```

```
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

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:

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:

```
[3]:
```

```
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):

```
[4]:
```

```
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):

```
[5]:
```

```
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):

```
[6]:
```

```
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):

```
[7]:
```

```
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

```
[8]:
```

```
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:

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:

```
[10]:
```

```
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:

```
[25]:
```

```
# 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.

```
[26]:
```

```
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.

```
[27]:
```

```
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

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:

```
[28]:
```

```
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:

```
[29]:
```

```
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:

```
[30]:
```

```
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:

```
[31]:
```

```
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:

```
[32]:
```

```
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

```
[33]:
```

```
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).

```
[34]:
```

```
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).