protoRT.rtcube

Created on Thu May 26 11:01:11 2021

@author: daniel

Module Contents

Classes

RadiativeTransferCube

Class for radiative transfer through a 3D dust density cube and for computing the mass excess.

Functions

particles_to_density(xp, yp, zp, x, y, z[, ...])	Deposit particles onto a 3D grid to generate a dust density field.
find_index_bisect(qpar, q)	Find the index of the element in a sorted list closest to a target value using binary search.
load_cube()	Load a 256×256×256 dust density cube from a simulation snapshot.

class protoRT.rtcube.RadiativeTransferCube(data: Optional[numpy.ndarray] = None, axis: Optional[numpy.ndarray] = None, code_rho: float = 1.0, code_cs: float = 2 * np.pi, code_omega: float = 2 * np.pi, column_density: float = 10.0, T: float = 30.0, H: float = 5 * const.au.cgs.value, stoke: Union[float, numpy.ndarray] = 0.3, grain_rho: Union[float, numpy.ndarray] = 1.675, wavelength: float = 0.1, include_scattering: bool = False, kappa: Optional[Union[float, numpy.ndarray]] = None, sigma: Optional[Union[float, numpy.ndarray]] = None, p: float = 2.5, npar: int = int(1000000.0), ipars: Optional[numpy.ndarray] = None, xp: Optional[numpy.ndarray] = None, yp: Optional[numpy.ndarray] = None, zp: Optional[numpy.ndarray] = None, xgrid: Optional[numpy.ndarray] = None, ygrid: Optional[numpy.ndarray] = None, zgrid: Optional[numpy.ndarray] = None, rhopswarm: Optional[numpy.ndarray] = None, particle_weight: Optional[Union[float, numpy.ndarray]] = 6.30813659928, grid_func: Optional[str] = 'linear', num_grid_points: Optional[int] = 262, num_interp_points: Optional[int] = 256, index_limits_1: Optional[int] = 3, index_limits_2: Optional[int] = 258, aps: Optional[numpy.ndarray] = None, eps_dtog: Optional[float] = None, init_var: Optional[numpy.ndarray] = None)

Class for radiative transfer through a 3D dust density cube and for computing the mass excess.

This class computes the optical depth and outgoing intensity through a protoplanetary disk model by performing radiative transfer along the z-axis of a 3D density cube. It supports both monodisperse and polydisperse dust distributions and can automatically compute dust opacities using the DSHARP model (Birnstiel et al. 2018) if they are not provided.

It also computes the mass excess — defined as the ratio of the true dust mass in the cube to the mass inferred from the observed intensity under the assumption of optically thin emission. This metric quantifies how much disk mass may be underestimated in observations due to optically thick regions.

Note

• The current version of the code assumes that the cube is symmetric along all axes (i.e., cubic domain).

• If no input data is provided, a default test dataset from a streaming instability simulation by Yang & Johansen (2014) is used. This is a snapshot at orbit 100.

Parameters:

• data (ndarray, optional) – 3D dust density cube (e.g., rhop field). Defaults to test data if None.

• axis (ndarray, optional) – 1D coordinate array along the integration axis (must be z), in units of gas scale height.

• code_rho (float, optional) – Midplane gas density in code units (rho0 in start.in). Default is 1.

• code_cs (float, optional) – Sound speed in code units (cs0 in start.in). Default is 2π.

• code_omega (float, optional) – Orbital frequency in code units (Omega in start.in). Default is 2π.

• column_density (float, optional) – Gas column density in cgs units (g/cm²). Default is 10.

• T (float, optional) – Isothermal box temperature (K). Default is 30.

• H (float, optional) – Pressure scale height (cm). Default is 7.5e13 (≈ 5 AU).

• stoke (float or ndarray) – Stokes number(s) for the dust population.

• grain_rho (float or ndarray, optional) – Internal grain density (g/cm³). Must correspond with the input Stokes number(s). If None, defaults to 1.675, which is from the DSHARP dust model and must be used if no opacities are input as DSHARP opacities are used by default.

• wavelength (float, optional) – Wavelength at which opacities are computed (cm). Range: 1e-5 to 10. Default is 0.1.

• include_scattering (bool, optional) – If True, include scattering opacity in radiative transfer. Default is False.

• kappa (float or ndarray, optional) – Dust absorption opacity (cm²/g). Must correspond with the input Stokes number(s). If None, computed from DSHARP.

• sigma (float or ndarray, optional) – Dust scattering opacity (cm²/g). Must correspond with the input Stokes number(s). If None, computed from DSHARP.

• p (float, optional) – Power-law index for grain size distribution (n(a) ∝ a^{-p}). Default is 2.5.

• npar (int, optional) – Number of particles in the simulation. Default is 1,000,000. Only required for multi-species simulations.

• ipars (ndarray, optional) – Particle species and positions. Only required for multi-species simulations.

• xp (ndarray, optional) – Particle species and positions. Only required for multi-species simulations.

• yp (ndarray, optional) – Particle species and positions. Only required for multi-species simulations.

• zp (ndarray, optional) – Particle species and positions. Only required for multi-species simulations.

• xgrid (ndarray, optional) – Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations.

• ygrid (ndarray, optional) – Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations.

• zgrid (ndarray, optional) – Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations.

• rhopswarm (ndarray, optional) – Local swarm densities used in density map conversion and proto-mass calculation. Only required for polydisperse and/or self-gravitating simulations.

• particle_weight (float or array_like, optional) – Particle density or weight. If a float, all particles have the same weight, if array, must correspond with the particle positions. Pencil Code stores this as the rhop_swarm attribute in read_param(). Default is 6.30813659928. Only required for multi-species simulations.

• grid_func ({'linear'}, optional) – Grid interpolation scheme for each axis. Only ‘linear’ is currently supported. Pencil Code stores these as attributes in read_param(). Default is ‘linear’. Only required for multi-species simulations.

• num_grid_points (int) – Number of grid points in x, y, and z directions. Current code assumes cubic domain. Pencil code stores these as attributes in read_dim(). Default is 262. Only required for multi-species simulations.

• num_interp_points (int) – Number of interpolation points in x, y, and z. Current code assumes cubic domain. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim(). Default is 256. Only required for multi-species simulations.

• index_limits_1 (int) – Grid index limits along x (l1), y (m1), and z (n1) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim(). Default is 3. Only required for multi-species simulations.

• index_limits_2 (int) – Grid index limits along x (l2), y (m2), and z (n2) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim(). Default is 258. Only required for multi-species simulations.

• aps (ndarray, optional) – Azimuthal positions of particles for proto-mass calculation. Only required for self-gravitating simulations. Must be input to enable protomass calculations.

• eps_dtog (float, optional) – Dust-to-gas ratio used to estimate proto-masses. Only required for self-gravitating simulations. Must be input to enable protomass calculations.

• init_var (ndarray, optional) – Initial density cube for mass excess calculation. Only required for self-gravitating simulations. In these cases the initial mass should be considered in the mass excess calculation.

axis

The 1D coordinate array corresponding to the integration axis (must be the z-axis), in units of gas scale height.

Type:

ndarray

data

The 3D dust density cube (rhop field), used for radiative transfer.

Type:

ndarray

filling_factor

Fraction of columns with τ ≥ 1 (optically thick).

Type:

float

grain_size

Derived grain size(s) based on Stokes number(s).

Type:

float or ndarray

H

Gas pressure scale height of the box (cm).

Type:

float

intensity

Outgoing intensity map computed along the z-axis.

Type:

ndarray

mass

Total dust mass in the cube (g).

Type:

float

mass_excess

Ratio of true dust mass to inferred dust mass (from the optically thin approximation).

Type:

float

proto_mass

Total mass in gravitationally bound clumps (g), if simulation is self-gravitating.

Type:

float

tau

Optical depth map computed along the z-axis.

Type:

ndarray

configure()

Initializes internal attributes and computes the opacity (if not provided), calculates the optical depth, performs the radiative transfer, and then computes the mass excess.

blackbody()

Computes blackbody intensity from Planck’s law at temperature T.

calc_tau()

Computes the 2D optical depth (τ) map.

calc_t(rhod, kappa, sigma)

Computes cumulative τ along one vertical column.

calc_intensity()

Solves the radiative transfer equation along the z-axis.

calc_mass_excess()

Calculates the dust mass underestimation factor.

calc_grain_size()

Converts Stokes number(s) into physical grain size(s).

extract_opacity()

Computes opacities using the DSHARP dust model.

get_proto_mass()

Estimates mass in bound clumps (protoplanets) using swarm density and azimuthal info.

configure()

Initialize and configure physical parameters for the density cube object.

This method performs unit conversions, sets up spatial dimensions, computes the total mass of the cube in cgs units, and prepares quantities required for radiative transfer and mass excess analysis.

It must be re-run if key inputs, such as the scale height H, the absorption/scattering opacities (kappa, sigma), or grid spacing (axis), are updated.

Raises:

ValueError – If self.data or self.axis are not set prior to calling this method.

Returns:

This method updates the following internal attributes in place:

• self.Lx, self.Ly, self.Lzfloat

  Physical dimensions of the cube (cm).

• self.dx, self.dy, self.dzfloat

  Grid cell spacing (cm).

• self.massfloat

  Total dust mass in the cube (g).

• self.unit_massfloat

  Conversion factor from code mass units to CGS (g).

• self.unit_densityfloat

  Conversion factor from code density units to CGS (g/cm³).

• self.unit_sigmafloat

  Conversion factor for dust surface density (g/cm²).

• self.datandarray

  Density cube converted to CGS units (g/cm³).

• self.frequencyfloat

  Frequency corresponding to the specified wavelength (Hz).

• self.mass_excessfloat

  Ratio of true dust mass to inferred mass assuming optically thin emission.

• self.taundarray

  Optical depth map (computed in later steps).

• self.proto_massfloat

  Total protoplanet mass (if aps and rhopswarm are available).

Return type:

None

blackbody()

Compute the blackbody spectral radiance using Planck’s law.

This method calculates the specific intensity (spectral radiance) of a blackbody at temperature T and frequency self.frequency, assuming pure absorption (i.e., no scattering contribution). The result is stored in self.B_nu.

Raises:

ValueError – If self.frequency or self.T is not set, or if T <= 0.

Returns:

Updates the self.B_nu attribute with the spectral radiance in units of erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹.

Return type:

None

calc_tau()

Compute the optical depth map (self.tau) by integrating along the z-axis of the density cube.

For monodisperse simulations, a uniform opacity coefficient is applied throughout. For polydisperse simulations, species-specific densities are reconstructed in 3D, and local weighted opacity coefficients are calculated at each cell using the particle data.

This method updates self.tau for all (x, y) columns and computes the self.filling_factor, defined as the fraction of columns with optical depth τ ≥ 1.

Notes

• self.data must be in cgs units (g/cm³).

• self.dz must be the vertical grid spacing in cm.

• self.kappa (absorption opacity) is required.

• If include_scattering is True, self.sigma (scattering opacity) must also be provided.

• For polydisperse simulations:

  • self.grain_size must be an array.

  • self.ipars, self.xp, and self.xgrid must be defined to compute density per species.

Raises:

ValueError – If required attributes (self.data, self.dz, self.kappa, etc.) are not defined. If any necessary polydisperse inputs are missing.

Returns:

The following attributes are updated in place:

• self.taundarray

  2D optical depth map integrated along the z-axis for each (x, y) column.

• self.filling_factorfloat

  Fraction of columns with τ ≥ 1 (optically thick).

• self.effective_kappandarray

  3D array of locally averaged absorption opacities (polydisperse only).

• self.effective_sigmandarray

  3D array of locally averaged scattering opacities (if include_scattering=True).

• self.density_per_speciesndarray

  4D array of species-separated dust density grids (polydisperse only).

Return type:

None

calc_t(rhod, effective_kappa, effective_sigma)

Compute the cumulative optical depth along a single vertical column.

This method integrates the optical depth from the bottom of the cube (z = 0) upward to each height z, returning a 1D array of cumulative τ values. Unlike calc_tau, which computes total τ for an entire column, this returns the depth-dependent progression of optical depth as a function of height.

Parameters:

• rhod (ndarray) – 1D array of dust mass densities along the vertical column (g/cm³).

• effective_kappa (ndarray) – 1D array of absorption opacity coefficients along the column (cm²/g).

• effective_sigma (ndarray) – 1D array of scattering opacity coefficients along the column (cm²/g).

Notes

• All input arrays must be 1D and of length self.Nz.

• The grid spacing self.dz is assumed to be in cm and must be positive.

Raises:

ValueError – If input arrays are not 1D or do not match the vertical grid size self.Nz, or if self.dz is not positive.

Returns:

t – 1D array of cumulative optical depth values at each vertical cell (unitless).

Return type:

ndarray

calc_intensity()

Compute the outgoing intensity along the z-axis for all (x, y) columns.

This method solves the radiative transfer equation using the emissivity term only, assuming no incident intensity at the base (i.e., no back-illumination). It integrates the column-wise emissivity weighted by the exponential attenuation factor exp(−τ) to determine the observed intensity at the top of the cube.

Scattering is included if self.include_scattering is True, in which case the mean intensity self.J is also used in the emissivity term. Supports both monodisperse and polydisperse models.

Notes

• Requires self.tau to be precomputed using calc_tau().

• Assumes that self.data contains dust density in cgs units (g/cm³).

• Assumes Planck function (self.B_nu) is constant along each column.

• The returned intensity is in units of erg s⁻¹ cm⁻² Hz⁻¹ sr⁻¹.

• If scattering is enabled, self.effective_sigma and self.J must be defined.

Raises:

ValueError – If required attributes such as self.data, self.tau, or opacity arrays are missing or have incorrect shapes.

Returns:

The following attribute is updated:

• self.intensityndarray

  2D map of outgoing specific intensity at the top of the cube for each (x, y) position.

Return type:

None

calc_mass_excess()

Calculate the mass underestimation factor (“mass excess”) for the dust density cube.

This method solves the radiative transfer equation along each column, calculates the outgoing intensity, and derives the observed dust mass assuming optically thin emission. The ratio of the true dust mass to the inferred dust mass is stored as self.mass_excess, quantifying how much mass is missed due to optically thick regions or scattering.

Scattering is handled via a two-stream approximation (Miyake & Nakagawa 1993; Zhu et al. 2019), which modifies the source function. Supports both monodisperse and polydisperse models.

Notes

• Requires self.kappa to be defined. If include_scattering=True, self.sigma is also required.

• The cube must be configured first with configure(), which sets self.mass and grid info.

• The outgoing intensity is computed using calc_intensity(), based on the full RT solution.

• Assumes thermal equilibrium (LTE) and no external illumination (I₀ = 0).

• For polydisperse simulations, the average opacity assumed in the optically thin approximation is taken to be the mean of the absorption opacities across species.

Raises:

ValueError – If required attributes such as density data, opacities, or mass are not initialized.

Returns:

The following attributes are updated:

• self.src_fnfloat or ndarray

  The effective source function used in radiative transfer integration.

• self.intensityndarray

  2D array of emergent intensity for each (x, y) column.

• self.observed_massfloat

  The dust mass inferred under the assumption of optically thin emission.

• self.mass_excessfloat

  The true mass divided by the observed mass. A value > 1 indicates underestimation.

Return type:

None

calc_grain_size()

Compute the grain size(s) from the Stokes number and gas column density.

This method uses the Stokes–Epstein drag relation to calculate the physical grain size a corresponding to a given Stokes number St, assuming all grains are in the Epstein regime. For monodisperse models, a single value is returned; for polydisperse models, an array of grain sizes is computed based on the per-species Stokes number and internal grain density.

The grain size is given by:

a = (St * 2 * Σ_g) / (π * ρ_g)

where:

• a is the grain size (cm),

• St is the Stokes number (dimensionless),

• Σ_g is the gas column density (g/cm²),

• ρ_g is the internal grain density (g/cm³).

Notes

• Assumes Epstein drag regime is valid.

• For polydisperse simulations, both stoke and grain_rho must be arrays of equal length.

• Requires self.column_density > 0.

Raises:

ValueError – If stoke or grain_rho is not set, non-positive, or if their lengths mismatch in the polydisperse case.

Returns:

The computed grain sizes are stored in the following attribute:

• self.grain_sizefloat or ndarray

  Computed grain size(s) in cm.

Return type:

None

extract_opacity()

Compute and assign dust opacity coefficients using the DSHARP model.

This method calculates the absorption (kappa) and scattering (sigma) opacity coefficients based on the dust grain sizes in the simulation. Opacities are derived using the DSHARP dust opacity tables (Birnstiel et al. 2018), which assume a fixed internal grain composition and density.

For polydisperse simulations (grain_size is an array), the method computes opacity values using a binned approximation consistent with the multi-species model.

Notes

• The grain_size and wavelength attributes must be set prior to calling this method.

• The method internally calls the compute_opacities module (compute_opacities.dsharp_model()).

Raises:

ValueError – If grain_size or wavelength is not set or invalid (implicitly by the underlying function in the compute_opacities module).

Returns:

Updates the following attributes in place:

• self.kappafloat or ndarray

  Absorption opacity coefficient(s) (cm²/g).

• self.sigmafloat or ndarray

  Scattering opacity coefficient(s) (cm²/g).

• self.grain_size_binsndarray

  Grain size bin edges used in the binned opacity approximation.

Return type:

None

get_proto_mass()

Estimate the masses of protoplanets (planetesimal clumps) from local particle densities.

This method calculates the mass of planetesimal clumps based on the local swarm density (rhopswarm), azimuthal particle positions (aps), and the total box dust mass. The particle mass is inferred from the total number of particles and the global dust-to-gas ratio.

Clumps are identified via non-zero azimuthal particle positions, and individual clump masses are computed as:

M_p = 10 ** [log10(N_clump) + log10(M_particle)]

where: - N_clump is the number of particles in a given clump (estimated from density), - M_particle is the dust mass per particle.

Notes

• Requires configure() to have been run (sets self.mass).

• Assumes that rhopswarm is proportional to the number of particles per grid cell.

• Returns clump masses in cgs units (g), sorted in ascending order.

• If no clumps are found, self.proto_mass is set to 0.

Raises:

ValueError – If required attributes (eps_dtog, mass, rhopswarm, aps, npar) are not defined.

Returns:

The resulting clump masses are stored in the self.proto_mass attribute:

• self.proto_massfloat or ndarray

  Sorted array of clump masses in grams, or 0 if no clumps are detected.

Return type:

None

protoRT.rtcube.particles_to_density(xp, yp, zp, x, y, z, rhop_swarm=6.30813659928, grid_func1='linear', grid_func2='linear', grid_func3='linear', mx=262, my=262, mz=262, nx=256, ny=256, nz=256, n1=3, n2=258, m1=3, m2=258, l1=3, l2=258, density=True)

Deposit particles onto a 3D grid to generate a dust density field.

This function interpolates particle mass or count contributions onto a 3D mesh using a quadratic weighting scheme adapted from an IDL routine by Anders Johansen. The method accounts for non-uniform grid spacing and allows interpolation weights to be computed using either particle mass (rhop_swarm) or uniform weighting.

Parameters:

• xp (array_like) – Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files.

• yp (array_like) – Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files.

• zp (array_like) – Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files.

• x (array_like) – Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid().

• y (array_like) – Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid().

• z (array_like) – Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid().

• rhop_swarm (float or array_like, optional) – Particle density or weight. If a float, all particles have the same weight; if an array, must match length of xp. Default is 6.30813659928. Pencil Code stores this as attribute in read_param().

• grid_func1 ({'linear'}, optional) – Grid interpolation scheme for each axis. Only ‘linear’ is currently supported. Pencil Code stores these as attributes in read_param().

• grid_func2 ({'linear'}, optional) – Grid interpolation scheme for each axis. Only ‘linear’ is currently supported. Pencil Code stores these as attributes in read_param().

• grid_func3 ({'linear'}, optional) – Grid interpolation scheme for each axis. Only ‘linear’ is currently supported. Pencil Code stores these as attributes in read_param().

• mx (int) – Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim().

• my (int) – Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim().

• mz (int) – Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim().

• nx (int) – Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim().

• ny (int) – Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim().

• nz (int) – Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim().

• l1 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• l2 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• m1 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• m2 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• n1 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• n2 (int) – Grid index limits along x (l), y (m), and z (n) directions for trimming ghost zones. Pencil code stores these as attributes in read_dim().

• density (bool, optional) – If True, return a mass-weighted density field. If False, return particle count field. Default is True.

Returns:

3D array of shape (n2 - n1 + 1, m2 - m1 + 1, l2 - l1 + 1), representing the deposited density or particle count on the grid.

Return type:

ndarray

Notes

• The interpolation uses a second-order accurate kernel.

• The resulting array excludes ghost zones as defined by index limits.

• Ensure xp, yp, zp fall within bounds of x, y, z to avoid edge errors.

References

• 23. Lyra, private communication (adapted from IDL code by A. Johansen)

protoRT.rtcube.find_index_bisect(qpar, q)

Find the index of the element in a sorted list closest to a target value using binary search.

This function returns the index of the element in the sorted array q that is nearest to the target value qpar, using an efficient bisection method.

Parameters:

• qpar (float) – The target value for which the nearest index is sought.

• q (array_like) – A 1D sorted array of floats (in ascending order) to search within.

Returns:

Index of the element in q closest to qpar.

Return type:

int

Notes

This method assumes that q is sorted in ascending order. If the array is not sorted, the result is undefined.

References

• 23. Lyra, private communication (adapted from IDL code by A. Johansen)

protoRT.rtcube.load_cube()

Load a 256×256×256 dust density cube from a simulation snapshot.

This function loads a 3D density cube corresponding to a single snapshot (orbit 100) of a 100-orbit streaming instability simulation in a protoplanetary disk, provided by Yang & Johansen (2014). The full cube is stored across two .npy files and is concatenated along the z-axis.

Returns:

• data (ndarray) – A 3D NumPy array of shape (256, 256, 256) representing the dust density cube.

• axis (ndarray) – A 1D NumPy array representing the spatial coordinates corresponding to each axis (z, y, x) in code units.