Calculations API Reference

Get variables or derived quantities from the dataset:

• overview the list of predefined quantities with: getinfo()
• select variable(s) and their unit(s)
• give the spatial center (with units) of the data within the box (relevant e.g. for radius dependency)
• relate the coordinates to a direction (x,y,z)
• pass a modified database
• pass a mask to exclude elements (cells/particles/...) from the calculation

getvar(gas, :r_sphere, :kpc)                  # about the box CORNER (reminder shown)
getvar(gas, :r_sphere, :kpc, center=[:bc])    # about the box centre
getvar(gas, :x, :kpc)                         # absolute coordinates — no reminder

getvar(gravity, hydro, :total_binding_energy, :erg)   # either object order works

getvar(   dataobject::DataSetType, var::Symbol;
        filtered_db::IndexedTables.AbstractIndexedTable=IndexedTables.table([1]),
        center::CenterType=[0.,0.,0.],
        center_unit::Symbol=:standard,
        direction::Symbol=:z,
        unit::Symbol=:standard,
        mask::MaskType=[false],
        ref_time::Real=dataobject.info.time)

return Array{Float64,1}

Arguments

Required:

• dataobject: needs to be of type: "DataSetType"
• var(s): select a variable from the database or a predefined quantity (see field: info, function getvar(), dataobject.data)

Predefined/Optional Keywords:

• filtered_db: pass a filtered or manipulated database together with the corresponding DataSetType object (required argument)
• center: in units given by argument center_unit; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
• center_unit: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
• direction: todo
• unit(s): return the variable in given units
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
• ref_time: reference zero-time for particle age calculation

Defined Methods - function defined for different arguments

• getvar( dataobject::DataSetType, var::Symbol; ...) # one given variable -> returns 1d array
• getvar( dataobject::DataSetType, var::Symbol, unit::Symbol; ...) # one given variable with its unit -> returns 1d array
• getvar( dataobject::DataSetType, vars::Array{Symbol,1}; ...) # several given variables -> array needed -> returns dictionary with 1d arrays
• getvar( dataobject::DataSetType, vars::Array{Symbol,1}, units::Array{Symbol,1}; ...) # several given variables and their corresponding units -> both arrays -> returns dictionary with 1d arrays
• getvar( dataobject::DataSetType, vars::Array{Symbol,1}, unit::Symbol; ...) # several given variables that have the same unit -> array for the variables and a single Symbol for the unit -> returns dictionary with 1d arrays

Examples

# read simulation information
julia> info = getinfo(420)
julia> gas = gethydro(info)

# Example 1: get the mass for each cell of the hydro data (1dim array)
mass1 = getvar(gas, :mass)  # in [code units]
mass = getvar(gas, :mass) * gas.scale.Msol # scale the result from code units to solar masses
mass = getvar(gas, :mass, unit=:Msol) # unit calculation, provided by a keyword argument
mass = getvar(gas, :mass, :Msol) # unit calculation provided by an argument


# Example 2: get the mass and |v| (several variables) for each cell of the hydro data
quantities = getvar(gas, [:mass, :v]) # in [code units]
returns: Dict{Any,Any} with 2 entries:
  :mass => [8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8…
  :v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  2.28274e-7, 2.…

quantities = getvar(gas, [:mass, :v], units=[:Msol, :km_s]) # unit calculation, provided by a keyword argument
quantities = getvar(gas, [:mass, :v], [:Msol, :km_s]) # unit calculation provided by an argument

# Example 3: get several variables in the same units by providing a single argument
quantities = getvar(gas, [:vx, :vy, :vz], :km_s)
...

Radiative transfer (RT) quantities

For an RT run (rt = getrt(info)) the stored variables are, per photon group g, the photon number density :Np<g> and the flux components :Fx<g>, :Fy<g>, :Fz<g> (code units). Derived RT quantities on the RT object:

• :Fmag<g> flux magnitude |F_g| = √(Fx²+Fy²+Fz²) [code units]
• :Np_total photon number density summed over all groups [code units]
• :reducedflux<g> reduced flux f = |Fg| / (c·Npg), dimensionless in 0,1
• :Np<g>_cgs physical photon number density = Npg · unitnp [photons cm⁻³]
• :Fmag<g>_cgs physical flux magnitude = |Fg| · unitpf [photons cm⁻² s⁻¹]
• :photon_energy_density<g> radiation energy density of group g = Npg · unitnp · egyg [erg cm⁻³] (egyg = mean photon energy from info.descriptor.rt[:group_egy])
• :rad_energy_density total radiation energy density summed over all groups [erg cm⁻³]

The per-group photon properties parsed from info_rt (mean energy, energy bounds, photoionization cross-sections, species→group map) are available in info.descriptor.rtPhotonGroups ([g][:egy_eV], [:csn_cm2], [:cse_cm2], plus [:L0_eV], [:L1_eV], [:spec2group]).

The ionization fractions are passive hydro scalars (located via the RT descriptor info.descriptor.rt[:iIons]); request them on the hydro object (gas = gethydro(info)):

• :xHII, :xHeII, :xHeIII ionization fractions (dimensionless)
• :xHI neutral atomic-hydrogen fraction (a stored scalar with H2 chemistry, else the closure 1 − xHII; dimensionless)
• :xH2 molecular-hydrogen fraction = (1 − xHI − xHII)/2 (H2-chemistry runs only; ½ ⇒ fully molecular)
• :n_HII, :n_HI, :n_e, :n_H2 HII / HI / free-electron / H₂ number density cm⁻³
• :em_recomb recombination-emissivity proxy ∝ nₑ·nHII ≈ nHII² cm⁻⁶
• :mu RT-aware mean molecular weight from the ionization (and, with H2 chemistry, molecular) state and metallicity: μ = 1/[XH·hₚ + (XHe/4)(1+xHeII+2xHeIII) + Z/AZ], where hₚ = 1+xHII (no H2) or xHI+2·xHII+xH2 (H2); XH/XHe scaled by the local metal mass fraction Z (the :metallicity scalar, 0 if absent) and AZ≈16. Metal free electrons are neglected (RT does not track metal ionization).

• :T_rt gas temperature [K] using the local μ (= (P/ρ)·scale.Tmu·μ). Plain :T (unit=:K) bakes in a constant μ (= scale.K/scale.Tmu ≈ 1/0.76 ≈ 1.32, the fixed primordial default — independent of the run's X), so it over-estimates the ionized-gas temperature by μconst/μlocal (e.g. ≈1.32/0.5 ≈ 2.6× for fully-ionized pure hydrogen). Prefer :T_rt for RT runs.

rt   = getrt(info)
gas  = gethydro(info)
f    = getvar(rt, :reducedflux1)                 # reduced flux of group 1
nphot = getvar(rt, :Np1_cgs)                     # physical photon density [cm^-3]
xHII = getvar(gas, :xHII)                        # ionization fraction (hydro scalar)
ne   = getvar(gas, :n_e)                         # free-electron density [cm^-3]
T    = getvar(gas, :T_rt)                        # RT-aware temperature [K] (local μ)

:mu and :T_rt also work on non-RT hydro runs: without tracked ionization fractions they fall back to the constant μ Mera's temperature scaling assumes (μ = scale.K/scale.Tmu), so there `:Trtequalsgetvar(:T, :K). The other RT quantities (:xHII,:xHI,:n*,:emrecomb`) require an RT run and error otherwise.

RAMSES-RT field reference (physical quantity → Mera accessor):

Radiation–matter rates (RT object; use the reduced light speed rt_c_frac·c):

• :Gamma_HI<g>, :Gamma_HI HI photoionization rate per group / total [s⁻¹]
• :photoheating_HI<g>, :photoheating_HI HI photoheating rate per HI atom [erg s⁻¹]

Recombination (hydro object): :recomb_rate = αB(T)·nₑ·nHII cm⁻³ s⁻¹.

Combined radiation+gas — request on the RT object with the matching hydro_data (both loaded over the same cells); Mera asserts the cell sets align:

• :photoionizations = ΓHI·nHI [cm⁻³ s⁻¹]
• :ionization_balance = photoionizations − recombinations (≈0 in local equilibrium)

rt  = getrt(info);  gas = gethydro(info)            # same lmax/ranges → aligned cells
Γ   = getvar(rt, :Gamma_HI)                          # photoionization rate [1/s]
bal = getvar(rt, :ionization_balance, hydro_data=gas)  # Γ·n_HI − α_B·nₑ·n_HII  [cm^-3 s^-1]

Most RT quantities are single-object (photon fields on rt, ionization state on gas); only the coupling terms above need both. For other hydro variables you may also pass hydro_data=gas to getvar(rt, …) to fetch them aligned to the RT cells.

Get gravity data with optional hydro data for advanced energy analysis

Arguments:

• dataobject: needs to be of type: "GravDataType"
• var: select a variable from the database or a predefined quantity

Keyword Arguments:

• hydro_data: optional hydro data object for energy calculations that require density/mass
• center: center position (default: [0.,0.,0.])
• direction: direction for cylindrical coordinates (default: :z)
• ...

Examples:

# Basic gravity analysis
grav_data = getvar(grav, :epot)

# Advanced energy analysis with hydro data (keyword syntax)
energy_density = getvar(grav, :gravitational_energy_density, hydro_data=hydro)
binding_energy = getvar(grav, :gravitational_binding_energy, hydro_data=hydro)

# NEW: Simplified positional syntax
jeans_mass = getvar(grav, hydro, :jeansmass, :Msol)
thermal_energy = getvar(grav, hydro, :etherm, :erg)
mixed_analysis = getvar(grav, hydro, [:epot, :T, :jeanslength], [:erg, :K, :pc])

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.

Calculate the total mass of any ContainMassDataSetType:

msum(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])

return Float64

Arguments

Required:

• dataobject: needs to be of type: "ContainMassDataSetType"

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units) :Msol, :Mearth, :Mjupiter, :g, :kg (of typye Symbol) ..etc. ; see for defined mass-scales viewfields(info.scale)
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
• periodic: treat the box as periodic when averaging positions. false (default) keeps the plain mass-weighted mean. true applies it to all three axes; a run that wraps in some directions only takes (x=true, y=true, z=false) or a 3-tuple. Needed whenever the structure touches a face: without it, a clump on the boundary averages to the middle of the box. getinfo reports whether the run is periodic (info.boundaries).

Sub-regions: the sum is boundary-aware

msum sums getvar(obj, :mass), which honours the per-cell :fraction attached by a value-type sub-region (subregion(gas, Sphere(10)), split=true by default): a boundary cell contributes fraction * rho * volume, i.e. only the part of it inside the region. So the result is the mass inside the boundary, and adjacent regions add up exactly. Interior cells carry fraction = 1 and are unaffected.

Cuts made any other way carry no :fraction and therefore count whole boundary cells: the loaders' xrange/yrange/zrange, the classic symbol subregion/shellregion, covering_grid. Particles and clumps are points and have no fraction at all. See subregion.

Calculate the joint center-of-mass of any HydroPartType:

center_of_mass(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])

return Tuple{Float64, Float64, Float64,}

Arguments

Required:

• dataobject: needs to be of type: "Array{HydroPartType,1}""

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
• mask: needs to be of type MaskArrayAbstractType which contains two entries with supertype of Array{Bool,1} or BitArray{1} and the length of the database (rows)

Calculate the center-of-mass of any ContainMassDataSetType:

com(dataobject::ContainMassDataSetType; unit::Symbol=:standard, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}

Arguments

Required:

• dataobject: needs to be of type: "ContainMassDataSetType"

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
• periodic: treat the box as periodic when averaging positions. false (default) keeps the plain mass-weighted mean. true applies it to all three axes; a run that wraps in some directions only takes (x=true, y=true, z=false) or a 3-tuple. Needed whenever the structure touches a face: without it, a clump on the boundary averages to the middle of the box. getinfo reports whether the run is periodic (info.boundaries).

Calculate the joint center-of-mass of any HydroPartType:

com(dataobject::Array{HydroPartType,1}, unit::Symbol; mask::MaskArrayAbstractType=[[false],[false]])

return Tuple{Float64, Float64, Float64,}

Arguments

Required:

• dataobject: needs to be of type: "Array{HydroPartType,1}""

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
• mask: needs to be of type MaskArrayAbstractType which contains two entries with supertype of Array{Bool,1} or BitArray{1} and the length of the database (rows)

Calculate the average velocity (w/o mass-weight) of any ContainMassDataSetType:

bulk_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}

Arguments

Required:

• dataobject: needs to be of type: "ContainMassDataSetType"

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units) :kms, :ms, :cm_s (of typye Symbol) ..etc. ; see for defined velocity-scales viewfields(info.scale)
• weighting: use different weightings: :mass (default), :volume (hydro), :no
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
• periodic: treat the box as periodic when averaging positions. false (default) keeps the plain mass-weighted mean. true applies it to all three axes; a run that wraps in some directions only takes (x=true, y=true, z=false) or a 3-tuple. Needed whenever the structure touches a face: without it, a clump on the boundary averages to the middle of the box. getinfo reports whether the run is periodic (info.boundaries).

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)
sig = projection(gas, :σlos; vcenter=f, center=:bc, direction=:edgeon)

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.

Calculate the average velocity (w/o mass-weight) of any ContainMassDataSetType:

average_velocity(dataobject::ContainMassDataSetType; unit::Symbol=:standard, weighting::Symbol=:mass, mask::MaskType=[false])

return Tuple{Float64, Float64, Float64,}

Arguments

Required:

• dataobject: needs to be of type: "ContainMassDataSetType"

Optional Keywords:

• unit: the unit of the result (can be used w/o keyword): :standard (code units) :kms, :ms, :cm_s (of typye Symbol) ..etc. ; see for defined velocity-scales viewfields(info.scale)
• weighting: use different weightings: :mass (default), :volume (hydro), :no
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)
• periodic: treat the box as periodic when averaging positions. false (default) keeps the plain mass-weighted mean. true applies it to all three axes; a run that wraps in some directions only takes (x=true, y=true, z=false) or a 3-tuple. Needed whenever the structure touches a face: without it, a clump on the boundary averages to the middle of the box. getinfo reports whether the run is periodic (info.boundaries).

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.

Calculate statistical values w/o weighting of any Array:

wstat(array::Array{<:Real,1}; weight::Array{<:Real,1}=[1.], mask::MaskType=[false])

WStatType(mean, median, std, skewness, kurtosis, min, max)

Arguments

Required:

• array: Array needs to be of type: "<:Real"

Optional Keywords:

• weight: Array needs to be of type: "<:Real" (can be used w/o keyword)
• mask: needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the Array

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.

gettime(output::Real; path::String="./", unit::Symbol=:standard)
gettime(dataobject::DataSetType; unit::Symbol=:standard)
gettime(dataobject::InfoType, unit::Symbol=:standard)

Get the physical simulation time in selected units. Returns a Float64.

For a cosmological run (see iscosmological) info.time is conformal time, so gettime instead returns the age of the universe at the snapshot's scale factor (from cosmology), converted to unit. Supported units for cosmological runs: :Gyr, :Myr, :yr, :s (:standard ⇒ seconds). For a non-cosmological run the behaviour is unchanged (info.time scaled to unit).

gettime(1, path="/path/to/sim", unit=:Myr)
gettime(gas, unit=:Gyr)
gettime(info, :Myr)

Arguments Function 1

Required:

• output: give the output-number of the simulation

Predefined/Optional Keywords:

• path: the path to the output folder relative to the current folder or absolute path
• unit: return the variable in given unit

Arguments Function 2

Required:

• dataobject: needs to be of type: "DataSetType"

Predefined/Optional Keywords:

• unit: return the variable in given unit

Arguments Function 3

Required:

• dataobject: needs to be of type: "InfoType"

Predefined/Optional Keywords:

• unit: return the variable in given unit

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)

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.

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

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.

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.

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.

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