:py:mod:`protoRT.rtcube` ======================== .. py:module:: protoRT.rtcube .. autoapi-nested-parse:: Created on Thu May 26 11:01:11 2021 @author: daniel Module Contents --------------- Classes ~~~~~~~ .. autoapisummary:: protoRT.rtcube.RadiativeTransferCube Functions ~~~~~~~~~ .. autoapisummary:: protoRT.rtcube.particles_to_density protoRT.rtcube.find_index_bisect protoRT.rtcube.load_cube .. py:class:: 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. :param data: 3D dust density cube (e.g., rhop field). Defaults to test data if None. :type data: ndarray, optional :param axis: 1D coordinate array along the integration axis (must be `z`), in units of gas scale height. :type axis: ndarray, optional :param code_rho: Midplane gas density in code units (`rho0` in `start.in`). Default is 1. :type code_rho: float, optional :param code_cs: Sound speed in code units (`cs0` in `start.in`). Default is 2π. :type code_cs: float, optional :param code_omega: Orbital frequency in code units (`Omega` in `start.in`). Default is 2π. :type code_omega: float, optional :param column_density: Gas column density in cgs units (g/cm²). Default is 10. :type column_density: float, optional :param T: Isothermal box temperature (K). Default is 30. :type T: float, optional :param H: Pressure scale height (cm). Default is 7.5e13 (≈ 5 AU). :type H: float, optional :param stoke: Stokes number(s) for the dust population. :type stoke: float or ndarray :param grain_rho: 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. :type grain_rho: float or ndarray, optional :param wavelength: Wavelength at which opacities are computed (cm). Range: 1e-5 to 10. Default is 0.1. :type wavelength: float, optional :param include_scattering: If True, include scattering opacity in radiative transfer. Default is False. :type include_scattering: bool, optional :param kappa: Dust absorption opacity (cm²/g). Must correspond with the input Stokes number(s). If None, computed from DSHARP. :type kappa: float or ndarray, optional :param sigma: Dust scattering opacity (cm²/g). Must correspond with the input Stokes number(s). If None, computed from DSHARP. :type sigma: float or ndarray, optional :param p: Power-law index for grain size distribution (n(a) ∝ a^{-p}). Default is 2.5. :type p: float, optional :param npar: Number of particles in the simulation. Default is 1,000,000. Only required for multi-species simulations. :type npar: int, optional :param ipars: Particle species and positions. Only required for multi-species simulations. :type ipars: ndarray, optional :param xp: Particle species and positions. Only required for multi-species simulations. :type xp: ndarray, optional :param yp: Particle species and positions. Only required for multi-species simulations. :type yp: ndarray, optional :param zp: Particle species and positions. Only required for multi-species simulations. :type zp: ndarray, optional :param xgrid: Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations. :type xgrid: ndarray, optional :param ygrid: Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations. :type ygrid: ndarray, optional :param zgrid: Grid coordinates used to bin particles into grid cells for opacity weighting. Only required for multi-species simulations. :type zgrid: ndarray, optional :param rhopswarm: Local swarm densities used in density map conversion and proto-mass calculation. Only required for polydisperse and/or self-gravitating simulations. :type rhopswarm: ndarray, optional :param particle_weight: 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. :type particle_weight: float or array_like, optional :param grid_func: 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. :type grid_func: {'linear'}, optional :param num_grid_points: 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. :type num_grid_points: int :param num_interp_points: 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. :type num_interp_points: int :param index_limits_1: 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. :type index_limits_1: int :param index_limits_2: 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. :type index_limits_2: int :param aps: Azimuthal positions of particles for proto-mass calculation. Only required for self-gravitating simulations. Must be input to enable protomass calculations. :type aps: ndarray, optional :param eps_dtog: Dust-to-gas ratio used to estimate proto-masses. Only required for self-gravitating simulations. Must be input to enable protomass calculations. :type eps_dtog: float, optional :param init_var: 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. :type init_var: ndarray, optional .. attribute:: axis The 1D coordinate array corresponding to the integration axis (must be the z-axis), in units of gas scale height. :type: ndarray .. attribute:: data The 3D dust density cube (rhop field), used for radiative transfer. :type: ndarray .. attribute:: filling_factor Fraction of columns with τ ≥ 1 (optically thick). :type: float .. attribute:: grain_size Derived grain size(s) based on Stokes number(s). :type: float or ndarray .. attribute:: H Gas pressure scale height of the box (cm). :type: float .. attribute:: intensity Outgoing intensity map computed along the z-axis. :type: ndarray .. attribute:: mass Total dust mass in the cube (g). :type: float .. attribute:: mass_excess Ratio of true dust mass to inferred dust mass (from the optically thin approximation). :type: float .. attribute:: proto_mass Total mass in gravitationally bound clumps (g), if simulation is self-gravitating. :type: float .. attribute:: tau Optical depth map computed along the z-axis. :type: ndarray .. method:: 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. .. method:: blackbody() Computes blackbody intensity from Planck's law at temperature T. .. method:: calc_tau() Computes the 2D optical depth (τ) map. .. method:: calc_t(rhod, kappa, sigma) Computes cumulative τ along one vertical column. .. method:: calc_intensity() Solves the radiative transfer equation along the z-axis. .. method:: calc_mass_excess() Calculates the dust mass underestimation factor. .. method:: calc_grain_size() Converts Stokes number(s) into physical grain size(s). .. method:: extract_opacity() Computes opacities using the DSHARP dust model. .. method:: get_proto_mass() Estimates mass in bound clumps (protoplanets) using swarm density and azimuthal info. .. py:method:: 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.Lz` : float Physical dimensions of the cube (cm). - `self.dx`, `self.dy`, `self.dz` : float Grid cell spacing (cm). - `self.mass` : float Total dust mass in the cube (g). - `self.unit_mass` : float Conversion factor from code mass units to CGS (g). - `self.unit_density` : float Conversion factor from code density units to CGS (g/cm³). - `self.unit_sigma` : float Conversion factor for dust surface density (g/cm²). - `self.data` : ndarray Density cube converted to CGS units (g/cm³). - `self.frequency` : float Frequency corresponding to the specified wavelength (Hz). - `self.mass_excess` : float Ratio of true dust mass to inferred mass assuming optically thin emission. - `self.tau` : ndarray Optical depth map (computed in later steps). - `self.proto_mass` : float Total protoplanet mass (if `aps` and `rhopswarm` are available). :rtype: None .. py:method:: 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⁻¹. :rtype: None .. py:method:: 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. .. rubric:: 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.tau` : ndarray 2D optical depth map integrated along the z-axis for each (x, y) column. - `self.filling_factor` : float Fraction of columns with τ ≥ 1 (optically thick). - `self.effective_kappa` : ndarray 3D array of locally averaged absorption opacities (polydisperse only). - `self.effective_sigma` : ndarray 3D array of locally averaged scattering opacities (if `include_scattering=True`). - `self.density_per_species` : ndarray 4D array of species-separated dust density grids (polydisperse only). :rtype: None .. py:method:: 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. :param rhod: 1D array of dust mass densities along the vertical column (g/cm³). :type rhod: ndarray :param effective_kappa: 1D array of absorption opacity coefficients along the column (cm²/g). :type effective_kappa: ndarray :param effective_sigma: 1D array of scattering opacity coefficients along the column (cm²/g). :type effective_sigma: ndarray .. rubric:: 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). :rtype: ndarray .. py:method:: 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. .. rubric:: 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.intensity` : ndarray 2D map of outgoing specific intensity at the top of the cube for each (x, y) position. :rtype: None .. py:method:: 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. .. rubric:: 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_fn` : float or ndarray The effective source function used in radiative transfer integration. - `self.intensity` : ndarray 2D array of emergent intensity for each (x, y) column. - `self.observed_mass` : float The dust mass inferred under the assumption of optically thin emission. - `self.mass_excess` : float The true mass divided by the observed mass. A value > 1 indicates underestimation. :rtype: None .. py:method:: 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³). .. rubric:: 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_size` : float or ndarray Computed grain size(s) in cm. :rtype: None .. py:method:: 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. .. rubric:: 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.kappa` : float or ndarray Absorption opacity coefficient(s) (cm²/g). - `self.sigma` : float or ndarray Scattering opacity coefficient(s) (cm²/g). - `self.grain_size_bins` : ndarray Grain size bin edges used in the binned opacity approximation. :rtype: None .. py:method:: 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. .. rubric:: 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_mass` : float or ndarray Sorted array of clump masses in grams, or 0 if no clumps are detected. :rtype: None .. py:function:: 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. :param xp: Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files. :type xp: array_like :param yp: Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files. :type yp: array_like :param zp: Particle positions in x, y, and z. Pencil Code stores these as attributes in the pvar files. :type zp: array_like :param x: Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid(). :type x: array_like :param y: Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid(). :type y: array_like :param z: Grid node coordinates along x, y, and z axes, respectively. Pencil Code stores these as attributes in read_grid(). :type z: array_like :param rhop_swarm: 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(). :type rhop_swarm: float or array_like, optional :param grid_func1: Grid interpolation scheme for each axis. Only `'linear'` is currently supported. Pencil Code stores these as attributes in read_param(). :type grid_func1: {'linear'}, optional :param grid_func2: Grid interpolation scheme for each axis. Only `'linear'` is currently supported. Pencil Code stores these as attributes in read_param(). :type grid_func2: {'linear'}, optional :param grid_func3: Grid interpolation scheme for each axis. Only `'linear'` is currently supported. Pencil Code stores these as attributes in read_param(). :type grid_func3: {'linear'}, optional :param mx: Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim(). :type mx: int :param my: Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim(). :type my: int :param mz: Number of grid points in x, y, and z directions. Pencil code stores these as attributes in read_dim(). :type mz: int :param nx: Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim(). :type nx: int :param ny: Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim(). :type ny: int :param nz: Number of interpolation points in x, y, and z. Controls smoothing kernel width. Pencil code stores these as attributes in read_dim(). :type nz: int :param l1: 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(). :type l1: int :param l2: 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(). :type l2: int :param m1: 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(). :type m1: int :param m2: 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(). :type m2: int :param n1: 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(). :type n1: int :param n2: 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(). :type n2: int :param density: If True, return a mass-weighted density field. If False, return particle count field. Default is True. :type density: bool, optional :returns: 3D array of shape (n2 - n1 + 1, m2 - m1 + 1, l2 - l1 + 1), representing the deposited density or particle count on the grid. :rtype: ndarray .. rubric:: 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. .. rubric:: References - W. Lyra, private communication (adapted from IDL code by A. Johansen) .. py:function:: 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. :param qpar: The target value for which the nearest index is sought. :type qpar: float :param q: A 1D sorted array of floats (in ascending order) to search within. :type q: array_like :returns: Index of the element in `q` closest to `qpar`. :rtype: int .. rubric:: Notes This method assumes that `q` is sorted in ascending order. If the array is not sorted, the result is undefined. .. rubric:: References - W. Lyra, private communication (adapted from IDL code by A. Johansen) .. py:function:: 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.