How Quantities Are Computed

This page documents the exact formulas behind Mera's derived quantities and aggregate statistics, for transparency and reproducibility. Every formula here is transcribed directly from the implementation; the code lives in src/functions/getvar/getvar_hydro.jl, getvar_gravity.jl and src/functions/basic_calc.jl. For the machine-readable dependency graph (which raw variables each derived field needs) and the add_field extension API, see Derived Fields & add_field.

Which data to load

Each section is marked with the data object(s) the quantity is defined on — load that type from the same info (getinfo) and call getvar on it:

Multi-type names. Geometry (:x/:y/:z, :r_cylinder, :r_sphere, :ϕ), velocities (:v, :vr_cylinder, …), angular momentum (:hx, :lz, …) and :ekin share one name across hydro, particles and clumps (Julia dispatches on the object). :cellsize/:volume exist for the AMR cell types (hydro/gravity/RT), not particles/clumps. The RT ionization quantities (:xHII, :mu, :T_rt, :n_*, …) are passive hydro scalars — request them on gethydro of an RT run — whereas the photon-group fields live on the getrt object.

Thermodynamics

Data: hydro (gethydro) — needs :rho, :p. (:ekin/:mass are also defined on particles and clumps.)

• Temperature :T $= (p/\rho)\cdot s_K$. The :K unit scale already folds in a constant mean molecular weight: $s_K = \tfrac{m_H}{k_B}\big(\tfrac{\mathrm{unit}_l}{\mathrm{unit}_t}\big)^2\,\mu$ with $\mu = 1/X = 1/0.76 \approx 1.32$ — the neutral-primordial value from RAMSES cooling_module.f90 (hydrogen mass fraction $X=0.76$). So :T is a physical temperature, but it assumes neutral gas everywhere. The raw "temperature per μ" (RAMSES $T/\mu$) is the separate scale scale.T_mu $= \tfrac{m_H}{k_B}(\mathrm{unit}_l/\mathrm{unit}_t)^2$ — i.e. the assumed $\mu$ is exactly scale.K/scale.T_mu. In ionized gas the true $\mu$ falls to $\approx 0.6$, so :T overestimates $T$ by up to $\sim\!2\times$ there; use the ionization-aware :T_rt (see the Radiative-transfer (RT) quantities section below) in RT runs.
• Sound speed :cs is adiabatic, with $\gamma$ taken from info.gamma (typically 5/3).
• Kinetic energy :ekin is bulk motion only, with $v^2 = v_x^2+v_y^2+v_z^2$; it does not include thermal/random energy (that is :etherm).
• Thermal energy :etherm is $p\,V$ over the cell volume $V$ (equivalently $E_\mathrm{therm} = \tfrac{p}{\gamma-1}V$ up to the EOS constant, in the RAMSES convention $p=(\gamma-1)\rho e$).

Entropy family

With $k_B$ Boltzmann's constant, $m_u$ the atomic mass unit and $\gamma$ the adiabatic index:

Velocities & geometry

Data: hydro or particles (velocities); the geometry names :r_cylinder, :r_sphere, :ϕ, :x/:y/:z also work on gravity, RT and clumps.

Positions $x,y,z$ are relative to center (pass center=[:bc] for the box centre). Components that divide by a radius are set to 0 where that radius is zero (on axis / at the centre), rather than returning NaN.

(The acceleration components :ar_cylinder, :aϕ_cylinder, :ar_sphere, :aθ_sphere, :aϕ_sphere and the magnitude :a_magnitude use the identical projections with $\mathbf a$ in place of $\mathbf v$; see the Gravity section below.)

Angular momentum

Data: hydro or particles — needs mass + velocity + position.

Specific angular momentum $\mathbf h = \mathbf r \times \mathbf v$ (per unit mass), and the total $\mathbf L = m\,\mathbf h$:

The cylindrical/spherical angular-momentum components follow the same pattern (mass × the corresponding specific component).

Mach numbers

Data: hydro. The magnetosonic Mach numbers (:mach_alfven, :mach_fast, :mach_slow) need an MHD run. RAMSES stores the field as 6 constrained-transport faces (:b{x,y,z}_left/right); Mera detects MHD from the hydro_file_descriptor (any descriptor version), names the variables canonically (:rho/:vx/:vy/:vz/:p at their true indices), and derives the cell-centred components :bx,:by,:bz = ½(B_left + B_right). The Alfvén speed is vₐ = |B|/√(4πρ) (Gaussian-CGS).

The magnetosonic numbers require an MHD run with :bx,:by,:bz; the magnetic field is taken in RAMSES code units and converted to Gaussian-CGS internally (hence the $4\pi$). They error if the field components are absent.

Magnetic quantities

Data: hydro (MHD). Built from the cell-centred field $\mathbf B = (B_x,B_y,B_z)$. Code-unit convention: RAMSES-MHD absorbs the Gaussian $4\pi$ into the field, so in code units the magnetic pressure is $P_\mathrm{mag} = B^2/2$ and the Alfvén speed is $v_A = |\mathbf B|/\sqrt{\rho}$; the factor reappears only in the physical-unit conversion $B_\mathrm{phys}[\mathrm{G}] = B_\mathrm{code}\cdot\texttt{scale.Gauss}$, $\texttt{scale.Gauss} = \sqrt{4\pi\,\rho_0 v_0^2}$. No new unit type is required.

All five reuse existing unit scales — magnetic-field strengths (:Gauss/:muG/:microG/:nG/:Tesla), pressure (:Ba/:g_cm_s2), velocity (:km_s/:cm_s) and energy (:erg) — so introducing MHD analysis needed essentially no new units. The quantities require an MHD run and error if the field components are absent. (:nG, nanogauss = scale.Gauss·10⁹, is the one unit added — via a versioned ScalesType003 so pre-existing mera files still load.)

Jeans & collapse

Data: hydro — needs :cs (:p) and :rho.

With $G$ the gravitational constant (info.constants.G), $\Delta x$ the cell size and $m$ the cell mass:

Gravity

Data: gravity (getgravity) — needs :epot and/or :ax,:ay,:az.

From the gravitational potential $\phi$ (:epot) and acceleration $\mathbf a$ (:ax,:ay,:az):

The $\max(\cdot,0)$ clamp on the escape speed avoids a negative argument where the potential is unbound ($\phi \ge 0$, possible near domain boundaries) — those cells return 0 rather than erroring.

Radiative-transfer (RT) quantities

Data: hydro of an RT run (gethydro) — the ionization fractions are passive hydro scalars. The photon-group fields (:Np, fluxes, :Gamma_HI, …) live on the getrt object.

These need an RT run: the ionization fractions are passive hydro scalars located via the RT descriptor (info.descriptor.rt, key :iIons), and each quantity errors with a clear message on a non-RT run. RAMSES-RT stores them in a fixed order — $[x_\mathrm{HI}$ (only with H₂ chemistry) $, x_\mathrm{HII}, x_\mathrm{HeII}, x_\mathrm{HeIII}$ (only with He) $]$ — but writes no isH2 flag, so Mera infers the layout from the species count: $n_\mathrm{Ions} = 1 + \mathtt{isH2} + 2\,\mathtt{isHe}$ ⇒ $\mathtt{isH2} = \mathrm{iseven}(n_\mathrm{Ions})$ ($\in\{2,4\}$) and $\mathtt{isHe} = n_\mathrm{Ions}\ge 3$, and remaps every species accordingly. The hydrogen number density used throughout is

\[n_H = \rho\,\cdot\,\texttt{scale.nH}\,\cdot\,\frac{X}{0.76},\]

i.e. scale.nH $= (0.76/m_H)\,\mathrm{unit}_d$ rescaled by the run's actual hydrogen mass fraction $X$ from the descriptor (so a pure-hydrogen $X=1$ Strömgren test is correct; the factor is 1 for the default $X=0.76$).

Mean molecular weight & RT temperature

where the hydrogen particle count per H nucleus is $h_p = 1 + x_\mathrm{HII}$ without H₂ chemistry, or $h_p = x_\mathrm{HI} + 2x_\mathrm{HII} + x_{\mathrm{H_2}}$ with it ($h_p \to 1,\,0.5,\,2$ for neutral-atomic, fully-ionized, fully-molecular pure H). $X_H = X(1-Z)/(X+Y)$ and $X_\mathrm{He} = Y(1-Z)/(X+Y)$ (so $X_H+X_\mathrm{He}+Z = 1$ per cell), with $X,Y$ the primordial fractions, $Z$ the local metal mass fraction (:metallicity; $0$ if absent), and $A_Z \approx 16$. He is taken neutral when not tracked; metal free electrons are neglected (sub-percent). So $\mu$ runs from $\approx 2$ (fully molecular) through $\approx 1.32$ (neutral atomic) to $\approx 0.5\text{–}0.6$ (ionized). On a non-RT run, :mu returns the constant scale.K/scale.T_mu and :T_rt reduces exactly to :T.

Ionization & molecular fractions, number densities

with $n_\mathrm{He} = n_H\,Y/(4X)$ (the He term in $n_e$ enters only when He is tracked; H₂ is neutral and contributes no electrons). :xH2/:n_H2 require an H₂-chemistry run (even $n_\mathrm{Ions}$) and error otherwise.

Recombination

:em_recomb projected with mode=:sum is a mock recombination-line (e.g. Hα) emission map of an HII region. :recomb_rate uses the RT-aware temperature :T_rt (clamped to $\ge 1\,\mathrm{K}$ in the $\alpha_B$ power law) and pairs with the RT photoionization rate for ionization-balance checks.

Cell size & volume

Data: any AMR cell type — hydro, gravity or RT (not particles/clumps).

For an AMR cell at refinement level (uniform-grid runs use lmax), with box length $L_\mathrm{box}$:

\[\Delta x = \frac{L_\mathrm{box}}{2^{\text{level}}}, \qquad V = (\Delta x)^3 .\]

Aggregate statistics

Data: any loaded type (hydro / particles / gravity / clumps), depending on the field requested.

These operate over a whole data object (with optional mask), and live in basic_calc.jl.

Total mass — msum

\[M_\mathrm{total} = \sum_i m_i .\]

Centre of mass — center_of_mass / com

Mass-weighted mean position (returned as a 3-tuple):

\[\mathbf r_\mathrm{cm} = \frac{\sum_i m_i\,\mathbf r_i}{\sum_i m_i} .\]

Bulk velocity — bulk_velocity

Mass-weighted by default; volume-weighted (hydro only) or unweighted on request:

\[\mathbf v_\mathrm{bulk}^{\text{(mass)}} = \frac{\sum_i m_i\,\mathbf v_i}{\sum_i m_i}, \qquad \mathbf v_\mathrm{bulk}^{\text{(vol)}} = \frac{\sum_i V_i\,\mathbf v_i}{\sum_i V_i}, \qquad \mathbf v_\mathrm{bulk}^{\text{(none)}} = \operatorname{mean}(\mathbf v) .\]

Weighted statistics — wstat

wstat returns a WStatType with the weighted mean, median, standard deviation, skewness, kurtosis, and extrema:

\[\bar{x} = \frac{\sum_i w_i x_i}{\sum_i w_i}, \qquad \sigma = \sqrt{\frac{\sum_i w_i (x_i-\bar{x})^2}{\sum_i w_i}} .\]

• The standard deviation is the population form (corrected=false — no Bessel $n/(n-1)$ correction).
• The weighted median uses StatsBase.median(x, Weights(w)); skewness and kurtosis use StatsBase evaluated at the weighted mean.
• Without weights it reduces to the ordinary mean/median/population-std.

Binned reductions — profile, phase, profile3d

Data: any 3-D data — hydro, particles, gravity or clumps.

profile (1-D), phase (2-D) and profile3d (3-D) bin cells/particles by one/two/three axis fields and reduce a target field $y$ in each bin (weighted — mass by default, or :volume, a field, or unweighted). Per bin, with members $i$, weights $w_i$, values $y_i$ and $S_w=\sum_i w_i$:

Also available: var, neff, optional bootstrap mean_ci/median_ci + median_se (nboot>0); bins :linear / :log / :equal (quantile-spaced, equal-population); geometry=:spherical/:cylindrical adds density = S_w/\text{shell volume}; cumulative adds cumsum (e.g. enclosed mass $M(<r)$); normalize=:pdf returns a normalised PDF. Two wrappers build on this:

Conceptual guide and worked examples: Profiles & Phase Diagrams.

Projection maps — projection

Data: hydro (and particles); gravity via the combined hydro+gravity interface.

projection deposits cells onto a 2-D pixel grid (mass-conservatively; see Off-axis Projection). Per pixel, with deposited weight $W=\sum w$ (mass by default) and field $q$:

The dispersion maps are built from two deposited maps, $\langle v\rangle$ and $\langle v^2\rangle$, so $\sigma$ is the spread about that pixel's own weighted-mean velocity — the local line-of-sight (or component) dispersion, with the per-pixel mean (the bulk + rotation seen down that column) removed by construction; the $\max(\cdot,0)$ guards round-off. The axis-aligned $:σ*$ are map-only and need direction=:x/:y/:z; :σlos works for any off-axis line of sight.

Velocity dispersion — which σ am I getting?

Mera never subtracts a single global bulk velocity from a dispersion: every $\sigma$ is a weighted variance about the local mean of the set it is computed over, so net bulk/rotation/streaming at that scale cancels automatically. Only the set differs:

velocitydispersion(…; thermal=true) and localdispersion(…) additionally fold in the thermal line width $\sigma_\mathrm{th}=\sqrt{k_B\langle T\rangle/(\mu m_H)}$ to give the total $\sigma_\mathrm{tot}=\sqrt{\sigma_\mathrm{turb}^2+\sigma_\mathrm{th}^2}$ an observer measures.

So profile(gas, :r_cylinder, :vϕ_cylinder) returns both the mean $\langle v_\phi\rangle(R)$ (the kinematic rotation curve — it keeps its sign) and the std $\sigma_\phi(R)$ (the spread about it). A projected σ (profile a per-pixel :σlos map vs. radius) and a 3-D per-bin σ answer different questions — see the σ note in Profiles & Phase Diagrams.

Worked example: Mach number end-to-end

The derived quantities compose, so you can reproduce any of them by hand. For the Mach number of a cell with $\rho$, $p$ and velocity $(v_x,v_y,v_z)$:

# from first principles
γ  = gas.info.gamma                      # adiabatic index (e.g. 5/3)
cs = sqrt.(γ .* getvar(gas, :p) ./ getvar(gas, :rho))   # = getvar(gas, :cs)
v  = getvar(gas, :v)                                      # √(vx²+vy²+vz²)
M  = v ./ cs                                              # = getvar(gas, :mach)

# the one-liner Mera gives you
M2 = getvar(gas, :mach)
M ≈ M2     # true (same computation)

Every entry above is computed exactly this way internally — getvar simply wires the raw stored variables through these formulas (and the dependency registry in Derived Fields & add_field records which raw variables each one needs).