pyFLD package

Submodules

pyFLD.FLD module

class pyFLD.FLD.RadiativeEnvironment(fluids, verbose=True, preconditioner='ilu', constant_fluxlimiter=False, implicit_viscosity=False, **kwargs)

Bases: _GenericFluid

Radiation hydrodynamics solver using Flux-Limited Diffusion (FLD).

Collects multiple fluids._RadiativeFluid objects (Gas and Dust) and solves the coupled radiation-hydrodynamics equations using FLD approximation. Handles irradiation, viscous heating, and frequency-dependent opacities. Assumes all fluids share the same grid, star, and temperature field.

Parameters:

fluids (list of _RadiativeFluid) – List of Dust and/or Gas fluid objects to solve
verbose (bool, optional) – Enable verbose output (default: True)
constant_fluxlimiter (bool, optional) – Use constant flux limiter (lambda=1/3) instead of Kley formula (default: False)
preconditioner (str, optional) – Preconditioning method for iterative solver: ‘ilu’, ‘jacobi’, ‘rescaled ilu’, ‘rescaled jacobi’, or False (default: ‘ilu’)
implicit_viscosity (bool, optional) – Treat viscosity implicitly, experimental (default: False)
**kwargs – Arguments passed to fluids._GenericFluid (debug, name)

E0

Reference radiation energy density

Type:: astropy.units.Quantity [erg/cm^3]

Qirr

Irradiation heating term

Type:: astropy.units.Quantity [erg/(s cm^3)]

Qvisc

Viscous heating term

Type:: astropy.units.Quantity [erg/(s cm^3)]

irradiation

Any fluid has irradiation enabled

Type:: bool

viscosity

Any fluid has viscosity enabled

Type:: bool

M

Radiation matrix (7 diagonals for 3D)

Type:: ndarray

B

Right-hand side vector

Type:: astropy.units.Quantity [erg/cm^3]

solve_x1, solve_x2, solve_x3

Whether to solve in each direction

Type:: bool

E0 = <Quantity 7.56573325e-11 erg / cm3>

Qirr = None

Qvisc = None

advance(dt=<Quantity 1. yr>, uf=0)

Advance radiation-hydrodynamics by one timestep.

Complete time stepping procedure: update opacities and effective fractions for all fluids, compute source terms (irradiation, viscous heating), build and solve FLD system, update temperature and synchronize fluids. This is the main function to call in time-stepping loops.

Parameters:

dt (astropy.units.Quantity [s], optional) – Timestep size (default: 1 yr)
uf (float [0-1], optional) – Under-relaxation factor for temperature update (default: 0)

Raises:

ValueError – If any fluid is not initialized

check_grids(): Verify all fluids share the same grid object.

check_source_terms(): Determine which heating source terms are active across all fluids.

check_stars(): Verify all fluids share the same star object.

check_temperatures(): Verify all fluids share the same temperature field.

compute_Qirr()

Compute irradiation heating source term via ray tracing.

Performs ray-tracing from star inward through domain to compute radiation flux at each cell accounting for local opacity. Supports both gray and frequency-dependent irradiation. Computes Q_irr = F x geom_factor where F is attenuated stellar flux.

Accumulates dtau and tau0 across all fluids with irradiation enabled. Ray tracing accounts for cumulative optical depth from inner boundary.

Updates

Qirrastropy.units.Quantity [erg/(s cm^3)]: Total irradiation heating rate
tauastropy.units.Quantity: Optical depth field used in calculation

raises ValueError:: If star object is not available

compute_Qvisc()

Compute viscous heating source term.

Integrates viscous dissipation across all fluids with viscosity enabled: Q_visc = (9/4) alpha sqrt(gamma)^2 P Omega_K. Only works in spherical geometry.

Updates

Qviscastropy.units.Quantity [erg/(s cm^3)]: Total viscous heating rate

raises NotImplementedError:: If grid geometry is not spherical
raises ValueError:: If star object is not available

compute_irradiation_geometry(): Precompute geometric factors for irradiation heating term.

configure_solver(solver='bicgstab', rtol=1e-08, maxiter=10000, guess=True): Configure iterative linear solver parameters.

declare_boundaries(bound, val)

Set boundary conditions for radiation field (FLD).

Specifies boundary type and value for each domain boundary. Arrays must have length 2*dims (6 for 3D): rmin, rmax, thetamin, thetamax, phimin, phimax.

Parameters:

bound (array-like, length 2*dims) – Boundary types: FIXEDVALUE, PERIODIC, REFLECTIVE, UNDEFINED
val (array-like, length 2*dims) – Boundary values for FIXEDVALUE boundaries [erg/cm^3] or Quantity

Raises:

ValueError – If arrays have wrong length or contain unrecognized types

enforce_hydrostatic_equilibrium(rhogas=None, uf=0)

Enforce hydrostatic equilibrium for all fluids.

Solves for vertical structure of all fluids under gravity with pressure and viscous support (dust also settles). Recomputes total density from equilibrium profiles and updates pressure/energy fields.

Parameters:

rhogas (astropy.units.Quantity [g/cm^3], optional) – Gas density field for dust settling. If None, computed from Gas fluids (default: None)
uf (float [0-1], optional) – Under-relaxation factor for convergence (default: 0, full update)

enforce_matrix_boundaries(): Apply boundary conditions to FLD matrix and RHS vector.

get_flux_limiter(R): Compute flux limiter from normalized gradient using Kley (1989) formula.

get_surface_density()

Compute vertically integrated surface density.

Integrates total density in polar angle direction to obtain surface density as function of radius and azimuth. Works only in spherical geometry.

Returns:: Vertically integrated surface density (phi, R)
Return type:: astropy.units.Quantity [g/cm^2]

irradiation = False

iters = -1

prepare_grid()

Initialize matrix structure and determine solution dimensions.

Sets up sparse matrix for FLD equations and determines which directions have resolution > 1 (for diffusion terms). Matrix has 7 diagonals for 3D (center + 6 off-diagonals).

refresh_fluids(): Call refresh() on all fluids to update opacities and effective fractions. Useful to call before visualizing or dumping snapshots.

setup()

Initialize all fluids and compute auxiliary quantities.

Ensures all fluids are initialized, computes total density (summing effective fractions), initializes heating terms (Qirr, Qvisc), and prepares for FLD solving. Shared temperature and grid from fluids[0].

update_fluids()

Synchronize all fluids with updated temperature and radiation energy.

Copies temperature and radiation energy from RadiativeEnvironment to each fluid and recomputes pressure and internal energy.

update_pressure_energy()

Update pressure and internal energy from density and temperature.

Computes pressure using ideal gas law and internal energy density using adiabatic relations. Total density already includes effective fractions (eps) from all fluids.

update_temperature(Enew, dt, uf=0)

Update temperature field from new radiation energy density.

Solves the coupled radiation-matter energy balance equation to obtain temperature. Includes source terms from irradiation and (optionally) implicit viscosity. Applies temperature floor at 3 K to prevent crashes.

Parameters:

Enew (astropy.units.Quantity [erg/cm^3]) – New radiation energy density from FLD solve
dt (astropy.units.Quantity [s]) – Timestep size
uf (float [0-1], optional) – Under-relaxation factor (default: 0, full update)
Updates –
------- –
tmp (astropy.units.Quantity [K]) – Updated temperature field (with under-relaxation)
Er (astropy.units.Quantity [erg/cm^3]) – Updated radiation energy density
prs (astropy.units.Quantity) – Updated pressure and internal energy from new temperature
eng (astropy.units.Quantity) – Updated pressure and internal energy from new temperature

viscosity = False

write_snapshot(dirc='./', name='REnv')

pyFLD.FLD.create_from_snapshot(dirc='./', name='REnv')

pyFLD.debug_state module

pyFLD.debug_state.depth()

pyFLD.debug_state.format_debug(prefix='', msg='', indent=0, log_char='| '): Formats a debug message with indentation.

pyFLD.debug_state.format_warning(prefix='', msg='', indent=0)

pyFLD.debug_state.pop()

pyFLD.debug_state.push()

pyFLD.fluids module

class pyFLD.fluids.Dust(grain_rho=<Quantity 1. g / cm3>, grain_size=<Quantity 1. um>, **kwargs)

Bases: _RadiativeFluid

Dust fluid for protoplanetary disk simulations.

Extends _RadiativeFluid to include dust-specific properties such as grain size and density. Implements dust settling in hydrostatic equilibrium via the Stokes number and dust-to-gas coupling.

Arguments (in addition to `_RadiativeFluid`)

grain_rhofloat [g/cm^3] or astropy.units.Quantity, optional: Material density of dust grains (default: 1 g/cm^3)
grain_sizefloat [cm] or astropy.units.Quantity, optional: Grain size assuming monodisperse distribution (default: 1 μm)

grain_rho

Material density of dust grains

Type:: astropy.units.Quantity [g/cm^3]

grain_size

Grain size

Type:: astropy.units.Quantity [cm]

isDust

Identifier flag (True for Dust)

Type:: bool

isGas

Identifier flag (False for Dust)

Type:: bool

compute_hydrostatic_equilibrium(rhogas=None, uf=0)

Compute dust density profile under hydrostatic equilibrium with settling.

Solves for the vertical dust density structure accounting for gravity, pressure support, and dust-gas interaction via the Stokes number. Dust settles toward the midplane proportional to the Stokes number. Uses a second-order predictor-corrector integration scheme.

Works only in spherical geometry. Enforces correct surface density and prevents density from increasing toward the poles (falling matter).

Parameters:

rhogas (array-like [g/cm^3] or astropy.units.Quantity) – Gas density field required to compute dust-gas interaction and settling. Must have same shape as grid (including ghost cells). (required)
uf (float [0-1], optional) – Under-relaxation factor for convergence. (default: 0)

Returns:

Updated dust density field with boundaries applied

Return type:

ndarray [g/cm^3]

Raises:

ValueError – If rhogas is None, geometry is not spherical, grid structure is invalid, or NaNs appear in result

class pyFLD.fluids.Gas(**kwargs)

Bases: _RadiativeFluid

Gas fluid for protoplanetary disk simulations.

Extends _RadiativeFluid with gas-specific methods. Implements gas hydrostatic equilibrium with vertical structure determined by pressure support and gravity.

isGas

Identifier flag (True for Gas)

Type:: bool

isDust

Identifier flag (False for Gas)

Type:: bool

compute_hydrostatic_equilibrium(rhogas=None, uf=0)

Compute gas density profile under hydrostatic equilibrium.

Solves for the vertical gas density structure from pressure gradient balance. Uses the sound speed and scale height to determine vertical structure. Integrates radial pressure gradients and uses a second-order predictor-corrector scheme for accurate vertical profiles.

Works only in spherical geometry. Enforces correct surface density and prevents density from increasing toward poles (prevents material flowing upward).

Parameters:

rhogas (array-like [g/cm^3] or astropy.units.Quantity, optional) – Not used for gas fluids. Included for API compatibility. (default: None)
uf (float [0-1], optional) – Under-relaxation factor for convergence. (default: 0)

Returns:

Updated gas density field with boundaries applied

Return type:

ndarray [g/cm^3]

Raises:

ValueError – If geometry is not spherical, grid structure is invalid, or NaNs appear in result

pyFLD.fluids.create_from_snapshot(dirc='./', name='fluid')

pyFLD.grid module

class pyFLD.grid.Grid(x1i, x2i, x3i, nghost=2, geometry='cartesian')

Bases: object

3D computational grid for the FLD solver.

Supports Cartesian, cylindrical, and spherical geometries with automatic handling of coordinates and metric factors. Arguments with length units (e.g., x1i) can be provided as astropy quantities or as plain arrays in au. Angles (x2i for cylindrical/spherical, x3i for spherical) can be provided as astropy quantities or as plain arrays in radians. The class will internally convert all length-related inputs to astropy quantities, and strip units from angle-related inputs if they are astropy quantities. This makes it easier to work with the grid class internally while still allowing flexible input formats.

Parameters:

x1i (array_like or astropy.units.Quantity) – 1D array of cell interfaces in direction 1. - For Cartesian: x-coordinate [au]. - For cylindrical: radial distance [au]. - For spherical: radial distance [au].
x2i (array_like or astropy.units.Quantity) – 1D array of cell interfaces in direction 2. For Cartesian: y-coordinate [au]. For cylindrical: azimuthal angle (radians). For spherical: polar angle (radians).
x3i (array_like or astropy.units.Quantity) – 1D array of cell interfaces in direction 3. For Cartesian: z-coordinate [au]. For cylindrical: z-coordinate [au]. For spherical: azimuthal angle (radians).
nghost (int, optional) – Number of ghost cells in each direction (default: 2)
geometry (str, optional) – Grid geometry type: ‘cartesian’, ‘cylindrical’, or ‘spherical’ (default: ‘cartesian’)

x1, x2, x3

Cell center coordinates in each direction

Type:: astropy.units.Quantity (ndarray)

dx1, dx2, dx3

Cell widths in each direction

Type:: astropy.units.Quantity (ndarray)

A1, A2, A3

3D arrays of cell face areas in each direction

Type:: astropy.units.Quantity (ndarray)

dV

3D array of cell volumes

Type:: astropy.units.Quantity (ndarray)

h1, h2, h3

Metric scale factors in each direction (equal to 1 for Cartesian)

Type:: astropy.units.Quantity (ndarray) or float

R, z

2D arrays of cylindrical coordinates (spherical geometry only)

Type:: astropy.units.Quantity (ndarray)

NX1, NX2, NX3

Number of active cells in each direction

Type:: int

NX1_TOT, NX2_TOT, NX3_TOT

Total number of cells (including ghost cells) in each direction

Type:: int

ibeg, jbeg, kbeg

Starting indices of active cells in each direction

Type:: int

iend, jend, kend

Ending indices of active cells in each direction

Type:: int

act1, act2, act3

Slice objects for active cells in each direction

Type:: slice

active

Slice object for all active cells in 3D

Type:: slice

write_snapshot(dirc='./', name='grid')

pyFLD.grid.add_ghost_cells(arr, nghost=2, log=False)

Add ghost cells to a 1D array of cell interfaces.

Ghost cells are appended on both sides of the array and spaced according to the spacing of the outermost cells in the original array.

Parameters:

arr (ndarray) – 1D array of cell interfaces
nghost (int, optional) – Number of ghost cells to add on either side (default: 2)
log (bool, optional) – Use logarithmic spacing for ghost cells (default: False)

Returns:

1D array of cell interfaces with ghost cells added (length Nx+1 + 2*nghost)

Return type:

ndarray

pyFLD.grid.create_from_snapshot(dirc='./', name='grid')

pyFLD.grid.get_domain_with_ghosts(xmin, xmax, Nx=20, log=False, nghost=2)

Create a 1D array of cell interfaces with ghost cells.

Generates Nx+1 cell interfaces from xmin to xmax with optional logarithmic spacing and appends ghost cells on both sides.

Parameters:

xmin (float) – Minimum value of the domain (left edge of first cell)
xmax (float) – Maximum value of the domain (right edge of last cell)
Nx (int, optional) – Number of cells in the domain (default: 20)
log (bool, optional) – Use logarithmic spacing (default: False)
nghost (int, optional) – Number of ghost cells to append on either side (default: 2)

Returns:

1D array of cell interfaces including ghost cells

Return type:

ndarray

pyFLD.radmc_helper module

Radiative transfer helper functions for RADMC3D integration.

A collection of tools to streamline setting up and analyzing RADMC3D models when used alongside radmc3dPy. Provides utilities for managing stars, grids, opacities, and density/temperature fields in a standardized format.

Typical Workflow

Set up a star and grid using radmc_helper.Star and radmc_helper.Grid.
Define a volume density field or use helper functions with surface density and temperature profiles.
Generate dust opacities using OpTool or equivalent.
Use write_*() functions to output RADMC3D-format files.
Run RADMC3D (mctherm/image) to compute radiative transfer.
Read results with radmc3dPy.analyze.readData() and process with helper functions.

See also

star.Star: Base star class with stellar parameters
write: Output star parameters to RADMC3D stars.inp file
read_from_file: Load star parameters from RADMC3D stars.inp file

read_from_file(filename=None)

Read star parameters from RADMC3D stars.inp file.

Parses a RADMC3D format star file and initializes stellar parameters (mass, radius, temperature, luminosity, spectral data).

Parameters:: filename (str, optional) – Input filename. If None, uses the filename from initialization. Default is None.
Returns:: Initializes parent class via super().__init__() with extracted parameters.
Return type:: None
Raises:: ValueError – If file format version is not 2. If number of stars is not 1.

Notes

Supports RADMC3D format version 2
Only single-star files are supported
Temperature is stored as negative in file, converted to positive
Wavelengths converted from microns to astropy units
Luminosity computed from Stefan-Boltzmann law

Examples

>>> rmc_star = radmc_helper.Star()
>>> rmc_star.read_from_file('stars.inp')

setup_from(star: Star)

Initialize this star from another Star object.

Copies all attributes (mass, radius, temperature, luminosity, wavelengths, etc.) from the source star to this instance.

Parameters:: star (star.Star) – Source star object to copy parameters from.
Returns:: Modifies instance in place.
Return type:: None

Examples

>>> source_star = star.Star(wl=np.logspace(-1, 4, 100)*u.um)
>>> radmc_star = radmc_helper.Star()
>>> radmc_star.setup_from(source_star)

update_wavelengths(wl)

Update wavelength grid and recompute spectral luminosity.

Replaces the wavelength array and recalculates luminosity in each wavelength bin based on stellar spectrum.

Parameters:: wl (array_like [cm] or astropy.units.Quantity) – New wavelength array. Can be astropy quantity or numpy array (interpreted as cm if plain array).
Returns:: Modifies instance attributes (wl, nbins, spectral_L) in place.
Return type:: None

Notes

Wavelengths are stored in cm internally
Spectral luminosity recomputed assuming Planck spectrum
Uses methods get_and_set_spectral_L() and _enforce_L_consistency()
Total luminosity should remain consistent

Examples

>>> rmc_star = radmc_helper.Star(wl=np.logspace(-1, 4, 50)*u.um)
>>> new_wl = np.logspace(-1, 4, 200)*u.um  # Finer grid
>>> rmc_star.update_wavelengths(new_wl)

write(flnm=None)

Write star parameters to RADMC3D stars.inp file.

Outputs stellar mass, radius, temperature, position, and spectral information in RADMC3D format. Currently supports single-star systems.

Parameters:: flnm (str, optional) – Output filename. If None, uses the filename from initialization. Default is None.
Returns:: Writes to file as a side effect.
Return type:: None

Notes

Format: RADMC3D format version 2
Temperature is written as negative value in file
Wavelengths are stored for each spectral element
Only one star per file is currently supported

Examples

>>> rmc_star = Star(M=1*u.Msun, T=5772*u.K, R=1*u.Rsun, wl=np.logspace(-1, 4, 100)*u.um)
>>> rmc_star.write('stars.inp')

pyFLD.radmc_helper.get_wavelength_array(wl_min=<Quantity 0.1 um>, wl_max=<Quantity 1. cm>, Npts=100, refine_10um=True, Tstar=<Quantity 5772. K>)

Generate a logarithmically-spaced wavelength array for radiative transfer.

Creates a wavelength grid suitable for RADMC3D simulations, with optional refinement around the 10 um silicate feature and validation of stellar peak coverage.

Parameters:

wl_min (float or astropy.units.Quantity, optional) – Minimum wavelength [um]. Default is 0.1 um.
wl_max (float or astropy.units.Quantity, optional) – Maximum wavelength [um]. Default is 1 cm.
Npts (int, optional) – Number of wavelength points (per segment if refine_10um=True). Default is 100.
refine_10um (bool, optional) – If True, refines grid around 10 um silicate feature using 2*Npts points in the range [7, 25] um. Default is True.
Tstar (float or astropy.units.Quantity, optional) – Stellar temperature [K] used to estimate emission peak. Default is 5772 K (Sun).

Returns:

wl – Wavelength array [um], logarithmically spaced.

Return type:

astropy.units.Quantity

Warns:

UserWarning – If wl_min is greater than 1/4 of the stellar emission peak wavelength, indicating insufficient resolution of stellar spectrum.

Notes

If refine_10um=True, grid has 3 segments:
- [wl_min, 7 um]: Npts//2 points
- [7 um, 25 um]: 2*Npts points (refined around 10 um feature)
- [25 um, wl_max]: Npts//2 points
Stellar emission peak computed using Wien displacement law: lambda_peak = 2898 um*K / T
Useful for dust continuum radiative transfer with silicate absorption

Examples

>>> wl = get_wavelength_array(wl_min=0.1*u.um, wl_max=10*u.mm, Npts=50)
>>> print(f"Total points: {len(wl)}")
>>> # Output: Total points: 175 (50//2 + 2*50 + 50//2)

>>> # Without refinement (simple log grid)
>>> wl_simple = get_wavelength_array(Npts=200, refine_10um=False)
>>> print(f"Total points: {len(wl_simple)}")
>>> # Output: Total points: 200

pyFLD.radmc_helper.read_and_format_radmc_data(dirc='./')

Read and format RADMC3D output files using radmc3dPy.

Parses RADMC3D mctherm output (density, temperature, grid) and returns formatted astropy quantities with proper dimensional units.

Parameters:

dirc (str, optional) – Directory containing RADMC3D output files. Default is './'.

Returns:

rho (astropy.units.Quantity) – Dust volume density array with shape (Nr_species, Nr, Ntheta, Nphi) Units: [g/cm3].
tmp (astropy.units.Quantity) – Dust temperature array with shape (Nr_species, Nr, Ntheta, Nphi) Units: [K].
Requires
——–
radmc3dPy (External package for RADMC3D data I/O) – Must be installed and importable.

Raises:

ImportError – If radmc3dPy is not available.

Notes

Expected files in directory: amr_grid.inp, dust_density.binp, dust_temperature.bdat
Array axis order transposed to (Nr_species, Nr, Ntheta, Nphi)
Current working directory temporarily changed during file read, then restored
Requires successful RADMC3D mctherm run

Examples

>>> # After RADMC3D mctherm run
>>> rho, tmp = read_and_format_radmc_data(dirc='./radmc_output/')
>>> print(f"Density shape: {rho.shape}")  # (Nr_species, Nr, Ntheta, Nphi)
>>> print(f"Temperature range: {tmp.min():.1f} to {tmp.max():.1f}")
>>> print(f"Density units: {rho.unit}")

pyFLD.radmc_helper.read_wavelength_array(dirc='./', flnm='wavelength_micron.inp')

Read wavelength array from RADMC3D format file.

Parses a wavelength_micron.inp file and returns the wavelength array.

Parameters:

dirc (str, optional) – Input directory. Default is './'.
flnm (str, optional) – Input filename. Default is 'wavelength_micron.inp'.

Returns:

wl – Wavelength array [um].

Return type:

astropy.units.Quantity

Notes

First line: number of wavelength points (Npts)
Following Npts lines: wavelength values in um
Expects space or newline-separated format

Examples

>>> wl = read_wavelength_array(dirc='./', flnm='wavelength_micron.inp')
>>> print(f"Wavelength range: {wl[0]:.2e} to {wl[-1]:.2e}")
>>> print(f"Number of points: {len(wl)}")

pyFLD.radmc_helper.write_control_file(dirc='./', flnm='radmc3d.inp', extra=True, **kwargs)

Write RADMC3D control file with simulation parameters.

Generates a radmc3d.inp configuration file with specified parameters. Includes sensible defaults for common simulations if extra=True.

Parameters:

dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'radmc3d.inp'.
extra (bool, optional) – If True, adds helpful default parameters for dust continuum radiative transfer. Default is True.
**kwargs – Custom parameters to write as key=value pairs. Overrides defaults if extra=True.

Returns:

Writes to file as a side effect.

Return type:

None

Notes

Default parameters (if extra=True):

incl_dust=1: Include dust continuum
incl_lines=0: Exclude line radiative transfer
nphot=1e9: Number of photon packages (mctherm/image)
nphot_scat=1e8: Photon packages for scattering
nphot_spec=1e8: Packages for spectrum
nphot_mono=1e8: Packages for monochromatic
rto_style=3: Use Monte Carlo (Bjorkman & Wood)
istar_sphere=1: Treat star as sphere
setthreads=4: OpenMP threads
tgas_eq_tdust=1: Gas temp = dust temp
modified_random_walk=1: Use modified random walk for absorption
scattering_mode_max=3: Maximum scattering order
iranfreqmode=1: Frequency-dependent mode
mc_scat_maxtauabs=5.0: Maximum optical depth for scattering

Examples

>>> # Default parameters (good for most simulations)
>>> write_control_file(dirc='./')

>>> # Custom parameters
>>> write_control_file(dirc='./', nphot=int(1e10), setthreads=8)

>>> # No extra defaults but custom parameters
>>> write_control_file(dirc='./', extra=False, nphot=int(1e8))

pyFLD.radmc_helper.write_density_array(rho=[], dirc='./', flnm='dust_density.binp')

Write dust density array(s) to RADMC3D binary format file.

Outputs volume density for one or more dust species to dust_density.binp in RADMC3D binary format.

Parameters:

rho (list of array_like [g/cm3] or astropy.units.Quantity) – List of dust density arrays (one per species). Each array can be astropy quantity or numpy array (interpreted as g/cm3 if plain array). Default is empty list.
dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'dust_density.binp'.

Returns:

Writes to binary file as a side effect.

Return type:

None

Notes

Binary format: C-style (little-endian double precision)
Header format: [iformat=1, precision=8, Ncells, Nr_species]
All density arrays must have same shape (Ncells)
Single species: pass as [rho], not rho
Units converted to g/cm3 before writing

Examples

>>> # Single dust species
>>> rho = 1e-11 * u.g/u.cm**3 * np.ones((100, 50, 30))
>>> write_density_array([rho], dirc='./', flnm='dust_density.binp')

>>> # Multiple species (same spatial grid)
>>> rho1 = 1e-11 * u.g/u.cm**3 * np.ones((100, 50, 30))
>>> rho2 = 5e-12 * u.g/u.cm**3 * np.ones((100, 50, 30))
>>> write_density_array([rho1, rho2], dirc='./')

pyFLD.radmc_helper.write_flat_opacities(wl, kabs, dirc='./', flnm='dustkappa_flat.inp')

Write dust opacity file in RADMC3D format.

Outputs a wavelength-dependent dust opacity to dustkappa_*.inp in RADMC3D format. Includes absorption and scattering opacities (scattering set to zero, phase function albedo to zero).

Note

RADMC3D expects wavelengths in microns in the output file. This function interprets plain arrays as cm by default (following pyFLD convention), but properly handles astropy quantities with any length unit.

All three calls produce output in microns but interpret input differently:

wl = np.array([1, 2, 10]) in dubious units
write_flat_opacities(wl) interprets wl as cm
write_flat_opacities(wl*u.um) uses wl as is (in um)
write_flat_opacities(wl*u.au) converts from au to um

Parameters:

wl (array_like [cm] or astropy.units.Quantity) – Wavelength array. Can be astropy quantity with any length unit, or numpy array (interpreted as cm if plain array). Output will be in um.
kabs (array_like [cm2/g] or astropy.units.Quantity) – Absorption opacity (mass absorption coefficient). Can be astropy quantity or numpy array (interpreted as cm2/g if plain array). Must have same shape as wl. Scattering and albedo set to zero.
dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'dustkappa_flat.inp'.

Returns:

Writes to file as a side effect.

Return type:

None

Notes

Format: RADMC3D opacity file format version 3
Each line: wavelength [um], k_abs [cm2/g], k_scat [cm2/g], phase_albedo
Scattering opacity set to 0 (absorption-only)
Phase function albedo set to 0
Compatible with OpTool and RADMC3D

Examples

>>> # Wavelength-dependent absorption (e.g., silicate)
>>> wl = np.logspace(-1, 4, 100)*u.um
>>> kabs = 1e-3 * u.cm**2/u.g * np.exp(-((wl.to_value('um') - 10)/5)**2)
>>> write_flat_opacities(wl, kabs, dirc='./', flnm='dustkappa_silicate.inp')

>>> # Constant opacity
>>> kabs_const = np.ones_like(wl) * 1e-2 * u.cm**2/u.g
>>> write_flat_opacities(wl, kabs_const, dirc='./')

pyFLD.radmc_helper.write_opacity_handle(names=[], dirc='./', flnm='dustopac.inp')

Write dust opacity file handles to RADMC3D format file.

Outputs opacity file references to dustopac.inp in RADMC3D format. References OpTool-generated or custom dustkappa_<name>.inp files.

Parameters:

names (list of str) – Opacity file handles (names after dustkappa_ prefix). Example: [‘silicate’, ‘graphite’] -> dustkappa_silicate.inp, dustkappa_graphite.inp Default is empty list.
dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'dustopac.inp'.

Returns:

Writes to file as a side effect.

Return type:

None

Notes

Format: RADMC3D dustopac.inp version 2
Each opacity assumes thermal grain (no aligned grains)
File names expanded to full dustkappa_<name>.inp format
Compatible with RADMC3D and OpTool

Examples

>>> names = ['1um', '10um', '100um']
>>> write_opacity_handle(names, dirc='./', flnm='dustopac.inp')
>>> # Creates references to:
>>> # dustkappa_1um.inp
>>> # dustkappa_10um.inp
>>> # dustkappa_100um.inp

pyFLD.radmc_helper.write_temperature_array(tmp=[], dirc='./', flnm='dust_temperature.bdat')

Write dust temperature array(s) to RADMC3D binary format file.

Outputs dust temperature for one or more dust species to dust_temperature.bdat in RADMC3D binary format. Compatible with RADMC3D mctherm output.

Parameters:

tmp (list of array_like [K] or astropy.units.Quantity) – List of dust temperature arrays (one per species). Each array can be astropy quantity or numpy array (interpreted as K if plain array). Default is empty list.
dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'dust_temperature.bdat'.

Returns:

Writes to binary file as a side effect.

Return type:

None

Notes

Binary format: C-style (little-endian double precision)
Header format: [iformat=1, precision=8, Ncells, Nr_species]
All temperature arrays must have same shape (Ncells)
Single species: pass as [tmp], not tmp
Units converted to K before writing

Examples

>>> # Single dust species at thermal equilibrium
>>> T_dust = 100 * u.K * np.ones((100, 50, 30))
>>> write_temperature_array([T_dust], dirc='./', flnm='dust_temperature.bdat')

>>> # Multiple species (with different temperatures)
>>> T_small = 150 * u.K * np.ones((100, 50, 30))
>>> T_large = 80 * u.K * np.ones((100, 50, 30))
>>> write_temperature_array([T_small, T_large], dirc='./')

pyFLD.radmc_helper.write_wavelength_array(wl, dirc='./', flnm='wavelength_micron.inp')

Write wavelength array to RADMC3D format file.

Outputs a wavelength grid to wavelength_micron.inp in RADMC3D format.

Parameters:

wl (array_like [um] or astropy.units.Quantity) – Wavelength array. Can be astropy quantity or numpy array (interpreted as um if plain array).
dirc (str, optional) – Output directory. Default is './'.
flnm (str, optional) – Output filename. Default is 'wavelength_micron.inp'.

Returns:

Writes to file as a side effect.

Return type:

None

Notes

First line: number of wavelength points
Following lines: wavelength values in um, one per line
Values written in scientific notation with 18 decimal places
Compatible with RADMC3D wavelength_micron.inp format

Examples

>>> wl = np.logspace(-1, 4, 100)*u.um
>>> write_wavelength_array(wl, dirc='./', flnm='wavelength_micron.inp')
>>> # File format:
>>> # 100
>>> # 1.000000000000000e-01
>>> # ...

pyFLD.star module

class pyFLD.star.Star(M=<Quantity 1. solMass>, L=<Quantity 1. solLum>, T=<Quantity 5772. K>, wl=None)

Bases: object

Stellar parameters for radiative irradiation calculations.

Holds stellar mass, luminosity, effective temperature, and optionally a frequency-dependent spectrum. Computes stellar radius from the Stefan-Boltzmann law.

Parameters:

M (float [Msun] or astropy.units.Quantity, optional) – Stellar mass (default: 1 Msun)
L (float [Lsun] or astropy.units.Quantity, optional) – Stellar luminosity (default: 1 Lsun)
T (float [K] or astropy.units.Quantity, optional) – Stellar effective temperature (default: 5772 K)
wl (array-like [cm] or astropy.units.Quantity, optional) – Wavelengths for spectrum sampling in cm (default: None, no spectrum)

M

Stellar mass in solar masses

Type:: astropy.units.Quantity

L

Stellar luminosity in solar luminosities

Type:: astropy.units.Quantity

T

Stellar effective temperature in Kelvin

Type:: astropy.units.Quantity

R

Stellar radius in solar radii (computed from L and T)

Type:: astropy.units.Quantity

GM

Gravitational parameter G * M

Type:: astropy.units.Quantity

wl

Wavelengths for spectrum sampling in cm (set if provided at initialization)

Type:: astropy.units.Quantity, optional

nbins

Number of wavelength bins (set if wl is provided)

Type:: int, optional

L_wl

Stellar luminosity at each wavelength such that sum(L_wl) = L (set if wl is provided)

Type:: astropy.units.Quantity, optional

L_wl = None

get_and_set_spectral_L()

Compute and store the frequency-dependent spectrum.

Calculates spectral luminosity from the Planck function and ensures consistency with the total luminosity through normalization.

get_spectral_L()

Compute frequency-dependent stellar luminosity from Planck function.

Calculates spectral luminosity at each wavelength using the Planck function and normalizes such that the integrated sum equals the total luminosity L.

Returns:: Stellar luminosity at each wavelength in solar luminosities
Return type:: astropy.units.Quantity

get_wavelength_spacing()

Compute wavelength spacing array for spectrum calculations.

Assumes wavelengths are logarithmically spaced and computes spacing in logarithmic space. Works with any spacing pattern.

Returns:: Wavelength spacing array in cm
Return type:: astropy.units.Quantity
Raises:: ValueError – If wavelength array wl is not defined for this star

nbins = None

wl = None

write_snapshot(dirc='./', name='star')

pyFLD.star.create_from_snapshot(dirc='./', name='star')

pyFLD.tools module

pyFLD.tools.BBflux(T=<Quantity 5780. K>, wl=<Quantity 0.0001 cm>)

Compute black body flux spectrum and its temperature derivative.

Calculates the Planck function and its derivative with respect to temperature using the standard black body formula:

f = exp(h*c/(λkT))
B(λ) = 2hc²/λ⁵ / (f - 1)
U(λ) = 2hc²/λ⁵ * (h*c/(λkT²)) * f / (f - 1)²

Parameters:

T (float [K] or astropy.units.Quantity, optional) – Temperature of the black body (default: 5780 K, solar temperature)
wl (float [cm] or astropy.units.Quantity, optional) – Wavelengths to sample (default: 1e-4 cm)

Returns:

B (astropy.units.Quantity [erg/(cm^3 s)]) – Black body flux at requested wavelengths
U (astropy.units.Quantity [erg/(cm^3 s K)]) – Derivative of black body flux with respect to temperature

pyFLD.tools.Er_to_T(Er)

Convert radiation energy density to temperature.

Inverts the relation Er = aR * T^4 to solve for T.

Parameters:: Er (float [erg/cm^3] or astropy.units.Quantity) – Radiation energy density
Returns:: Temperature
Return type:: astropy.units.Quantity [K]

pyFLD.tools.T_to_Er(T)

Convert temperature to radiation energy density.

Uses the relation Er = aR * T^4 where aR is the radiation constant.

Parameters:: T (float [K] or astropy.units.Quantity) – Temperature
Returns:: Radiation energy density
Return type:: astropy.units.Quantity [erg/cm^3]

pyFLD.tools.assign_dictionary(obj, d)

pyFLD.tools.atleast_nd(a, d=3)

Extend array to at least d dimensions by appending new axes.

Similar to np.atleast_1d/2d/3d but works for arbitrary dimensions. New dimensions are always appended at the end.

Parameters:

a (array-like or scalar) – Input array or scalar
d (int, optional) – Desired number of dimensions (default: 3)

Returns:

Input array reshaped to have at least d dimensions

Return type:

ndarray

pyFLD.tools.compute_mean_opacities(wl, kabs, ksca, g, T=None, Tmin=<Quantity 0.1 K>, Tmax=<Quantity 100000. K>, nT=1000)

Compute Rosseland and Planck mean opacities from wavelength-dependent opacities.

Integrates absorption and scattering opacities weighted by the Planck and inverse-mean Rosseland functions across wavelength to obtain temperature-dependent mean opacities.

Parameters:

wl (array-like [cm] or astropy.units.Quantity) – Wavelengths to sample on
kabs (array-like [cm^2/g] or astropy.units.Quantity) – Absorption opacity at each wavelength
ksca (array-like [cm^2/g] or astropy.units.Quantity) – Scattering opacity at each wavelength
g (array-like) – Asymmetry parameter at each wavelength
T (array-like [K] or astropy.units.Quantity, optional) – Temperature array at which to sample mean opacities. If provided, overrides Tmin, Tmax, and nT (default: None)
Tmin (float [K] or astropy.units.Quantity, optional) – Minimum temperature for mean opacity computation (default: 0.1 K)
Tmax (float [K] or astropy.units.Quantity, optional) – Maximum temperature for mean opacity computation (default: 1e5 K)
nT (int, optional) – Number of temperature points to compute (default: 1000)

Returns:

T (astropy.units.Quantity [K]) – Array of temperatures at which mean opacities were computed
kR (astropy.units.Quantity [cm^2/g]) – Rosseland mean opacity at each temperature
kP (astropy.units.Quantity [cm^2/g]) – Planck mean opacity at each temperature

pyFLD.tools.extract_and_store_params(obj, dirc='./', skipped_keys=[])

pyFLD.tools.extract_dictionary(obj, skipped_keys=[])

pyFLD.tools.flux_limiter_kley(Rrad)

Compute radiation flux limiter using the Kley (1989) formula.

The flux limiter is a function of the normalized radiation gradient R that transitions between the free-streaming limit (large R) and diffusion limit (small R).

Parameters:: Rrad (float or array-like) – Normalized radiation gradient (dimensionless)
Returns:: Flux limiter value(s) at requested R
Return type:: float or ndarray

pyFLD.tools.get_Wien_temperature(wl)

Compute temperature from peak wavelength (Wien’s displacement law).

Inverts Wien’s displacement law to determine the black body temperature corresponding to a given peak emission wavelength.

Parameters:: wl (float [cm] or astropy.units.Quantity) – Peak emission wavelength
Returns:: Black body temperature
Return type:: astropy.units.Quantity [K]

pyFLD.tools.get_Wien_wavelength(T)

Compute peak wavelength of black body emission (Wien’s displacement law).

Parameters:: T (float [K] or astropy.units.Quantity) – Temperature
Returns:: Peak wavelength at which black body emission is maximum
Return type:: astropy.units.Quantity [um]

pyFLD.tools.get_alpha_Cecil_Flock_2024(fluid)

Compute effective alpha viscosity parameter following Cecil & Flock (2024).

Calculates a temperature-dependent alpha viscosity parameter that transitions smoothly from 1e-3 at low temperatures to 0.1 at high temperatures, with a transition centered around 900 K and a width of 25 K.

Parameters:: fluid (GenericFluid or higher class) – Fluid object with temperature field
Returns:: Effective alpha viscosity parameter across the domain
Return type:: ndarray

pyFLD.tools.get_capped(arr, side='both', log=False)

Return a copy of input array with values capped to prevent overflow/underflow.

Parameters:

arr (array-like) – Input array to cap
side (str, optional) – Which side to cap: ‘both’, ‘min’, or ‘max’ (default: ‘both’)
log (bool, optional) – Whether input array is in log space (default: False)

Returns:

Copy of input array with values capped

Return type:

ndarray

Raises:

ValueError – If ‘side’ argument is not ‘both’, ‘min’, or ‘max’

pyFLD.tools.get_col0_Flock_2019(fluid)

Compute fluid column density at inner boundary following Flock et al. (2019).

Calculates the column density at the inner edge of the domain using the relation Σ = ρ * ε * (R - 3*R_*).

Parameters:: fluid (GenericFluid or higher class) – Fluid object with density and effective fraction fields
Returns:: Fluid column at innermost part of domain (φ, θ, ibeg)
Return type:: ndarray

pyFLD.tools.get_fraction_Isella_Natta_2005(fluid)

Compute dust sublimation fraction following Isella & Natta (2005).

Calculates the effective dust fraction accounting for sublimation at high temperatures using a smooth transition function. The sublimation temperature depends on dust density via T_subl = 2000 * ρ^0.0195 K.

Parameters:: fluid (GenericFluid or higher class) – Fluid object with temperature and density fields. Ideally Dust.
Returns:: Dust sublimation fraction (0 to 1) across the domain
Return type:: ndarray

pyFLD.tools.get_gas_Planck_opacity_Freedman_2008(fluid)

Compute gas Planck mean opacity following Freedman et al. (2008).

Simple temperature-dependent parameterization of gas opacity based on fitting functions from Freedman et al. (2008).

Parameters:: fluid (GenericFluid or higher class) – Fluid object with temperature field. Ideally Gas.
Returns:: Planck mean opacity at each location in the domain
Return type:: astropy.units.Quantity [cm^2/g]

pyFLD.tools.get_gas_Rosseland_opacity_Freedman_2008(fluid)

Compute gas Rosseland mean opacity following Freedman et al. (2008).

Density-dependent parameterization of gas opacity based on fitting functions from Freedman et al. (2008).

Parameters:: fluid (GenericFluid or higher class) – Fluid object with density field. Ideally Gas.
Returns:: Rosseland mean opacity at each location in the domain
Return type:: astropy.units.Quantity [cm^2/g]

pyFLD.tools.get_one(fluid)

Return a unit value with appropriate shape and units for a fluid.

Helper function for enrolling custom functions that compute fluid properties (e.g., opacity, viscosity) with proper dimensionality.

Parameters:: fluid (RadiativeFluid) – Fluid object (used for context)
Returns:: Unit value (shape and units handled by broadcasting)
Return type:: float

pyFLD.tools.get_parameters(obj, skipped_keys=[])

pyFLD.tools.get_surface_density(grid: Grid, rho)

Compute vertically integrated surface density.

Integrates volume density in the polar angle direction to obtain surface density as a function of radius and azimuth.

Parameters:

grid (pyFLD.grid.Grid) – Computational grid object (spherical geometry only)
rho (astropy.units.Quantity [g/cm^3]) – Volume density field with shape matching grid (including ghost cells)

Returns:

Vertically integrated surface density (φ, R)

Return type:

astropy.units.Quantity [g/cm^2]

Raises:

NotImplementedError – If grid geometry is not spherical

pyFLD.tools.get_tau0_Flock_2019(fluid)

Compute optical depth to inner boundary following Flock et al. (2019).

Calculates the optical depth at the inner edge of the domain using the relation τ = ρ * ε * κ_P,* * (R - 3*R_*), where κ_P,* is the Planck mean opacity at stellar temperature.

Parameters:: fluid (RadiativeFluid or higher class) – Fluid object with density, effective fraction, and opacity fields
Returns:: Optical depth at innermost part of domain (φ, θ, ibeg)
Return type:: ndarray

pyFLD.tools.get_zero(fluid)

Return a zero value with appropriate shape and units for a fluid.

Helper function for enrolling custom functions that compute fluid properties (e.g., opacity, viscosity) with proper dimensionality.

Parameters:: fluid (RadiativeFluid) – Fluid object (used for context)
Returns:: Zero value (shape and units handled by broadcasting)
Return type:: float

pyFLD.tools.have_same_class_contents(cls1, cls2)

Compare instance attributes and methods of two class objects.

Checks if two class instances have identical instance-level attributes and methods, excluding built-in methods. Useful for validating that objects in a collection have consistent state.

Parameters:

cls1 (object) – First object to compare
cls2 (object) – Second object to compare

Returns:

True if both objects have identical instance contents, False otherwise

Return type:

bool

pyFLD.tools.isQuantity(q)

Check if input is an astropy Quantity.

Parameters:: q (any) – Input object to check
Returns:: True if q is an astropy.units.Quantity, False otherwise
Return type:: bool

pyFLD.tools.is_halfdisk(grid: Grid)

Determine if grid covers only one hemisphere in spherical geometry.

Checks if the grid’s polar angle range is [0, π/2] or [π/2, π], indicating a half-disk domain.

Parameters:: grid (pyFLD.grid.Grid) – Computational grid object
Returns:: True if grid covers one hemisphere, False otherwise
Return type:: bool
Raises:: NotImplementedError – If grid geometry is not spherical

pyFLD.tools.mkdir(name)

pyFLD.tools.parse_optool_file(filename)

Read and parse an optool opacity file.

Extracts wavelength-dependent absorption opacity, scattering opacity, and asymmetry parameter from an optool-formatted file and converts to astropy Quantity objects.

Parameters:

filename (str) – Path to the optool file (typically named “dustkappa*.inp”)

Returns:

wl (astropy.units.Quantity [cm]) – Wavelength sampling points
kabs (astropy.units.Quantity [cm^2/g]) – Absorption opacity at each wavelength
ksca (astropy.units.Quantity [cm^2/g]) – Scattering opacity at each wavelength
g (ndarray) – Asymmetry parameter at each wavelength (Henyey & Greenstein 1941)

Raises:

FileNotFoundError – If the file does not exist
ValueError – If the file format is invalid or cannot be parsed

pyFLD.tools.read_and_assign_params(obj, dirc='./')

pyFLD.tools.read_dictionary(dirc='./')

pyFLD.tools.toArray(q, unit='')

Convert astropy Quantity to numpy array (removes units).

If input is a Quantity, converts to specified unit and returns array. If input is not a Quantity, converts to numpy array.

Parameters:

q (array-like or astropy.units.Quantity) – Input Quantity or array
unit (str or astropy.units.Unit, optional) – Target unit for conversion (default: ‘’)

Returns:

Input as numpy array without units

Return type:

ndarray

Raises:

astropy.units.UnitConversionError – If input Quantity has incompatible units

pyFLD.tools.toQuantity(arr, unit)

Convert array-like input to an astropy Quantity with specified unit.

If input is already a Quantity, converts to the target unit. If input is not a Quantity, multiplies by the specified unit.

Parameters:

arr (array-like or astropy.units.Quantity) – Input array or Quantity
unit (str or astropy.units.Unit) – Target unit for the result

Returns:

Input converted to specified unit

Return type:

astropy.units.Quantity

Raises:

astropy.units.UnitConversionError – If arr is a Quantity with incompatible units

pyFLD.tools.try_reading(filename, unit=None, shape=None, force=False)

pyFLD.tools.try_writing(q, filename, unit=None)

pyFLD.tools.write_dictionary(d, dirc='./')

pyFLD package

Submodules

pyFLD.FLD module

Updates

Updates

pyFLD.debug_state module

pyFLD.fluids module

Arguments (in addition to _RadiativeFluid)

pyFLD.grid module

pyFLD.radmc_helper module

Typical Workflow

pyFLD.star module

pyFLD.tools module

Arguments (in addition to `_RadiativeFluid`)