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 is AMR-only (hydro/gravity/RT). :volume is defined for those too, and additionally on GADGET/AREPO gas particles, where the reader stores $V = m/\rho$ per particle, so volume weighting works on that family as well. 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. GADGET/AREPO gas particles carry their own thermodynamics from :u, see below.)

• 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:

Gas particle thermodynamics (GADGET/AREPO family)

Data: particles (getparticles) on a GADGET or AREPO gas snapshot, needs the specific internal energy :u. Requesting these without :u raises an ArgumentError.

These codes store gas as particles rather than AMR cells, so the thermodynamics is computed from $u$ instead of $p/\rho$, with $\gamma = 5/3$. The formulas differ from the hydro table above, in particular $\mu$ is not the constant 1.32 used for RAMSES:

$X_H = 0.76$. The electron abundance $n_e$ comes from the :ne column when the snapshot carries one; without it Mera falls back to neutral primordial gas, $\mu = 4/(1+3X_H) \approx 1.22$. If you are writing a methods section for an AREPO analysis, cite this $\mu$, not the RAMSES constant above.

The magnetic quantities in Magnetic quantities apply unchanged to these gas particles when the snapshot carries :bx/:by/:bz (AREPO/TNG MHD): the columns follow the same code convention as RAMSES-MHD, so the formulas carry over verbatim.

The two center arguments

Everything with a radius, an azimuth or an axis in it is measured about an origin, and Mera has two separate center keywords that are easy to confuse.

They are independent arguments, a region placed at one point can perfectly well be asked for quantities measured about another, and they default differently, for historical reasons:

The rule of thumb: pass center explicitly whenever a quantity's name contains a geometry, and pass the same origin you gave the region.

Why the corner default is kept. For absolute positions, :x, :y, :z, the corner is the right origin: it returns the simulation's own coordinates. Changing that default would silently shift every existing script's positions by half a box. The same argument applies to :cuboid, whose ranges are absolute box coordinates.

What happens if you forget. Nothing is refused, because every origin is a well-defined one, you get a plausible number rather than an error. Mera therefore says so, once per session:

• asking for a frame-relative quantity (:r_sphere, :r_cylinder, :ϕ, and the v*/a*/l*/mach_* sphere- and cylinder-frame families) about the corner prints a [Mera] Hint: once per quantity;
• placing a distance-based region (:sphere, :cylinder, either shell) at the corner prints one once per shape: that region is valid, but only the part inside the box is kept, so a corner-placed sphere keeps an octant.

Absolute positions and :cuboid never trigger it. verbose(false) silences the reminders along with every other Mera message.

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; see The two center arguments above for why the default is the corner). 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), which stores the potential $\phi$ (:epot) and the acceleration $\mathbf a$ (:ax, :ay, :az).

info.namelist_content["&POISSON_PARAMS"]["gravity_type"]

What gravity can and cannot answer today

Mera pairs gravity with hydro. Everything below is the gas measured against the field the run produced, and that is the whole of what a snapshot supports without re-solving Poisson.

Available now

These are correct as they stand. You can use the field without knowing which mass made it: the potential and the acceleration already contain every source.

Not available from a snapshot

profile, pdf and phase on a gravity object. They work, but they weight by :mass unless told otherwise, and gravity carries no density, so there is no mass on that object. Weight by the cells instead:

profile(gravity, :r_cylinder, :epot; center=[:bc], weight=:volume)
phase(gravity, :epot, :a_magnitude; weight=:volume)
pdf(gravity, :epot; weight=:volume, logbins=false)   # phi is negative, so no log bins

hydro and particles need none of this: both carry a mass, so the default weight works. To bring a hydro column onto the gravity object for getvar itself, pass hydro_data=gethydro(info) loaded over the same cells.

From gravity alone

These need nothing but the gravity object. Every component measured about an axis or a centre depends on center, exactly as the velocity components do.

Energy and force, which need the cell mass

A potential is per unit mass, so an energy or a force needs the mass of the cell, and the mass lives on the hydro object rather than the gravity one. Pass both:

getvar(gravity, hydro, :Fg, :dyne)              # per cell
projection(gravity, hydro, :Fg, :dyne)          # as a map

Either object order works in both calls, so getvar(hydro, gravity, …) and projection(hydro, gravity, …) do the same thing. The Greek component names have ASCII spellings too: :Fphi_cylinder is the same quantity as :Fϕ_cylinder, exactly as :vphi_cylinder is for velocity.

Called on gravity alone these raise an error naming the fix, rather than guessing a mass. Load both over the identical lmax and ranges: Mera compares the cell indices of the two objects, not just how many there are, so a mass can never be paired with another cell's potential. Two different cuts holding the same number of cells are refused as well.

On a subregion this matters twice over, because the cell fraction that weights boundary cells comes from the hydro object. Cut both with the same region value:

R  = Sphere(10.)
getvar(subregion(gravity, R), subregion(hydro, R), :total_binding_energy, :erg)

See Subregions for what the fraction is applied to.

Each force is the mass times the acceleration component of the same name, taken from that component rather than recomputed, so the two share one definition and one treatment of center.

:gravitational_energy is negative where the cell is bound, following $\phi$. :total_binding_energy is its negative, positive where bound, which is the sign binding energies are usually quoted in.

Projecting a gravity field

Gravity carries no mass and no density, so it cannot weight its own line-of-sight average. That is why there is no single-argument projection(gravity, ...): the hydro object supplies the weight.

projection(hydro, gravity, :epot)                    # column-mass-weighted mean potential
projection(hydro, gravity, :epot; weighting=[:volume])

The result is a weighted mean along each ray, not a column integral, which is the right treatment for an intensive field: a sum would simply scale with the depth of the box. Note weighting takes a vector. Mass and volume weighting answer different questions and give different numbers, so state which one you used.

Velocity dispersion and frames

Data: hydro or particles. The dispersions are map quantities: they are computed per pixel from the spread of velocities along the ray, so they exist in projection, not per cell.

Each is $\sqrt{\langle v^2\rangle - \langle v\rangle^2}$ over the mass in a pixel.

Thermal broadening

:σ_thermal is a per-cell quantity, the 1D thermal width of the mean gas particle:

\[\sigma_\mathrm{thermal} = \sqrt{\frac{k_B T}{\mu m_H}} = \sqrt{P/\rho} = \frac{c_s}{\sqrt{\gamma}}\]

The middle form is what Mera computes, and it needs no μ: $P/\rho$ is $k_B T/(\mu m_H)$ by the ideal gas law, whatever the ionization state. That matters on an RT run, where μ varies per cell and the constant μ behind plain :T is wrong by a large factor. Use :T_rt and :mu there.

It is isotropic, so there is no directional version. Combine it with a directional bulk dispersion by adding the variances, because a line profile is a convolution:

\[\sigma_\mathrm{total}^2 = \sigma_\mathrm{bulk}^2 + \sigma_\mathrm{thermal}^2\]

This is the width for a particle of the mean mass $\mu m_H$. A line is broadened by the mass of the emitting species, so for a species of mass $m_X$ scale it by $\sqrt{\mu m_H/m_X}$: CO is 28 times heavier than hydrogen and so 5.3 times narrower.

Subtracting an ordered flow

A dispersion already subtracts the mean inside each pixel, so a constant boost cannot change it. What it cannot see is an ordered gradient along the ray: an edge-on sightline crosses many radii rotating at different speeds, so :σlos measures the rotation curve rather than turbulence.

restframe returns an object with a velocity frame subtracted, which is how a frame reaches projection (that function has no vcenter keyword; getvar does).

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

vcenter accepts a 3-vector, :auto for the mass-weighted bulk velocity, or a function f(x, y, z) returning an ordered velocity field. rotation_frame builds that function by binning cells in cylindrical radius and taking the mass-weighted mean $v_\varphi$ per bin.

Magnetic field

Data: hydro of an MHD run. RAMSES stores the field on cell faces as :bx_left/:bx_right and so on; Mera averages the two faces to the cell centre for :bx, :by, :bz.

B is a vector, so it decomposes into cylindrical and spherical components the same way velocity and acceleration do. These are measured about center, and warn if none is given.

ASCII spellings (:bphi_cylinder, :btheta_sphere, :bphi_sphere) resolve to the same quantities. A cell at $r=0$ has no defined direction, so those entries come back as 0.

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. :cellsize is AMR-only; :volume is also available on GADGET/AREPO gas particles, by a different route (below).

Every position Mera reports is a cell centre, not a corner. For integer cell indices $(c_x, c_y, c_z)$ at refinement level:

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

The half-cell offset matters whenever you compare a Mera position against something computed from raw indices, or against another tool's convention, a whole-cell error at the finest level is small, but it is systematic.

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 .\]

For GADGET/AREPO gas particles there is no refinement level, so the reader stores a per-particle volume from the density instead:

\[V = m/\rho .\]

It is NaN where $\rho$ is absent or zero, non-gas particle types, and empty cells, so mask those out before using it as a projection weight.

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).