Calculations API Reference

getvar(data, var::Symbol, unit::Symbol=:standard; kwargs...) -> Vector{Float64}
getvar(data, vars::Vector{Symbol}, units::Vector{Symbol}; kwargs...) -> Dict
getvar(data, vars::Vector{Symbol}, unit::Symbol=:standard; kwargs...) -> Dict
getvar(grav_or_rt, gas::HydroDataType, var_or_vars, unit_or_units; kwargs...)
getvar()

Return a stored variable or a derived quantity for every cell, particle or clump of data, in any unit.

Stored columns (e.g. :rho, :vx) are returned as they are; derived quantities (e.g. :mass, :T, :r_sphere, :vϕ_cylinder, :lz) are computed on request. getvar() without arguments prints the list of quantities; the formulas are on the Computation Reference page.

Arguments

• data: any Mera data object, e.g. from gethydro, getparticles, getgravity, getrt or getclumps, or a selection made from one.
• var, vars: one quantity, or several.
• unit, units: the unit of the result, one for all or one per quantity, e.g. :Msol, :km_s, :K, :g_cm3. :standard (the default) means code units. They can also be given as the keywords unit and units.
• gas: for gravity energies and forces, and for the RT rates that need the gas state, the hydro object of the same cells; see the notes below. The two objects can be given in either order.

Keywords

Frame:

• center = [0, 0, 0]: the origin for quantities measured about a point or an axis, such as radii, cylindrical and spherical velocities and angular momenta. [:bc] is the box centre.
• center_unit::Symbol = :standard: the unit of center; :standard means fractions of the box.
• direction::Symbol = :z: the axis of cylindrical quantities; it also comes last in Cartesian output, as for getpositions.
• vcenter = nothing: the velocity of the frame for velocity quantities: a vector [vx, vy, vz] in vunit, :auto for the bulk velocity of data (bulk_velocity), or a function f(x, y, z) of the position (code units, from center) for an ordered flow such as rotation. Without it, velocities are in the frame of the box.
• vunit::Symbol = :standard: the unit of vcenter.

Selection:

• mask: a Bool vector with one entry per row; only rows marked true are returned.
• filtered_db: a table derived from data.data, used instead of it.

Particles:

• ref_time = info.time: the time from which stellar ages are counted, in code units.

Returns

A Vector{Float64} for one quantity, a Dict of vectors, keyed by quantity, for several.

Examples

julia> gas = gethydro(info);

julia> getvar(gas, :mass, :Msol)                         # cell masses in solar masses

julia> d = getvar(gas, [:mass, :T], [:Msol, :K]);        # several quantities

julia> d[:T]

julia> getvar(gas, [:vx, :vy, :vz], :km_s);              # the same unit for all

julia> getvar(gas, :vϕ_cylinder, :km_s; center=[:bc], vcenter=:auto)   # rotation about the centre

julia> getvar(grav, gas, :total_binding_energy, :erg)    # needs gravity and hydro

See also getpositions, getvelocities, msum, subregion, bulk_velocity.

Get mass-array from the dataset (cells/particles/clumps/...):

getmass(dataobject::HydroDataType)
getmass(dataobject::PartDataType)
getmass(dataobject::ClumpDataType)

return Array{Float64,1}

add_field(name::Symbol, compute::Function; depends_on=Symbol[], datatypes=:hydro,
          unit::Symbol=:standard, description::String="")

Register a user-defined derived field that then behaves like any built-in getvar quantity — it works in getvar, and therefore in projection, profile, phase, etc.

• compute(dataobject, deps) — your kernel. deps is a Dict{Symbol,Vector} holding the arrays of depends_on (already centered / masked consistently). Return the field in code units; the requested unit (or this field's default unit) is applied for you.
• depends_on — the variables your kernel needs (built-in or other user fields). These are also recorded in the dependency graph so getvar_requirements (and the read-only-what-you-need logic in project/quicklook) cover your field.
• datatypes — a kind symbol or collection of them: :hydro, :gravity, :rt, :particle, :clump.
• unit — default unit symbol (must be a field of info.scale, or :standard).

add_field(:vmag2, (o, d) -> d[:vx].^2 .+ d[:vy].^2 .+ d[:vz].^2; depends_on=[:vx,:vy,:vz])
getvar(gas, :vmag2)
projection(gas, :vmag2)

See also delete_field, list_fields.

optional lists columns that change the result when they are loaded but are not needed, as described for getvar_optional; variants is a short text, shown in the field listings, that says how the result differs with and without them.

msum(data, unit::Symbol=:standard; mask=[false]) -> Float64

Return the total mass of the cells, particles or clumps in data.

The sum of getvar(data, :mass), so it counts only the part of each cell inside a region: a selection made with a region value, subregion(gas, sphere(10)), stores the fraction of each boundary cell that lies inside, and the mass uses it. Neighbouring regions therefore add up exactly. Selections by the loaders' ranges or the older subregion(gas, :sphere; ...) keep whole cells, and particles and clumps are points.

Arguments

• data: any Mera object with masses, e.g. from gethydro or getparticles, or a selection.
• unit::Symbol = :standard: the unit of the result, e.g. :Msol or :g; :standard means code units. It can also be given as the keyword unit.

Keywords

• mask: a Bool vector with one entry per row; only rows marked true are counted.

Examples

julia> msum(gas, :Msol)                               # all the gas

julia> msum(subregion(gas, sphere(10; center=[:bc])), :Msol)   # inside 10 kpc

julia> msum(gas, :Msol; mask=getvar(gas, :T, :K) .< 1e4)       # the cold gas

See also center_of_mass, subregion, getvar.

center_of_mass(data, unit::Symbol=:standard; mask=[false], periodic=false) -> (x, y, z)
center_of_mass([gas, part], unit::Symbol=:standard; mask=[[false], [false]]) -> (x, y, z)

Return the mass-weighted mean position of the cells, particles or clumps in data, measured from the box corner. com is a short name for it.

Given a vector of objects, e.g. gas and stars, it returns their joint centre of mass.

Arguments

• data: any Mera object with masses, or a vector of hydro and particle objects.
• unit::Symbol = :standard: the unit of the result, e.g. :kpc; :standard means code length. It can also be given as the keyword unit.

Keywords

• mask: a Bool vector with one entry per row; only rows marked true count. For a vector of objects, one mask per object.
• periodic = false: average on a periodic box, with the circular mean (Bai & Breen 2008), so a structure that crosses a face of the box is not averaged to the middle of the box. true for all three axes, or per axis as (x=true, y=true, z=false) or a 3-tuple. getinfo reports whether a run is periodic (info.boundaries). Single objects only.

Examples

julia> center_of_mass(gas, :kpc)

julia> center_of_mass([gas, part], :kpc)              # gas and particles together

julia> center_of_mass(clump, :kpc; periodic=true)     # an object on the box boundary

See also bulk_velocity, msum.

com(data, unit::Symbol=:standard; kwargs...) -> (x, y, z)

A short name for center_of_mass, with the same arguments and keywords (mask, periodic).

bulk_velocity(data, unit::Symbol=:standard; weighting=:mass, mask=[false]) -> (vx, vy, vz)

Return the mean velocity of the cells, particles or clumps in data, weighted by mass by default.

This is the velocity of the frame that getvar(...; vcenter=:auto) subtracts. average_velocity is another name for it.

Arguments

• data: any Mera object with masses and velocities, or a selection.
• unit::Symbol = :standard: the unit of the result, e.g. :km_s; :standard means code units. It can also be given as the keyword unit.

Keywords

• weighting::Symbol = :mass: :mass, :volume (hydro only) or :no for a plain mean. Any other value raises an error.
• mask: a Bool vector with one entry per row; only rows marked true count.

Examples

julia> bulk_velocity(gas, :km_s)

julia> v0 = bulk_velocity(stars, :km_s; mask=getvar(stars, :r_sphere, :kpc; center=[:bc]) .< 10);

julia> getvar(stars, :vϕ_cylinder, :km_s; center=[:bc], vcenter=v0, vunit=:km_s)

See also center_of_mass, getvar.

restframe(dataobject; vcenter, vunit=:standard, center=[0.,0.,0.], center_unit=:standard, mask=[false])

Return a copy of dataobject with a velocity frame subtracted, so every later call sees the boosted velocities. getvar takes vcenter= directly, but projection does not, so this is how a frame reaches a projected quantity such as :σlos.

vcenter accepts what getvar accepts: a 3-vector, :auto for the mass-weighted bulk velocity, or a function f(x, y, z) giving an ordered velocity field, which is the form that removes a rotation curve.

f       = rotation_frame(gas; center=:bc)          # measured from the data
gas_rot = restframe(gas; vcenter=f, center=:bc)
projection(gas_rot, :σlos, :km_s; direction=:edgeon, center=:bc)

A constant boost leaves any dispersion unchanged, because a dispersion already subtracts the mean in each pixel. A varying field does not: it removes an ordered gradient the pixel mean cannot see, which is why an edge-on :σlos drops once the rotation is taken out.

rotation_frame(dataobject; nbins=100, center=[:bc], rmax=nothing) -> Function

Build a vcenter function that subtracts the object's own mean rotation at each cell's radius, leaving the residual motion.

This is the local bulk velocity in the sense that matters for a dispersion: it is measured from the data rather than assumed. Cells are binned by cylindrical radius, the mass-weighted mean azimuthal velocity is taken per bin, and the returned closure evaluates that curve at any position and hands back the corresponding ordered velocity vector.

f       = rotation_frame(gas; center=:bc)
gas_rot = restframe(gas; vcenter=f, center=:bc)   # projection has no vcenter; apply it first
# take the camera from the gas that still rotates, then reuse it: removing the rotation
# drives the net L of `gas_rot` to nearly zero, so asking it for :edgeon points nowhere
a = projection(gas,     :σlos, :km_s; center=:bc, direction=:edgeon)
b = projection(gas_rot, :σlos, :km_s; center=:bc, los=a.los)

Without it, an edge-on :σlos is dominated by ordered rotation along the sightline, because one ray crosses many radii. With it, what is left is the genuine spread about the rotation curve.

nbins sets the radial resolution of the measured curve. Radii beyond the outermost bin reuse the outermost value rather than extrapolating.

average_velocity(data, unit::Symbol=:standard; weighting=:mass, mask=[false]) -> (vx, vy, vz)

Another name for bulk_velocity, with the same arguments and keywords (weighting, mask).

average_mweighted(dataobject, var::Symbol; mask=[false]) -> Float64

Mass-weighted mean of var over the cells or particles in dataobject: $\langle q \rangle_m = \sum m_i q_i / \sum m_i$.

The mass weight is the natural one for intensive quantities — it follows the dense gas, whereas a volume weight follows the diffuse. Use wstat when you also want the median, spread or higher moments, or to weight by something other than mass.

average_mweighted(gas, :T)              # mass-weighted mean temperature, code units
average_mweighted(gas, :T, mask=hot)    # over a subset only

See also wstat, center_of_mass, bulk_velocity.

wstat(values::Vector{<:Real}, weight::Vector{<:Real}=[1.0]; mask=[false]) -> WStatType

Return summary statistics of values, optionally weighted.

With weights, the mean, median and standard deviation are weighted; the skewness and kurtosis are computed about the weighted mean without weights. The standard deviation is the population value (divided by N, not N - 1), and the kurtosis is the excess kurtosis, 0 for a Gaussian.

Arguments

• values: the numbers, e.g. from getvar.
• weight: one weight per value, e.g. the cell masses. It can also be given as the keyword weight; without it all values count equally.

Keywords

• mask: a Bool vector with one entry per value; only values marked true count.

Returns

A WStatType with the fields mean, median, std, skewness, kurtosis, min and max.

Examples

julia> s = wstat(getvar(gas, :T, :K), getvar(gas, :mass));   # mass-weighted temperature

julia> s.mean, s.median, s.std

See also average_mweighted, getvar.

clumping(dataobject; weight=:volume, grid=nothing, grid_unit=:standard,
         mask=[false], verbose=true) -> NamedTuple

The clumping factor $C = \langle n^2 \rangle / \langle n \rangle^2$ of the gas density.

C = 1 is a uniform medium; larger means the mass is concentrated into a smaller fraction of the volume. It is the usual way to quantify how much unresolved structure a simulation has, and the usual way to get burned:

Cell-by-cell (the default) weights each cell by weight:

• :volume — the volume-weighted mean density, the convention C is normally defined with.
• :mass — mass-weighted. This answers a different question — how clumpy is the gas the mass is actually in — and it does not simply give a larger number. For a two-phase medium whose mass sits overwhelmingly in the dense phase the weighted distribution is narrow, so C → 1: equal volumes at n = 1 and 99 give C = 1.96 volume-weighted but 1.01 mass-weighted. Neither is the other's approximation.
• :none — a plain average over cells, which weights a tiny cell the same as a huge one and is almost never what you want on an unstructured mesh.

On a fixed grid (grid = L) the cells are first deposited onto a cube of side L, and C is computed over those equal-volume bins — the form that compares to published values:

clumping(gas)                                    # cell-by-cell, volume-weighted
clumping(gas; grid=13.7, grid_unit=:kpc)         # fixed-grid, comparable to papers
clumping(gas; mask=getvar(gas,:cellsize,:pc) .< 500)   # resolution-restricted

Returns (C, mean_n, mean_n2, n_cells, weight, grid, grid_cells, empty_fraction) with mean_n in cm⁻³.

See also getvar for :cellsize, and the Zoom Simulations page.

max_bins = 5e7 is a guard against a grid that would not fit in memory: a request for more bins raises an error that names the likely cause (often a grid_unit mix-up).

gettime(output::Real, path::String="./", unit::Symbol=:standard) -> Float64
gettime(data, unit::Symbol=:standard) -> Float64

Return the time of a snapshot, in any time unit.

For a run without cosmology this is the simulation time. For a cosmological run (see iscosmological) RAMSES stores a conformal time, so gettime returns the age of the universe at the snapshot's scale factor instead, from cosmology.

Arguments

• output: an output number, read from path (a keyword or the second argument, default "./").
• data: an InfoType from getinfo, or any Mera data object.
• unit::Symbol = :standard: the unit of the result, e.g. :Myr or :Gyr; also a keyword. :standard means code units, and seconds for a cosmological run, which supports only :Gyr, :Myr, :yr and :s.

Examples

julia> gettime(info, :Myr)

julia> gettime(gas, :Gyr)

julia> gettime(300, "path/to/simulation", :Myr)    # without loading anything first

See also getinfo, cosmology.

printtime(text::String="", verbose::Bool=verbose_mode)

Print a Mera timestamp with optional text message to the screen when verbose mode is enabled.

Arguments

• text::String="": Optional text message to display before the timestamp
• verbose::Bool=verbose_mode: Control output display (uses global verbose_mode by default)

Examples

# Print timestamp with default message
printtime()

# Print timestamp with custom message
printtime("Starting calculation")

# Override verbose setting
printtime("Debug info", true)

Returns nothing.

cosmic_time(info::InfoType, a; unit=:Gyr)

The age of the universe at scale factor a, for this run's own cosmology. Accepts a scalar or any array, and is vectorised over it.

gettime answers for the snapshot; this answers for any epoch, which is what tracking across a set of outputs needs — an output tells you its scale factor, and a alone is not a time.

cosmic_time(info, 0.2276)                     # Gyr at the snapshot's scale factor
cosmic_time(info, [gc.aexp for gc in cats])   # …across every catalogue

See also lookback_time, age_of_universe.

Returns a Float64 in unit for a single a, or a vector for a vector of scale factors.

lookback_time(info::InfoType, a; unit=:Gyr, from=1.0)

Time elapsed between scale factor a and from — by default from = 1.0, i.e. the standard lookback time from the present day. Pass from = info.aexp to measure back from the snapshot instead. Scalar or array.

lookback_time(info, 0.2276)                  # Gyr since z = 3.39, from today
lookback_time(info, 0.1,  from=info.aexp)    # …from this snapshot instead

Returns a Float64 in unit for a single a, or a vector for a vector of scale factors.

age_of_universe(info::InfoType; unit=:Gyr)

The age of the universe today (a = 1) for this run's cosmology — cosmic_time(info, 1).

A natural sanity check on a run's cosmological parameters, and one people otherwise compute by hand: a Planck-like ΛCDM gives ≈ 13.8 Gyr, so a value far from it means the parameters read from the header are not what you assumed.

Returns a Float64 in unit.

stellar_age(info::InfoType, birth; unit::Symbol=:Gyr)

Physical age of star particle(s) for a cosmological RAMSES run, from their super-conformal :birth time(s) (scalar or array, as returned by getvar(particles, :birth)). The age is t_proper(a_snap) − t_proper(a_birth) obtained from the Friedmann table (see cosmology); info.time provides the snapshot's conformal time.

unit is a time unit symbol like elsewhere in Mera: :Gyr (default), :Myr, :yr, :s (:standard ⇒ seconds). Non-star sentinels (birth = 0) and any birth time ≥ the snapshot return age 0. This is the conversion used internally by getvar(particles, :age) on cosmological runs.

part = getparticles(info)                 # cosmological run
ages = stellar_age(info, getvar(part, :birth))          # [Gyr]
ages = stellar_age(info, getvar(part, :birth), unit=:Myr)

age_from_aform(info::InfoType, aform; unit::Symbol=:Gyr)

Stellar age(s) for data whose formation time is stored as the scale factor a_form (AREPO / IllustrisTNG GFM_StellarFormationTime, exposed by Mera as :aform), rather than as the RAMSES super-conformal :birth used by stellar_age.

Wind particles — which TNG marks with a_form < 0 — return NaN, as does any non-finite entry.

unit accepts :Gyr (default), :Myr, :yr, :s.

stars = getparticles(info; families=[4])
ages  = age_from_aform(info, getvar(stars, :aform))          # [Gyr]

getunit(dataobject, quantity::Symbol, vars, units; uname=false) -> Real
getunit(info::InfoType, unit::Symbol; uname=false) -> Real

Return the numerical factor that converts quantity from Mera's internal code units into the requested unit, or 1.0 when :standard (code units) is asked for.

This is the conversion every getvar(gas, :rho, :g_cm3)-style call performs internally. Reach for it directly when you hold raw arrays and need the same factor Mera would apply — multiplying by it is exactly what the unit argument does.

getunit(info, :Msol)                       # grams-per-code-mass -> Msol factor
rho_cgs = getvar(gas, :rho) .* getunit(info, :g_cm3)   # equivalent to getvar(gas, :rho, :g_cm3)

With uname=true the unit's name is returned alongside the factor, which is what the plotting helpers use to label axes.

Mixing a code-unit array with a CGS constant is the classic source of silently wrong answers — this function is how you avoid it.

See also createscales, createconstants, getvar.

center_of(data; method=:com, unit=:standard, mask=[false])

Find the centre of an object and return [x, y, z] in unit.

• method=:com — mass-weighted centre of mass (delegates to center_of_mass).
• method=:densest (:peak) — position of the densest hydro cell (needs hydro data).

mask (a Bool/BitArray over the cells/particles) restricts the calculation.

face_on(data; center=:com, aperture=nothing, range_unit=:standard)

Return a GalaxyFrame oriented face-on: the line of sight is the object's angular-momentum (spin) axis, so a projection with los=fr.los sees the disk from above.

• center — :com (default), :densest, or an explicit [x,y,z] in range_unit.
• aperture — optional sphere radius (in range_unit) around the centre; measure the spin only from gas inside it, to isolate the disk from the halo/outskirts.

fr = face_on(gas)                         # whole object, centred on the CoM
fr = face_on(gas; aperture=10, range_unit=:kpc)   # spin from the inner 10 kpc
projection(gas, :sd; los=fr.los, up=fr.up, center=fr.center, range_unit=fr.center_unit)

fr = face_on(gas; center=:densest, aperture=30, range_unit=:kpc)   # the densest galaxy
fr = face_on(gas; center=[x,y,z],  aperture=30, range_unit=:kpc)   # a catalogued halo

edge_on(data; center=:com, aperture=nothing, range_unit=:standard)

Return a GalaxyFrame oriented edge-on: the line of sight lies in the disk plane (perpendicular to the spin axis) and the spin axis points up in the image. Same arguments as face_on.

angmom(data; center=[:bc], range_unit=:standard, radius=nothing, mask=nothing, unit=:standard,
       vcenter=nothing)
angmom([gas, stars, dm]; ...) -> Vector{Float64}

Net angular momentum L = Σ mᵢ rᵢ × vᵢ about center, as a 3-vector.

This is the vector behind direction=:faceon, direction=:edgeon and axis=:angmom. Those presets compute it from one object; calling angmom yourself lets you decide what defines the orientation, then hand the answer to any view:

L = angmom(gas; center=[:bc], radius=[10, :kpc])
projection(gas, :sd; axis=L, inclination=60, center=[:bc])   # same camera for every dataset

Several components at once

A galaxy's plane is set by gas, stars and dark matter together, and each of the three has its own axis. Pass them as a vector and the momenta are added:

L = angmom([gas, stars, dm]; center=[:bc], radius=[20, :kpc])

The objects must share a snapshot; they need not share a resolution or a particle count. The sum is over mass, so a component contributes in proportion to how much of it there is: in most discs the dark matter dominates the total even though the gas defines the visible plane. Decide which tracer you mean, rather than combining out of habit.

Arguments

• center, range_unit: the point L is measured about. Change this and the vector changes, because r is measured from it.

• radius: restrict the sum to a sphere about center, as [10, :kpc] or a bare number in range_unit (or in radius_unit, when the centre is expressed in a different unit). Without it the sum runs over everything the object holds, where the outer material has the longest lever arm and can tilt a "face-on galaxy" by the halo.

• mask: a boolean per row, so the sum can be restricted to, say, cold gas or one particle family (see getparticlemask). It combines with radius: cold gas inside 10 kpc.

• unit: unit for the returned components. The direction does not depend on it, and only the direction is used for orientation.

• vcenter: the velocity frame, as in getvar: :auto for the object's own bulk velocity, or a vector [vx, vy, vz] in the velocity unit you pass (e.g. vunit=:km_s).

Velocity frame

By default velocities are taken as they are stored, in the box frame, and no hint is printed about it. For an object streaming through the box this changes |L|, and tilts it, because r × v picks up r × v_bulk. Measuring about the object's own centre of mass largely cancels it (Σ mᵢ rᵢ = 0 there), which is what face_on does with center=:com. To remove the bulk motion outright, pass vcenter. With several components, give one explicit vector rather than :auto, so that all of them are measured in the same frame.

See also face_on, edge_on, center_of, projection.

sfr_snapshot(p::PartDataType; windows=[5.0, 10.0, 100.0], time_unit=:Myr, mass=:auto,
             mask=[false], eta_sn=:auto, t_sn_delay=5.0) -> NamedTuple

Star-formation rate from a single snapshot, from the star particles (birth ≠ 0; cosmological birth times are converted to ages via stellar_age). Two complementary measures:

• Instantaneous (recent window). For each look-back window Δt in windows, SFR(Δt) = M_*(age ≤ Δt) / Δt — the standard observational "current SFR" (Hα ≈ 5–10 Myr, FUV ≈ 100 Myr). Returned in M⊙/yr.
• Lifetime mean. total stellar mass / age of the oldest star, in M⊙/yr.

mass selects the mass field; :auto (default) prefers a stored initial-mass column and falls back to current :mass — SFR should use the initial stellar mass (current mass is reduced by post-formation mass loss). Returns (; windows, time_unit, sfr, sfr_mean, n_stars, stellar_mass_Msol, oldest_age, mass_field, eta_sn) where sfr is the per-window vector aligned to windows. With no star particles every rate is 0.0.

eta_sn, t_sn_delay work exactly as in sfr: when only the current :mass is stored, the default :auto rebuilds birth masses with the fraction the run recorded, and the returned eta_sn says which value was applied (0.0 when nothing was). The shortest window is usually younger than t_sn_delay, so it is unaffected; the 100 Myr window and the lifetime mean move the most.

s = sfr_snapshot(parts)            # default 5/10/100 Myr windows + mean (auto initial-mass)
s.sfr                              # [SFR(5 Myr), SFR(10 Myr), SFR(100 Myr)]  M⊙/yr
s.sfr_mean                         # lifetime-averaged SFR  M⊙/yr
s.mass_field                       # which mass field was used (e.g. :minit or :mass)
s.eta_sn                           # the SN mass-loss fraction applied (0.0 if none)

See also sfr for the full star-formation history SFR(t).

depletion_time(dataobject, sfr_Msol_yr; mass=:mass, mask=[false]) -> NamedTuple

Gas depletion time and star-formation efficiency per free-fall time. Given a gas region and its star-formation rate sfr_Msol_yr [M⊙/yr], returns

• t_depl_Gyr = M_gas / SFR — how long the present SFR would take to consume the gas;
• t_ff_mw_Myr — the mass-weighted mean free-fall time ⟨√(3π/32Gρ)⟩ (per-cell :freefall_time);
• eps_ff = SFR · ⟨t_ff⟩ / M_gas — the dimensionless star-formation efficiency per free-fall time (Krumholz–McKee), typically ~0.01–0.1;
• M_gas_Msol, sfr.

Use a mask (e.g. dense star-forming gas getvar(gas,:rho,:nH) .> 27) to measure the efficiency of the gas actually forming stars. Works on any grid data with :mass and :freefall_time.

_, s = sfr_snapshot(stars).sfr[2], 0   # or any SFR estimate in M⊙/yr
d = depletion_time(gas, 1.5; mask = getvar(gas,:rho,:nH) .> 27)
d.t_depl_Gyr, d.eps_ff

timeseries(path, reducer; kwargs...)

Run reducer on every snapshot of a simulation and collect the results into a single table — one row per output, ordered by output number.

reducer receives the loaded data object of one snapshot and returns either a scalar or a NamedTuple. The returned table always carries an output column and a time column (physical time in Myr by default, from gettime; a cosmological run additionally gets redshift and aexp columns). A scalar reducer value lands in a value column, a NamedTuple is expanded into one column per field.

Snapshots are loaded strictly one at a time and released before the next, so memory stays bounded — suited to a laptop with limited RAM. Loading respects JULIA_NUM_THREADS (cap it at 4 on a laptop); snapshots are processed sequentially.

Arguments

• path::String : simulation directory (folder holding output_xxxxx/ for RAMSES, or output_xxxxx.jld2 mera files when mera_files=true).
• reducer : function data -> scalar | NamedTuple.

Keywords

• datatype::Symbol = :hydro : :hydro, :gravity, :particles, :clumps, or :rt.
• outputs = :all : :all, a range (1:10), or an explicit vector of output numbers. Numbers not present on disk are skipped.
• mera_files::Bool = false : load mera (JLD2) files via loaddata instead of raw RAMSES outputs.
• loader = nothing : custom info -> data to fully control loading (overrides datatype/ranges/lmax). Use e.g. loader = info -> gethydro(info, [:rho]; lmax=6).
• lmax = nothing : max AMR level to read (hydro/gravity); nothing uses info.levelmax.
• xrange, yrange, zrange, center, range_unit : spatial selection passed to the loader — cutting the region is the main lever to reduce RAM per snapshot.
• time_unit::Symbol = :Myr : unit for the time column — physical by default (:Myr/:Gyr/…); pass :standard for code units (see gettime). A cosmological run also gets redshift and aexp columns automatically.
• verbose::Bool = true : print per-snapshot progress.
• notify::Bool = false : call notifyme when finished (a no-op unless ~/email.txt / ~/zulip.txt is configured).

Examples

# evolution of total gas mass and peak density across all outputs
ts = timeseries("/data/sim/timeseries_sedov3d", d -> (
        mass    = msum(d, :Msol),
        rho_max = maximum(getvar(d, :rho)),
     ))

# same, but from mera files and only every region of interest (less RAM)
ts = timeseries("/data/sim/timeseries_sedov3d_mera",
                d -> msum(d, :Msol);
                mera_files=true, xrange=[0.4,0.6], yrange=[0.4,0.6], zrange=[0.4,0.6])

See also checkoutputs, gettime, gethydro, loaddata.

smallr is a density floor passed to gethydro for each snapshot; 0 switches it off.

pdf(dataobject, quantity; weight=:mass, norm=:density, logbins=true, bins=60,
    valrange=nothing, unit=:standard, mask=[false]) -> NamedTuple

Probability distribution function of a getvar quantity over the cells/particles of dataobject — works on hydro, particle, gravity, and RT data (any quantity/weight getvar supports; for signed fields like the potential :epot use logbins=false). The classic use is the density PDF — the log-normal (with a power-law high-density tail) signature of supersonic turbulence and star formation. To take the PDF of a projected 2D map (the column-density / N-PDF), pass a projection result instead (see below).

Returns (centers, edges, pdf, logbins, norm, quantity, unit, weight).

Keywords

• weight — :mass (default), :volume, or :cells/:count (number-weighted).
• norm — how pdf is normalised:
  • :density (default) — a probability density on the binning axis (log10(quantity) when logbins, so per dex); unit area, sum(pdf .* diff(logbins ? log10.(edges) : edges)) == 1. The bin-width-independent proper PDF.
  • :probability — per-bin probability mass, sum(pdf) == 1.
  • :peak — shape only, scaled so maximum(pdf) == 1.
  • :count (:none) — raw weighted counts, sum(pdf) == total weight.
• logbins — log-spaced bins over log10(quantity) (default; quantity must be > 0).
• bins — number of bins; valrange — (min, max) of the quantity (default: data range).
• unit — unit of quantity; mask — restrict to selected cells/particles.

P  = pdf(gas, :rho)                          # mass-weighted density PDF (area = 1)
Pv = pdf(gas, :rho; weight=:volume)          # volume-weighted (turbulence log-normal)
Pp = pdf(gas, :rho; norm=:probability)       # bins sum to 1
Pk = pdf(gas, :rho; norm=:peak)              # peak = 1 (compare shapes)
# plot: lines(log10.(P.centers), P.pdf)

P = pdf(gas, :rho)
dlog10 = diff(log10.(P.edges))
sum(P.pdf)              # 10.6376  — NOT 1, this is a density
sum(P.pdf .* dlog10)    # 1.000000

pdf(m::DataMapsType, var; weight=:area, norm=:density, logbins=true, bins=60, valrange=nothing)

PDF of the pixel values of a projected 2D map m.maps[var] (a projection result). With var=:sd and the default area weighting this is the column-density PDF (N-PDF) — a standard observational diagnostic. weight is :area (every pixel counts equally), :value (weight by the pixel value), or another map key Symbol in m.maps (weight pixel-by-pixel by that map). Other keywords are as for the cell/particle method.

pdf(map2d::AbstractMatrix; weights=nothing, norm=:density, logbins=true, bins=60, valrange=nothing)

PDF of the values of a raw 2D array (e.g. a mock_observe image or any matrix). weights is an optional matrix of the same size (default: equal per pixel).

getvar_optional(kind::Symbol, vars) -> Vector{Symbol}

The optional columns of vars — ones that change the result when present but are not needed for it to work. Empty for most fields.

The case this exists for is AREPO gas temperature: getvar(gas, :T) needs :u, but takes μ from :ne when that was loaded and otherwise falls back to a neutral-primordial μ ≈ 1.22. The two differ by up to a factor ~2 for ionised gas, so which variant ran depends on the vars= used at load time. getvar_requirements deliberately does not report these, so they never make a valid load look insufficient.

getvar_requirements(:particles, :T)   # [:u]        — what it needs
getvar_optional(:particles, :T)       # [:ne]       — what would improve it