Profiles & Phase Diagrams

profile(dataobject, xvar [, yvar]; weight=:mass, nbins=50, xrange=nothing, scale=:linear,
        edges=nothing, geometry=:none, cumulative=:none, normalize=:none, statistic=nothing,
        xunit=:standard, unit=:standard, center=[:bc], center_unit=:standard,
        quantiles=[0.16, 0.5, 0.84], mask=[false],
        bootstrap=0, ci=:percentile, confidence_level=0.95, bootstrap_seed=20240601) -> NamedTuple

1D profile from 3D data. Bin by xvar and reduce yvar (or the weight itself) per bin. Any getvar fields, for hydro / gravity / RT / particle data; the classic use is a radial profile with xvar a length (:r_cylinder, :r_sphere, :z, …) about a physical center. (Clumps expose only their own catalogue columns — they have no :mass/:rho/:r_sphere getvar field — so radial/mass profiles are not available for clump data.)

• weight — :mass, :volume, :none (equal weights), or any field; should be non-negative (negative summed weights make the weighted mean/std/quantiles ill-defined → NaN). Mass/volume weighting works for any data type that has that field (:volume is grid-only; particles use :mass).
• center/center_unit — reference point and its unit (e.g. center=[24,24,24], center_unit=:kpc). range_unit is accepted as a back-compat alias of center_unit.
• xrange=(lo,hi) (end point) in xunit (the binning axis can be any quantity, so its unit is xunit, not center_unit/range_unit — unlike the spatial xrange of projection). xunit/unit set the x-axis / yvar field units.
• scale — :linear (default), :log (log-spaced bins), or :equal (quantile-spaced adaptive bins so every bin holds ~the same number of points — robust for sparse outer radii / noisy data).
• bootstrap=N (>0) — also return bootstrap confidence intervals for the per-bin mean and median by resampling each bin N times (default off). ci selects the interval: :percentile (default), :basic, or :bca (bias-corrected & accelerated). confidence_level (default 0.95) and a fixed bootstrap_seed (deterministic). Adds mean_ci/median_ci (nbins×2: lower, upper) + median_se. Cost is O(bootstrap·N_bin) per bin (:bca adds an O(N_bin²) jackknife) — meant for modest counts.
• quantiles — percentile levels for the per-bin weighted quantiles (default 16/50/84%).
• edges — an explicit vector of bin edges (overrides nbins/xrange/scale); given in xunit.
• geometry — :spherical (shell volume 4/3·π·Δr³) or :cylindrical (annulus area π·Δr²) adds shell_volume and a density = sum / shell_volume (e.g. a radial ρ(r) or surface density), in weight-unit per xunit³ (xunit²). Pair with xvar=:r_sphere/:r_cylinder.
• cumulative — :forward (low→high) or :reverse adds cumsum/cumcount (e.g. enclosed mass M(<r)).
• normalize — :sum adds fraction = sum/Σsum (bins sum to 1); :pdf also adds pdf = fraction/Δedge (integral = 1).
• statistic — a function applied per bin (called as f(yview, wview) if it accepts weights, else f(yview)) returning a scalar; result in custom. Needs yvar.

yvar may be a vector of fields ([:T, :rho]): the data is binned once and each field is reduced in the same pass — the per-field statistics are returned under fields (e.g. p.fields[:T].mean) with the order in yvars.

Returns x (centres), edges, count, sum (Σ weight), sumw2 (Σ weight²); with a single yvar also weighted mean, std, sem (standard error on the mean via Kish neff), min, max, median, a quantiles matrix (nbins × length(qlevels)), and the weighted shape moments skewness and kurtosis (excess kurtosis — 0 for a Gaussian). (bootstrap adds the CI fields.)

profile(m::DataMapsType, var; xvar=:r, weight=:none, nbins=50, xrange=nothing, scale=:linear,
        edges=nothing, geometry=:none, cumulative=:none, statistic=nothing,
        xunit=:standard, center=nothing, center_unit=:standard, quantiles=[0.16,0.5,0.84])
    -> NamedTuple

1D profile from a projected 2D map (a projection result). Bin the pixels of m.maps[var] by xvar and reduce per bin with the same statistics as the data method. xvar is

• :r — image-plane radius from center (a surface-brightness / Σ(R) profile),
• :x / :y — an image-plane coordinate, or
• another map key — bin one map against another (any map vs any map).

weight is :none/:area (equal pixels) or another map key (e.g. weight=:sd for a column-weighted profile of :vlos). For the :r/:x/:y cases center (default: the map centre) is in center_unit (range_unit accepted as an alias) while xrange is in xunit (the radius is converted from code units via the map's scale). edges, geometry (:cylindrical annulus area for a 2-D map), cumulative and statistic work as in the data method.

Returns

The same per-bin statistic fields as the data profile method, plus var (the binned map variable), xvar, weight, xunit, unit (the mapped variable's unit, from the projection), and source=:map.

phase(dataobject, xvar, yvar [, cvar]; weight=:mass, nbins=(100,100), xrange=nothing,
      yrange=nothing, xscale=:linear, yscale=:linear, xunit=:standard, yunit=:standard,
      cunit=:standard, weight_unit=:standard, center=[:bc], center_unit=:standard,
      mask=[false]) -> NamedTuple

2D phase diagram — the summed weight in bins of (xvar, yvar) (e.g. a mass-weighted density–temperature diagram, phase(gas, :rho, :T; xscale=:log, yscale=:log)). With a third field cvar, also returns the per-bin weighted mean of cvar. H[i,j] is the total weight in x-bin i, y-bin j. normalize=:sum/:pdf adds fraction/pdf (2-D). With cvar, cstat selects a richer per-bin reduction of cvar in addition to mean: :std, :median, :min, :max, :full (all four), or a function f(cview, wview) (→ custom). Any non-:mean cstat also returns a per-bin weighted quantiles array (nbx × nby × length(quantiles)) at the quantiles levels (with qlevels).

ph = phase(gas, :rho, :T; xscale=:log, yscale=:log, weight_unit=:Msol)
ph.weight, ph.weight_unit          # (:mass, :Msol)  -> label the colour bar from these

Returns

A NamedTuple with xedges, yedges (bin edges), H (the nbx × nby summed-weight grid), xvar/yvar, weight, weight_unit, xunit/yunit, and source=:data. With normalize, also fraction/pdf. With a cvar: mean (and, for a non-:mean cstat, the corresponding std/median/min/max/quantiles/custom grids), plus cvar/cunit.

phase(m::DataMapsType, xvar, yvar [, cvar]; weight=:none, nbins=(100,100), …) -> NamedTuple

2D phase diagram from a projected map — bin two map variables against each other, weighted by :none/:area or another map key, optionally colouring by the weighted mean of a third map cvar.

Returns

The same fields as the data phase method (xedges, yedges, H, the cstat grids when a cvar is given, and fraction/pdf when normalized), with source=:map. weight_unit is :count for :none/:area (H counts pixels) and :unknown for a map key, whose unit the map does not record.

profile3d(dataobject, xvar, yvar, zvar [, cvar]; weight=:mass, nbins=(50,50,50),
          xrange=nothing, yrange=nothing, zrange=nothing, xscale=:linear, yscale=:linear,
          zscale=:linear, xunit=:standard, yunit=:standard, zunit=:standard, cunit=:standard,
          center=[:bc], center_unit=:standard, normalize=:none, mask=[false]) -> NamedTuple

3D profile — the summed weight in bins of three fields (xvar, yvar, zvar), the 3-D generalization of phase. H[i,j,k] is the total weight in the (i,j,k) cell, and sum(H) is the total weight (e.g. a ρ–T–Z mass cube). With a fourth field cvar, also returns the per-bin weighted mean of cvar. normalize=:sum adds fraction = H/ΣH; :pdf also adds pdf (÷ the 3-D cell volume, so the integral is 1). Each axis takes its own range/scale/unit, and nbins is an integer or a 3-tuple. Marginalizing one axis (sum(H; dims=3)) recovers the corresponding phase.

Returns xedges, yedges, zedges, H (and mean with cvar; fraction/pdf if normalized), plus weight/weight_unit. As in phase, H is the summed weight in weight_unit, which defaults to :standard = code units — pass e.g. weight_unit=:Msol to sum in solar masses instead of relabelling afterwards.

rotationcurve(dataobject; center=[:bc], center_unit=:standard, rvar=:r_sphere, nbins=50,
              xrange=nothing, scale=:linear, xunit=:kpc, mask=[false]) -> NamedTuple

Circular-velocity (rotation) curve from the enclosed mass. Bin all cells/particles by radius rvar (:r_sphere or :r_cylinder) about a physical center, form the enclosed mass M(<r) = Σ mass(< r), and return the Newtonian circular velocity v_circ = √(G·M(<r)/r) and the gravitational acceleration g = G·M(<r)/r². Needs a :mass field (hydro / particles); mass-bearing components can be combined by concatenating their radius/mass (see the component-split example in the profiles tutorial).

Returns x (the outer bin-edge radius in xunit, where the enclosed mass is complete), edges, count, m_enc (enclosed mass, M⊙), v_circ (km/s) and g (cm/s²). The cumulative mass Σ mass(< r) is the mass interior to each bin's upper edge, so the velocity/acceleration are evaluated at that same radius (edges[2:end]) — pairing them with the bin centre would mix a half-bin of mass against a smaller radius and overestimate the inner curve.

velocitydispersion(dataobject; rvar=:r_cylinder, components=(:vr_cylinder,:vϕ_cylinder,:vz),
                   weight=:mass, vunit=:km_s, nbins=50, xrange=nothing, scale=:linear, binsize=nothing,
                   center=[:bc], center_unit=:standard, xunit=:kpc, mask=[false],
                   thermal=false, mu=1.0, Tvar=:T) -> NamedTuple

Radial velocity-dispersion profile. Bins by rvar and returns the per-bin weighted standard deviation of each velocity component — the rest-frame dispersion (about the per-bin mean, so net rotation/streaming does NOT inflate it) — plus the total σ = √(σ₁²+σ₂²+σ₃²) and each component's mean. The default cylindrical triplet gives σR / σφ / σz (use `(:vrsphere,:vθsphere,:vϕsphere)for the spherical decomposition). Each σ is exactly thestdof [profile`](@ref) on that component — this is a convenience wrapper; see the profiles tutorial for the manual recipe.

Returns x (bin centres), edges, count, sigma (total), sigma_components (nbins × n) and mean_components (nbins × n) with the components order, and neff (Kish — small ⇒ noisy σ).

Orbital anisotropy. Pass the spherical triad — exactly (:vr_sphere, :vθ_sphere, :vϕ_sphere), with Greek letters, though the ASCII spellings :vtheta_sphere / :vphi_sphere are accepted — and the result also carries sigma_r, sigma_t = √(σ_θ² + σ_φ²) and the anisotropy parameter beta:

vd = velocitydispersion(part; components=(:vr_sphere, :vθ_sphere, :vϕ_sphere),
                        rvar=:r_sphere, center=cen)
vd.beta            # 1 − σ_t²/(2σ_r²):  0 isotropic, →1 radial, <0 tangential

sigma alone is direction-averaged and hides this. The fields appear only for the spherical triad: β is defined against the spherical radius, and computing it from the cylindrical R would give a different quantity that still looks like a plausible β.

Thermal & total dispersion (thermal=true). The kinematic σ above is turbulent (bulk) motion only. Set thermal=true to also fold in the gas thermal motion and report the quantity an observer measures as a line width. It adds, per bin:

• sigma_turb_1d = √(Σσ_i²/n) — the 1-D turbulent dispersion (the sigma field is the 3-D √(Σσ_i²));
• sigma_thermal = √(k_B⟨T⟩/(μ m_H)) — the 1-D thermal speed at mean molecular weight mu, from the mass-weighted ⟨T⟩. mu=:auto (default) follows getvar(:mu): the per-cell value from the ionization state on an RT run, mass-weighted per bin, and the constant Mera's temperature scaling assumes on a run without RT. Give a number when the thermal speed should be that of an observed tracer rather than of the mean gas (2.33 molecular, 1.0 atomic H, 0.6 ionized);
• sigma_total = √(sigma_turb_1d² + sigma_thermal²) — the total (turbulent ⊕ thermal) 1-D dispersion;
• mach = sigma_turb_1d / ⟨c_s⟩ — the turbulent Mach number (⟨c_s⟩ the mass-weighted sound speed);
• cs, T, mu — the per-bin mass-weighted sound speed / temperature and the mu used.

Tvar selects the temperature field. The default :T_rt is the temperature computed with the local μ, which is the right one in ionized gas and is equal to :T to machine precision on a run without RT, so it is safe either way; pass Tvar=:T for the constant-μ temperature. For a local, patch-de-streamed turbulent ⊕ thermal dispersion (TIGRESS/SILCC-style, removing bulk flow on a chosen length scale rather than per-radial-bin), see localdispersion.

localdispersion(dataobject; patchsize=[500,:pc], components=(:vr_cylinder,:vϕ_cylinder,:vz),
                weight=:mass, vunit=:km_s, thermal=true, mu=:auto, Tvar=:T_rt,
                center=[:bc], range_unit=:standard, mask=[false], min_cells_per_patch=20,
                quantiles=[0.16,0.5,0.84]) -> NamedTuple

Local (patch-de-streamed) velocity dispersion — turbulent ⊕ thermal. The TIGRESS/SILCC-style turbulence measure: the field of view is tiled into square patchsize patches in the (x, y) plane, the per-patch mass-weighted mean velocity is subtracted from every cell (so bulk rotation, shear and spiral streaming on scales larger than patchsize are removed), and the residual dispersion is the genuine small-scale turbulence below patchsize.

This differs from velocitydispersion, which removes the mean per radial bin: here the de-streaming is on a spatial grid, so it does not absorb non-circular streaming into "turbulence". Restrict the input first (e.g. shellregion to a solar annulus, and/or a midplane mask) — this returns aggregate scalars over the whole supplied region, not a radial profile.

Per the supplied region it returns the mass-weighted aggregate:

• sigma_r,sigma_phi,sigma_z (the components order), sigma_turb_3d = √(Σσ_i²), sigma_turb_1d = √(Σσ_i²/n);
• with thermal=true: sigma_thermal = √(k_B⟨T⟩/(μ m_H)), sigma_total = √(sigma_turb_1d²+sigma_thermal²);
• mach = sigma_turb_1d/⟨c_s⟩, anisotropy = sigma_z/σ_in-plane (for a 3-component triplet, z last);
• n_cell, n_eff (Kish), n_patch, and the patch-to-patch weighted percentiles of σtotal / Mach / anisotropy at quantiles (`sigmatotalq,machq,anisotropyq) — the physical region-to-region scatter, suitable as error bars. Only patches with ≥mincellsperpatch` cells enter the percentiles.

mu is the tracer mean molecular weight in H-atom masses (2.33 molecular, 1.0 atomic, 0.6 ionized).

profiletimeseries(loadfn, outputs, xvar [, yvar]; field=nothing, time_unit=:Myr, kwargs...)
    -> NamedTuple

Stack a profile across snapshots into a (nbins × n_snapshots) matrix — a radius-vs-time map (e.g. the evolution of a rotation/density/temperature profile). loadfn(output) returns a loaded data object for that snapshot, e.g. out -> gethydro(getinfo(out, path), verbose=false). For each output, profile(loadfn(output), xvar, yvar; kwargs...) is computed and the field field (default :mean with a yvar, else :sum) is stacked as a column; the snapshot time comes from gettime in time_unit.

For a vector yvar ([:rho, :T]) the per-field statistics live under pr.fields[field], so pass field=(:fieldname, :stat) (e.g. field=(:rho, :mean)) to pick which field/statistic to stack; the default is the first field's :mean.

Pass a fixed radius axis (xrange+nbins or explicit edges) so the columns align — otherwise the per-snapshot bin counts differ and an error is raised. Extra kwargs are forwarded to profile.

Returns x (bin centres), edges, t (times in time_unit), outputs, M (nbins × n_snap matrix of field) and field.

getparticlemask(dataobject::PartDataType, select; verbose=true) -> Vector{Bool}

Boolean mask selecting a subset of particles, for use as the mask= argument of profile, phase, rotationcurve (and projection, getvar, …). select may be

• a named type — :all, :dm (:dark_matter), :stars (:star), :clouds (:sink), :debris, :other, :tracer (family ≤ 0), :gas (gas tracer, family 0);
• a family code Int, or a vector of codes (matched against the RAMSES :family column);
• a NamedTuple combining family and/or tag, e.g. (family=2,), (tag=3,), (family=2, tag=1).

On the new RAMSES particle format the :family/:tag columns are used (DM=1, star=2, cloud=3, debris=4, other=5, tracers ≤ 0). On the legacy format (no :family column) only :stars (birth ≠ 0) and :dm (birth == 0) are available via the :birth field; other selections raise an error. The required column must be among the loaded particle variables.

parts = getparticles(getinfo(1, "spiral_ugrid"))
profile(parts, :r_cylinder; weight=:mass, mask=getparticlemask(parts, :stars),
        center=[:bc], range_unit=:kpc, xunit=:kpc)
rotationcurve(parts; mask=getparticlemask(parts, :dm), center=[:bc], range_unit=:kpc)