Radiation

closure_error(F) -> Vector

1 - Σⱼ Fᵢⱼ per element: the fraction of each element's emission that leaves the geometry entirely (zero for every element of a closed enclosure, up to discretization error).

view_factors(target; near = 3.0, level = 1) -> Matrix{Float64}

Dense view-factor matrix F[i, j]: the fraction of diffuse radiation leaving element i that arrives at element j. Reciprocity Aᵢ Fᵢⱼ = Aⱼ Fⱼᵢ holds exactly by construction, and rows sum to one for a closed enclosure (counting apertures); see closure_error.

Flat elements use the exact differential-area-to-polygon view factor onto the other element, integrated over the emitter: over 4^level sub-triangle points for pairs closer than near element sizes (√2 × the largest centroid-to-vertex distance, the side of a square element), which stays accurate for elements sharing an edge such as the wall–cap corner of a hohlraum, and from the centroid beyond that, averaged over the two directions.

Elements of a component with a CurvedSurface are integrated over the curved patch: quadrature points on the true surface with its normals and areas, point to point between two curved elements (4^level sub-triangle points per triangle when near, one per triangle out to three times that distance, the projected centroid beyond), and onto the exact polygon of a planar partner. Elements of a DiskSurface are exact polar sectors, planar receivers whose arcs are divided finely. A curved concave element also sees itself, which gives the diagonal Fᵢᵢ. Flat facets standing in for a curved wall misrepresent the exchange between nearby elements by an amount of order one relative to its size, which matters next to a concentrated source such as a laser spot.

Elements of more than two facets (remeshed CAD, remesh_mesh) are not integrated facet by facet: a flat one receives through its outline (exactly the sum over its facets, to round-off; piece by piece, each piece its facets in one plane, where the point's tangent plane cuts it), and near pairs integrate over cells that lump the emitter's facets (side size / (2^(level+2) √2) before refining toward the receiver, 2^(level+1) on a curved surface), so a pair costs a few times a two-facet pair whatever the facet count — as long as the receiver is: each point's factor sums its outline, or its pieces when the point's tangent plane cuts it, so an element on a finely faceted curved wall (hundreds of sliver facets, tens of pieces and a jagged outline of hundreds of edges) costs tens of times more per near pair. Give such parts their fitted surface (curved_components): they are then integrated on the curved surface instead.

Visibility is a single centroid-to-centroid ray between the flat facets; the facets of a curved surface (cylinder, sphere or cone) do not block a ray between two points inside it, since the true surface cannot. On curved elements the centroid facing test is relaxed by the spread of the true normals across the flat facet, so pairs the curved patches exchange across are kept (a general triangle on a cylinder has its plane tilted axially by about its size over 2R), and the point kernels apply the exact facing per point.

Enclosure

Result of enclosure: the reduced target, and how it was cut from source.

• target: the components of source cut to their kept patches (components left empty are dropped).
• index[i]: the element of source that element i of target is.
• kept[e]: whether element e of source passed the selection. A two-sided patch keeps both sides when either passes, so target can hold sides with kept false.
• flux[e]: the relative flux on element e of source that decided it (NaN for elements beyond radius, which cast no rays); see enclosure.
• escape: the fraction of the rays from the point that leave the bounding sphere (or reach infinity): the solid-angle fraction of the openings seen from the point.
• open_area: Σ over kept elements of area × fraction of the element's rays that escape (mm²). By reciprocity this estimates the open area of the cavity (its entrance holes, where no aperture or window closes them).
• leak_area: the same sum over rays that end on dropped elements (mm²): radiation the kept set sends to surfaces the reduced target no longer has.
• return_area: what those surfaces would send back, Σ over kept elements of area × the mean over the element's rays of albedo × flux at the dropped element hit (mm², in units of the reference flux): the error of the cut, to compare with open_area (what leaves through the openings) and the cavity's area.
• back_area: the sum over rays that end on one-sided surfaces seen from behind (mm²), which absorb without radiating in the full target too.
• back: the fraction of the rays from the point whose first hit was a one-sided surface seen from behind (1 for a point inside the capsule; otherwise a sign of wrongly oriented surfaces).
• rays: the rays cast per element.

cad_target(parts; orient = true, toward = (0, 0, 0), inside = nothing,
           remesh = nothing, coarsen = nothing, fit = false, material = gold(),
           two_sided = false, enclosure_kw = (;)) -> Target

A radiation-ready target from CAD meshes in one call: what the Radiation in an imported CAD hohlraum tutorial does by hand, and what visrad_target does with the CAD a VisRad workspace bundles. parts holds name => mesh pairs, the CAD parts (as read_stl_solids returns them), and Components, taken as they are (primitives, apertures, parts prepared by hand): those take part in the cavity selection but are not oriented, remeshed or fitted. The target's components keep the order of parts.

Each CAD part becomes a component of material radiating from its front (from both sides when two_sided), through these steps:

1. orient = true repairs its windings with orient_mesh: closed solids face out of their material, open sheets face the point toward (mm).
2. With inside, a point (mm) in the cavity, enclosure keeps only the elements that bound the free space around it, among the facets of all the parts as oriented (enclosure_kw are its keywords). The selection needs consistent windings.
3. remesh (mm) divides the kept facets into compact, near-planar elements of about that size (remesh_mesh); coarsen (mm) instead merges whole facets into elements (coarsen_mesh). Give at most one. Only the kept facets are divided, which is much cheaper when most of an assembly lies outside the cavity. With inside, the selection is then repeated on the elements: an element straddling the edge of the lit region averages its flux over its own area.
4. fit = true (or a named tuple of keywords of curved_components and fit_surfaces, such as minfraction = 1e-3 to leave small curved faces flat) splits each part into its cylinders, cones, spheres and the rest, fitted on the kept CAD facets, whose vertices lie on the true surfaces, so the radiation model integrates over the curved walls.

material, two_sided, remesh and coarsen take one value for every CAD part, a function of the part's name (a Symbol), or a dictionary keyed by name (Symbol or String). A part missing from a remesh or coarsen dictionary keeps its facets, one missing from a two_sided dictionary is one-sided, and one missing from a material dictionary is an error. For example, remesh = n -> n === :capsule ? 0.15 : 0.5 remeshes a capsule finer than the walls. For the account of what the cavity selection drops, call enclosure on the target built without inside.

enclosure(target, point = (0, 0, 0); threshold = 0.05, albedo = 0.9, rays = 64,
          points = 8, point_rays = 4096, radius = Inf, seed = 1) -> Enclosure
enclosure(components, point; kwargs...)

The elements of target that bound the free-space region containing point (mm), as a new Target (enc.target) with each component cut to those patches, with diagnostics (Enclosure; show prints elements and area kept per component). Imported CAD assemblies model a target as closed solids, with outer faces, end faces and mounting structure that cannot exchange radiation with the cavity, while view factors cost O(n²) in the elements. Only element sides that face the free space are found (fronts of one-sided surfaces, either side of two-sided ones), so CAD windings must be consistent first (orient_mesh, visrad_target(...; orient = true)).

The selection is the flux an isotropic source at the point would put on each element, directly and after diffuse re-emission from the surfaces with the uniform albedo albedo (a geometric weight, not the materials'). Every element within radius of the point casts rays rays from random points on it in cosine-weighted directions, stratified in angle, against the target's facets (apertures are transparent, as everywhere), and the flux is gathered along them:

Eⱼ = E⁰ⱼ + albedo × (mean of E over the element sides j's rays hit first),

with rays that escape (leave the sphere of radius about the point, or hit nothing) or meet a one-sided surface from behind contributing nothing. E⁰ⱼ is the direct flux, cos θ/(4π d²) averaged over points random points on the element, zero where the line of sight to the point is blocked; lines of sight (and the rays from the point) pass through surfaces seen from behind, so the point may lie inside the capsule. The iteration (Gauss–Seidel in element order, to 10⁻⁶) starts from E⁰. Elements with E ≥ threshold × E_ref are kept, where E_ref is the mean flux on the element sides that point_rays rays from the point (a Fibonacci lattice) first hit: the cavity's own scale. enc.flux holds E/E_ref. On a NIF CAD workspace (155k elements) this takes about 3 s on 8 threads.

Each element gathers through its own rays, which sample its whole view whatever its size or shape, so a small patch or a sliver deep in a recess or behind a shield is found through the surfaces around it, which see the cavity. Surfaces outside, such as outer walls and the far sides of mounts, get no flux at all. A surface that sees the cavity only through an opening, such as mounting structure beyond an entrance hole, gets about albedo × its view factor to the cavity, plus what its neighbors outside gather in turn, attenuated by albedo at each step, so structure outside does not carry itself in by seeing itself. The default threshold 0.05 drops a surface facing a hole of radius a from more than about 4a away on its axis; the rim of a hole in a thick wall spends about half its view on the cavity and is kept. threshold = 0 keeps every element with any flux at all (whatever is connected to the cavity by lines of sight, however weakly), which radius bounds. The random numbers come from seed and the element index, so the result does not depend on the thread count; an element whose flux is within the sampling noise of the threshold can fall either way.

Removing elements the cavity cannot see changes nothing in the radiosity of the rest; removing a surface seen through an opening turns the radiation it would have returned into an escape, which return_area estimates against open_area (lower threshold when it matters). Thin films that CAD models as solids, such as a capsule tent or an entrance-hole window, bound the free space like walls: with the point inside a tent the result is the tent's pocket around the capsule (open_area near zero). Leave such parts out of the target (they transmit x-rays at their thickness) or give a point on the side wanted.

Stefan–Boltzmann constant in TW mm⁻² eV⁻⁴ (σT⁴ with T in eV gives TW/mm²).

RadiationModel(target; F = view_factors(target), closed = false, kwargs...)

Target plus view factors: everything needed to solve for x-ray fluxes given laser deposition. kwargs pass through to view_factors. Set closed = true when the elements (apertures included) form a closed enclosure — e.g. a hohlraum with its LEHs — to remove discretization closure error with enforce_closure!; radiation that would otherwise leak through the mesh is amplified by multiple bounces at high albedo.

The model caches LU factorizations of the radiosity matrix I - diag(α) F keyed on the albedo vector (the two most recently used), so solve and flux_response with unchanged albedos skip the factorization. The cache is thread-safe. Do not mutate F after the first solve (or call Briareus.empty_cache!(model) if you do).

RadiationSolution

Fields: laser (absorbed laser power per element, TW), source (x-ray source flux ηP/A), B (radiosity), H (incident flux), all per element in TW/mm²; albedo per element (the albedos the radiosity was solved with, so B = source + albedo .* H exactly); intact, the albedo each element would have intact (a thin-foil element's before albedo, the albedo elsewhere); transmission, the fraction of its incident flux each element transmits out of the enclosure (nonzero only on a burned-through BurnthroughAlbedo foil; see energy_balance); time.

albedos(target, time = 0.0, H = nothing; burnthrough = nothing) -> Vector

Albedo per element at time; temperature-dependent models use each element's radiation temperature from incident flux H (or their reference temperature when H is nothing). burnthrough (ns, per element; NaN where not given) sets the burnthrough time of BurnthroughAlbedo elements.

energy_balance(sol) -> NamedTuple

Global x-ray power accounting (TW): source (Σ ηPᵢ), absorbed by each component ((1-α)HA), escaped through each aperture, transmitted through each burned-through thin-foil component (BurnthroughAlbedo) and each open or closing gap (GapClosure, (1-s)HA while a fraction s is closed), and lost to open space (what leaves the geometry without reaching any element). A thin-foil component's loss (1-α)HA splits into what it transmits, τHA with τ its sol.transmission, and what it absorbs, (1-α-τ)HA: with the default transmission T, τ = sT and α = (1-sT)α_before (s the transition's progress), so it absorbs (1-sT)(1-α_before)HA; with an explicit after, τ = α_before-α and it absorbs (1-α_before)HA, as if intact. source ≈ Σabsorbed + Σescaped + Σtransmitted + lost.

enforce_closure!(F, areas; tol = 1e-12, maxiter = 1000) -> F

Rescale view factors so every row sums to one while keeping reciprocity Aᵢ Fᵢⱼ = Aⱼ Fⱼᵢ: the exchange matrix Gᵢⱼ = Aᵢ Fᵢⱼ becomes sᵢ Gᵢⱼ sⱼ with s from symmetric (Sinkhorn) balancing. Only valid for closed enclosures. Rows that see nothing are left untouched.

flux_response(model, observed; time = 0.0, albedo = nothing, weights = nothing) -> Matrix

Linear operator R with H[observed] = R * P: incident x-ray flux on the observed elements (TW/mm²) per unit absorbed laser power on every element (TW). It is fixed for a given geometry and albedo, so optimizers that only change the laser source apply it as a matrix product instead of re-solving. Pass the converged sol.albedo of a reference solution to linearize about a temperature-dependent albedo state. Shares the model's factorization cache with solve.

With weights (a matrix with one column per observed element, e.g. legendre_weights), returns weights * R directly, solving for its rows only: far cheaper than R when only a few moments of the observed flux are wanted.

radiation_temperature(H)

Radiation temperature (eV) of incident flux H (TW/mm²): (H/σ)^(1/4).

solve(model, dep::Deposition; time = 0.0, albedo = nothing, burnthrough = nothing,
      tol = 1e-4, maxiter = 50, guess = nothing)
solve(model, laser_power::AbstractVector; ...)

Solve the radiosity system for the given absorbed laser power per element. Albedos come from the materials at time unless an explicit per-element albedo vector is passed. Temperature-dependent albedos (WallLossAlbedo) are iterated to self-consistency with each element's incident-flux temperature (to tol in albedo). burnthrough (ns, per element, NaN where not given) sets the burnthrough times of BurnthroughAlbedo elements; without it a thin-foil element uses its model's explicit burnthrough, or is intact (drive_history tracks it from the flux history instead).

The factorization of the system for the initial albedo (the materials' at time, or the explicit albedo) is cached on the model, so repeated solves with unchanged albedos cost only a pair of triangular solves and a matrix–vector product. Iterates of a temperature-dependent albedo differ from call to call and are factored afresh.

guess, a previous solution of the same model (e.g. the last step of a history), starts the albedo iteration from its incident-flux temperatures instead of the materials' reference temperatures: fewer iterations when it is close. The result then differs from a cold start by less than the iteration's tolerance. The first albedo is not cached in that case.

Convert a flux in TW/mm² to W/cm².

PulseShape(times, powers)

Piecewise-linear power history (ns, TW per beam), zero outside the table. Callable: pulse(t).

drive_history(model, op, beams_at, times; observe = :capsule, lmax = 6,
              wall = (:wall,), lpi = nothing, holes = nothing, ylm = false,
              transition = nothing, flux_change = nothing, min_step = 1e-3,
              tol = 1e-4, keep_albedo = false) -> NamedTuple

Quasi-static drive history: at each time, deposit beams_at(t) with op (after the laser–plasma model lpi, evaluated at that time), solve the radiosity at that time (time- and temperature-dependent albedos included), and record

• time — the times solved (ns): times, plus any steps inserted by transition or flux_change
• delivered — laser power (TW)
• flux, Tr — mean incident flux on observe (TW/mm²) and its radiation temperature (eV)
• moments — (lmax+1) × nt Legendre moments of the observed flux (aₗ/a₀ rows are moments[l+1, :] ./ moments[1, :])
• ylm — with ylm = true, the (lmax+1)² × nt real spherical-harmonic coefficients of the observed flux (ylm_weights, in mode_index order; columns work with mode_amplitude), otherwise nothing
• wall_Tr, wall_albedo — area-weighted mean incident-flux temperature and albedo of the wall components
• escaped — power (TW) through each aperture, and transmitted through each thin-foil (BurnthroughAlbedo) and closing-gap (GapClosure) component, as in energy_balance (by default the fraction s T of the flux on each burned-through element, and 1 - s on a gap a fraction s closed)
• window_loss — power lost through diagnostic windows: transmitted plus escaped through every aperture except the entrance holes (default: the apertures not on a CylinderSurface, i.e. a hohlraum's LEHs)
• front — per thin-foil component, its area-mean heat-front areal mass as a fraction of the foil's (≥ 1 once burned through), from its own flux history even when the model has an explicit burnthrough
• burnthrough — per thin-foil component, each element's burnthrough time (ns; Inf if it has not burned through)
• albedo — with keep_albedo = true, each step's converged albedo per element (the state a flux response linearizes about, e.g. for retune_drive), otherwise nothing.

Each step's temperature-dependent albedos are iterated to tol (see solve), starting from the previous step's temperatures.

Each thin-foil element without an explicit burnthrough time tracks its own heat front: the drive integral of heat_front over that element's incident-flux temperature history (trapezoidal between the steps, from times[1]), so its burnthrough is consistent with the flux the windowed hohlraum puts on it. A step is solved with the front extended from the previous step at constant temperature, and solved again until the step's own temperatures no longer move a burnthrough that affects its albedo.

A burnthrough transition takes the foil's width (0.4 ns by default), and a step much longer than that jumps across it. Two optional criteria insert steps between times:

• transition — no thin-foil element's transition (the smooth step s of BurnthroughAlbedo, 0 intact to 1 burned through), and no closing gap's (GapClosure, 0 open to 1 closed), advances by more than this in one step. Steps are planned from each element's known or predicted burnthrough time, and halved if the step's own temperatures move a burnthrough so that the advance exceeds 1.5 transition.
• flux_change — the observed flux relative to its mean, H/⟨H⟩ on each element of observe, changes by no more than this in one step; a step that exceeds it is halved and solved again.

Neither shortens a step below min_step (ns). Without them (the default) the steps are exactly times.

Laser energy (kJ = TW·ns) of one beam following the pulse.

with_pulses(beams, pulses, t; by = :group) -> Vector{Beam}

Beams with power set to pulses[key](t), where key is each beam's group (pulses given as a Dict or NamedTuple keyed by group), its quad (by = :quad), its name (by = :beam) or by(beam) for a function, or a single PulseShape for all.

cone_fraction_pulses(beams, times, total, fraction; inner = (:inner23, :inner30))
    -> Dict{Symbol,PulseShape}

Per-group pulses (power per beam, TW) from a total laser power history total (TW, at times ns) and a cone fraction history fraction (inner-cone share of the total), the way NIF pulses are specified: the inner groups' beams share fraction × total equally, the other groups' beams the rest. Keyed by each beam's group, for with_pulses.

read_pulse_table(path) -> Dict{Symbol,PulseShape}

Pulses per group, quad or beam from a CSV table in the layout of write_pulse_table: a time_ns column and one column per key (power per beam, TW). Read a reference shot's measured or requested powers this way.

retune_drive(model, op, beams, reference; by = :group, limits = nothing,
             modes = (1 => 1.0, 2 => 1.0, 3 => 1.0, 4 => 1.0), components = (),
             drive_weight = 1.0, smooth = 0.0, nominal = nothing, stay = 0.0,
             observe = :capsule) -> NamedTuple

Power histories for beams (grouped by :group, :quad, :beam or a function beam -> key; one power per beam of a group at each time) that reproduce a reference drive: reference is a drive_history of the reference shot run with ylm = true and keep_albedo = true, whose capsule-flux coefficients c_lm(t) and albedos are the target. op is the deposition operator of beams (e.g. the reference's restricted to the beams that remain, op[keep]), whose pointing holds through the pulse.

At each reference time the capsule flux is linear in the groups' powers, with the albedos fixed at the reference's converged ones (flux_response(...; albedo)). The powers minimize

Σ_(l => w) w Σ_m ((c_lm - ŝ_lm c_00)/r_00)² + Σ_((l,m) => w) w ((c_lm - ŝ_lm c_00)/r_00)²
    + drive_weight ((c_00 - r_00)/r_00)² + smooth Σ_g ((P_g - P_g,prev)/P_scale)²
    + stay Σ_g ((P_g - P_g,nominal)/P_scale)²,

where ŝ = r/r_00 is the reference shape, within the bounds limits (TW per beam): nothing (from 0, unbounded above), a Dict key => (lo, hi) (missing keys: (0, Inf)), or a function (key, t) -> (lo, hi), e.g. a quad's power capability, or (0, 0) for a dropped quad. smooth ties each step's powers to the previous step's, and stay to nominal powers (a Dict key => PulseShape, or a function (key, t) -> P, e.g. the reference shot's): among the many power sets that restore the drive it prefers the one closest to them, instead of one that runs many quads to their limits. This is a box-bounded linear least-squares problem per time; the drive term keeps the level, so drive_weight must be positive.

Returns (; time, keys, power, pulses, shape, reference_shape, drive, residual): the per-beam powers (length(keys) × nt), the same as PulseShapes for with_pulses (by the same key), the predicted and reference c_lm/c_00 ((lmax+1)² × nt), the predicted c_00/r_00 and the weighted residual norm per time. The prediction is linear about the reference's albedos: confirm with a drive_history of the pulses, where the albedos follow the new drive.

write_pulse_table(path, pulses) -> path

Pulses (power per beam, TW, against time, ns) per group, quad or beam as CSV: a time_ns column and one column per key, on the union of the pulses' times. Read back with read_pulse_table.