The voxel solver separates three processes that are easy to conflate:
interception removes radiation from a ray according to voxel geometry, PAD, the projection coefficient G, and path length;
scattering re-emits a fraction of intercepted radiation into new directions;
absorption is the remaining terminal fraction retained by vegetation.
This distinction matters because intercepted energy can be intercepted again after scattering.
For each waveband, voxel foliage has leaf reflectance \(\rho\) and leaf transmittance \(\tau\). The single-scattering albedo and absorptance are
All four quantities are dimensionless. ArchimedLight enforces \(0 \leq \rho \leq 1\), \(0 \leq \tau \leq 1\), and \(\rho + \tau \leq 1\) at construction. Set reflectance and transmittance explicitly in VoxelMaterial and supply those materials through models_for.
Leaf transmittance is not canopy porosity. Gaps are already represented by PAD through Beer–Lambert attenuation:
An empty voxel therefore intercepts nothing, irrespective of values stored in its optical arrays. Conversely, optical transmittance acts only after plant material has intercepted radiation.
VoxelOpticalProperties accepts scalar coefficients without allocating dense fields, or arrays aligned exactly with grid.pad. For a plain VoxelGrid, these are effective properties of the material represented in each cell.
The individual/organ workflow instead resolves explicit VoxelMaterial definitions and component proxy densities. With constant component projection coefficients, the existing isotropic solver is reused through total extinction and extinction-weighted band optics:
This is exact for the declared homogeneous mixture and constant projection model; it does not resolve organ positions or directional phase functions. Component interception and absorption are recovered with the same extinction weights, rather than by weighting the final absorbed total by area alone. See Interception And Assumptions.
The optical convention is:
PAR: 400–700 nm;
NIR: 700–2500 nm.
Voxel radiation outputs use energy, power or irradiance units, so spectra must be reduced with incident spectral irradiance weighting. For a band \(b\):
Do not mix photon-weighted PAR coefficients with energy-weighted model outputs. If the source data are photon spectral densities or PPFD, convert them to energy spectral density before averaging for this solver.
Meteorological NIR must match the optical convention
When only total shortwave radiation is available, ArchimedLight can derive its operational NIR forcing as the non-PAR fraction of shortwave radiation. That broad forcing is not guaranteed to be an exact 700–2500 nm integral. For calibrated spectral work, supply Ri_PAR_f and Ri_NIR_f derived over the same intervals used for the optical coefficients.
Useful coefficient sources include species-specific integrating-sphere measurements, the LOPEX93 database, the ANGERS leaf optical database, and spectra simulated with PROSPECT-D. Store the species, leaf state, measurement geometry, wavelength grid, averaging spectrum, and date with calibrated inputs. The convenience constructor voxel_generic_green_leaf_optics (0.17 PAR, 0.85 NIR) supplies teaching values, not species calibrations.
First-order rays begin at the canopy boundary. Scattered radiation instead starts inside an occupied voxel, so it uses a dedicated transport operator. For each occupied source voxel, ArchimedLight traces the full-sphere angular quadrature from the voxel centre. This centre-source approximation is deterministic, includes the remaining path through the emitting voxel, and therefore permits within-cell re-interception. It does not integrate over all possible sub-voxel emission positions; convergence with voxel refinement should be checked when this approximation matters.
The angular quadrature pairs every diffuse turtle direction with its opposite. The normalized weights sum to one, and each hemisphere receives one half of an isotropic foliage source. Sun sectors are excluded because scattered emission is independent of the instantaneous direct-beam direction.
The matrix-free operator applies Beer–Lambert attenuation along cached DDA paths. The step and series workflows prepare compact paths: empty-voxel segments are discarded, and each retained segment stores its linear voxel index and the same interception/transmission fractions used by geometric tracing. Their original order and terminal soil or escape destination are preserved. These fractions depend on PAD, G and path length, so the existing cache validation rejects changed extinction fields while allowing new leaf albedos or soil reflectance. For every application:
with ground absorption replacing bottom escape when a ground boundary is present. A tiny dense exchange matrix exists only as a test oracle; production code never materializes an \(N_\mathrm{voxel}^2\) matrix.
For scattering order \(n\) and one band:
The solver stops when the newly intercepted order is smaller than scattering_stop_ratio times first-order foliage interception, or at the floating-point floor. scattering_max_iter is a hard limit. If that limit is reached first, converged is false and the latest unprocessed order remains in unresolved_energy; it is never silently treated as absorption. The convergence decision also includes unresolved soil-to-soil reflection after the terrain transport's 128-bounce limit. Exhausting foliage orders alone does not establish convergence on a highly reflective concave terrain.
The terminal energy invariant is checked independently for PAR and NIR:
total_intercepted_energy and scattered_energy are intentionally absent from the right-hand side: they describe internal transfers that would otherwise be double-counted.
The result retains the scalar intercepted total for every processed order in order_intercepted_energy. This \(O(N_\mathrm{iteration})\) diagnostic is kept in normal results because it makes convergence auditable without retaining per-order voxel arrays. Reusing mutable grid-sized workspaces across returned band or timestep results is deferred: the current implementation favors clear result ownership and thread-safe read-only caches. A future in-place workspace API may reduce the remaining warmed allocations without changing equations.
AbstractVoxelTerrain represents a continuous planar, height-field, or triangulated soil receptor. A ray path stops at its nearest terrain intersection. Energy is split into absorbed and reflected fractions, then reflected over the existing full-sphere directions with weights proportional to quadrature solid angle times \(\max(0,n\cdot\omega)\) around the patch's local normal. Set its band reflectances with opaque(par=..., nir=...): reflectance is supplied explicitly, transmission is zero and absorption is 1 - reflectance. Triangle ground uses the same one-sided optical contract. Triangle plant models using translucent retain their equal diffuse reflection/transmission approximation. See Voxel Terrain And Soil for the geometry, AMAPVox/DTM import, boundary, and cache contracts.
Without terrain, lower-boundary radiation is reported as bottom escape. With a black soil it is reported as soil absorption. With a reflective soil it can be intercepted by foliage again, strike another terrain patch, or leave through the top or an open side. Attach the ground geometry with add_ground! and assign its optical properties through models_for(; ground=soil, ...).
The CPU reference solver supports constant projection coefficients, isotropic vegetation scattering, open or horizontally periodic regular grids, scalar or dense effective optical fields, point-proxy individual/organ mixtures, and Lambertian planar, height-field, or triangle terrain. It does not yet model specular reflection, polarized light, fluorescence, wavelength-resolved transport, anisotropic leaf phase functions, non-spherical angular leaf distributions, adaptive voxels, or GPU execution.
The isotropic assumption uses only \(\rho+\tau\) during redistribution. The separate coefficients are retained so the API and input provenance do not erase measured information and a later anisotropic model can distinguish the two hemispheres.
Transport caches are constructed locally and treated as read-only after creation. Independent tasks may safely read one cache concurrently only if their output arrays are distinct; callers must not mutate the cache dictionaries or their stored paths. The series API owns its cache within one call and does not expose shared mutable work buffers.
For related numerical context, see the voxel exchange treatment in RATP and the iterative discrete-direction transport in DART.
Analytical one- and two-voxel cases plus a tiny dense deterministic exchange matrix provide the current independent numerical oracles. A Monte Carlo oracle was not needed for the isotropic centre-source implementation; it remains a useful future challenge test if sub-voxel source sampling or anisotropic phase functions are introduced.