4. Basic Calculations
The following functions process different types of DataSetType:
- ContainMassDataSetType,
- HydroPartType,
- HydroDataType,
- PartDataType,
- ClumpDataType.
Load The Data
using Mera
info = getinfo(400, "/Volumes/FASTStorage/Simulations/Mera-Tests/manu_sim_sf_L14");
gas = gethydro(info, [:rho, :vx, :vy, :vz], lmax=8);
particles = getparticles(info, [:mass, :vx, :vy, :vz])
clumps = getclumps(info);[Mera]: 2025-06-30T09:30:13.374
Code: RAMSES
output [400] summary:
mtime: 2018-09-05T09:51:55
ctime: 2025-06-29T20:06:45.267
=======================================================
simulation time: 594.98 [Myr]
boxlen: 48.0 [kpc]
ncpu: 2048
ndim: 3
-------------------------------------------------------
amr: true
level(s): 6 - 14 --> cellsize(s): 750.0 [pc] - 2.93 [pc]
-------------------------------------------------------
hydro: true
hydro-variables: 7 --> (:rho, :vx, :vy, :vz, :p, :var6, :var7)
hydro-descriptor: (:density, :velocity_x, :velocity_y, :velocity_z, :thermal_pressure, :passive_scalar_1, :passive_scalar_2)
γ: 1.6667
-------------------------------------------------------
gravity: true
gravity-variables: (:epot, :ax, :ay, :az)
-------------------------------------------------------
particles: true
- Npart: 5.091500e+05
- Nstars: 5.066030e+05
- Ndm: 2.547000e+03
particle-variables: 5 --> (:vx, :vy, :vz, :mass, :birth)
-------------------------------------------------------
rt: false
-------------------------------------------------------
clumps: true
clump-variables: (:index, :lev, :parent, :ncell, :peak_x, :peak_y, :peak_z, Symbol("rho-"), Symbol("rho+"), :rho_av, :mass_cl, :relevance)
-------------------------------------------------------
namelist-file: false
timer-file: false
compilation-file: true
makefile: true
patchfile: true
=======================================================
[Mera]: Get hydro data: 2025-06-30T09:30:13.766
Key vars=(:level, :cx, :cy, :cz)
Using var(s)=(1, 2, 3, 4) = (:rho, :vx, :vy, :vz)
domain:
xmin::xmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
ymin::ymax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
zmin::zmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
Reading data...
Progress: 100%|█████████████████████████████████████████| Time: 0:02:47
Memory used for data table :51.839996337890625 MB
-------------------------------------------------------
[Mera]: Get particle data: 2025-06-30T09:33:01.147
Key vars=(:level, :x, :y, :z, :id)
Using var(s)=(1, 2, 3, 4) = (:vx, :vy, :vz, :mass)
domain:
xmin::xmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
ymin::ymax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
zmin::zmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
Progress: 100%|█████████████████████████████████████████| Time: 0:00:24
Found 5.089390e+05 particles
Memory used for data table :31.064148902893066 MB
-------------------------------------------------------
[Mera]: Get clump data: 2025-06-30T09:33:26.262
domain:
xmin::xmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
ymin::ymax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
zmin::zmax: 0.0 :: 1.0 ==> 0.0 [kpc] :: 48.0 [kpc]
Read 12 colums:
[:index, :lev, :parent, :ncell, :peak_x, :peak_y, :peak_z, Symbol("rho-"), Symbol("rho+"), :rho_av, :mass_cl, :relevance]
Memory used for data table :61.58203125 KB
-------------------------------------------------------Note
Many functions can provide the results in selected units. The internal scaling from code to physical units makes use of the following defined relations:
viewfields(info.scale)[Mera]: Fields to scale from user/code units to selected units
=======================================================================
Mpc = 0.0010000000000006482
kpc = 1.0000000000006481
pc = 1000.0000000006482
mpc = 1.0000000000006482e6
ly = 3261.5637769461323
Au = 2.0626480623310105e23
km = 3.0856775812820004e16
m = 3.085677581282e19
cm = 3.085677581282e21
mm = 3.085677581282e22
μm = 3.085677581282e25
Mpc3 = 1.0000000000019446e-9
kpc3 = 1.0000000000019444
pc3 = 1.0000000000019448e9
mpc3 = 1.0000000000019446e18
ly3 = 3.469585750743794e10
Au3 = 8.775571306099254e69
km3 = 2.9379989454983075e49
m3 = 2.9379989454983063e58
cm3 = 2.9379989454983065e64
mm3 = 2.937998945498306e67
μm3 = 2.937998945498306e76
Msol_pc3 = 0.9997234790001649
Msun_pc3 = 0.9997234790001649
g_cm3 = 6.76838218451376e-23
Msol_pc2 = 999.7234790008131
Msun_pc2 = 999.7234790008131
g_cm2 = 0.20885045168302602
Gyr = 0.014910986463557083
Myr = 14.910986463557084
yr = 1.4910986463557083e7
s = 4.70554946422349e14
ms = 4.70554946422349e17
Msol = 9.99723479002109e8
Msun = 9.99723479002109e8
Mearth = 3.329677459032007e14
Mjupiter = 1.0476363431814971e12
g = 1.9885499720830952e42
km_s = 65.57528732282063
m_s = 65575.28732282063
cm_s = 6.557528732282063e6
nH = 30.987773856809987
erg = 8.551000140274429e55
g_cms2 = 2.9104844143584656e-9
T_mu = 517028.3199143136
K_mu = 517028.3199143136
T = 680300.4209398864
K = 680300.4209398864
Ba = 2.910484414358466e-9
g_cm_s2 = 2.910484414358466e-9
p_kB = 2.1080995598777838e7
K_cm3 = 2.1080995598777838e7Total Mass
The function msum calculates the total mass of the data that is assigned to the provided object. For the hydro-data, the mass is derived from the density and cell-size (level) of all elements. info.scale.Msol (or e.g.: gas.info.scale.Msol) scales the result from code units to solar masses:
println( "Gas Mtot: ", msum(gas) * info.scale.Msol, " Msol" )
println( "Particles Mtot: ", msum(particles) * info.scale.Msol, " Msol" )
println( "Clumps Mtot: ", msum(clumps) * info.scale.Msol, " Msol" )Gas Mtot: 2.6703951073850353e10 Msol
Particles Mtot: 5.804426008528444e9 Msol
Clumps Mtot: 1.3743280681841675e10 MsolThe units for the results can be calculated by the function itself by providing an unit-argument:
println( "Gas Mtot: ", msum(gas, :Msol) , " Msol" )
println( "Particles Mtot: ", msum(particles, :Msol) , " Msol" )
println( "Clumps Mtot: ", msum(clumps, :Msol) , " Msol" )Gas Mtot: 2.670395107385036e10 Msol
Particles Mtot: 5.804426008528436e9 Msol
Clumps Mtot: 1.3743280681841673e10 MsolThe following methods are defined on the function msum:
methods(msum)2 methods for generic function <b>msum</b> from Mera:<ul><li> msum(dataobject::<b>ContainMassDataSetType</b>; <i>unit, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L22" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:22</a></li> <li> msum(dataobject::<b>ContainMassDataSetType</b>, unit::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L18" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:18</a></li> </ul>
Center-Of-Mass
The function center_of_mass or com calculates the center-of-mass of the data that is assigned to the provided object.
println( "Gas COM: ", center_of_mass(gas) .* info.scale.kpc, " kpc" )
println( "Particles COM: ", center_of_mass(particles) .* info.scale.kpc, " kpc" )
println( "Clumps COM: ", center_of_mass(clumps) .* info.scale.kpc, " kpc" );Gas COM: (23.32748735447764, 23.83541991952592, 24.04172014803584) kpc
Particles COM: (22.891354761211332, 24.174147282680273, 24.003205056545575) kpc
Clumps COM: (23.135765457064576, 23.741712325649267, 24.005012718586205) kpcThe units for the results can be calculated by the function itself by providing a unit-argument:
println( "Gas COM: ", center_of_mass(gas, :kpc) , " kpc" )
println( "Particles COM: ", center_of_mass(particles, :kpc) , " kpc" )
println( "Clumps COM: ", center_of_mass(clumps, :kpc) , " kpc" );Gas COM: (23.32748735447764, 23.83541991952592, 24.04172014803584) kpc
Particles COM: (22.891354761211332, 24.174147282680273, 24.003205056545575) kpc
Clumps COM: (23.135765457064576, 23.741712325649267, 24.005012718586205) kpcA shorter name for the function center_of_mass is defined as com :
println( "Gas COM: ", com(gas, :kpc) , " kpc" )
println( "Particles COM: ", com(particles, :kpc) , " kpc" )
println( "Clumps COM: ", com(clumps, :kpc) , " kpc" );Gas COM: (23.32748735447764, 23.83541991952592, 24.04172014803584) kpc
Particles COM: (22.891354761211332, 24.174147282680273, 24.003205056545575) kpc
Clumps COM: (23.135765457064576, 23.741712325649267, 24.005012718586205) kpcThe result of the coordinates (x, y, z) can be assigned e.g. to a tuple or to three single variables:
# return coordinates in a tuple
com_gas = com(gas, :kpc)
println( "Tuple: ", com_gas, " kpc" )
# return coordinates into variables
x_pos, y_pos, z_pos = com(gas, :kpc); #create variables
println("Single vars: ", x_pos, " ", y_pos, " ", z_pos, " kpc")Tuple: (23.32748735447764, 23.83541991952592, 24.04172014803584) kpc
Single vars: 23.32748735447764 23.83541991952592 24.04172014803584 kpcCalculate the joint centre-of-mass from the hydro and particle data. Provide the hydro and particle data with an array (independent order):
println( "Joint COM (Gas + Particles): ", center_of_mass([gas,particles], :kpc) , " kpc" )
println( "Joint COM (Particles + Gas): ", center_of_mass([particles,gas], :kpc) , " kpc" )Joint COM (Gas + Particles): (23.24961513830681, 23.895900266222746, 24.03484321295537) kpc
Joint COM (Particles + Gas): (23.24961513830681, 23.895900266222746, 24.03484321295537) kpcUse the shorter name com that is defined as the function center_of_mass :
println( "Joint COM (Gas + Particles): ", com([gas,particles], :kpc) , " kpc" )
println( "Joint COM (Particles + Gas): ", com([particles,gas], :kpc) , " kpc" )Joint COM (Gas + Particles): (23.24961513830681, 23.895900266222746, 24.03484321295537) kpc
Joint COM (Particles + Gas): (23.24961513830681, 23.895900266222746, 24.03484321295537) kpcmethods(center_of_mass)4 methods for generic function <b>centerofmass</b> from Mera:<ul><li> centerofmass(dataobject::<b>Vector{HydroPartType}</b>; <i>unit, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L106" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:106</a></li> <li> centerofmass(dataobject::<b>Vector{HydroPartType}</b>, unit::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L102" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:102</a></li> <li> centerofmass(dataobject::<b>ContainMassDataSetType</b>; <i>unit, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L50" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:50</a></li> <li> centerofmass(dataobject::<b>ContainMassDataSetType</b>, unit::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L46" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:46</a></li> </ul>
methods(com)4 methods for generic function <b>com</b> from Mera:<ul><li> com(dataobject::<b>Vector{HydroPartType}</b>; <i>unit, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L165" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:165</a></li> <li> com(dataobject::<b>Vector{HydroPartType}</b>, unit::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L161" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:161</a></li> <li> com(dataobject::<b>ContainMassDataSetType</b>; <i>unit, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L79" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:79</a></li> <li> com(dataobject::<b>ContainMassDataSetType</b>, unit::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L75" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:75</a></li> </ul>
Bulk Velocity
The function bulk_velocity or average_velocity calculates the average velocity (w/o mass-weight) of the data that is assigned to the provided object. It can also be used for the clump data if it has velocity components: vx, vy, vz. The default is with mass-weighting:
println( "Gas: ", bulk_velocity(gas, :km_s) , " km/s" )
println( "Particles: ", bulk_velocity(particles, :km_s) , " km/s" )Gas: (-1.4418303105424646, -11.708719305767852, -0.539324349686299) km/s
Particles: (-11.623422700314535, -18.440572802490234, -0.3291927731417527) km/sprintln( "Gas: ", average_velocity(gas, :km_s) , " km/s" )
println( "Particles: ", average_velocity(particles, :km_s) , " km/s" )Gas: (-1.4418303105424646, -11.708719305767852, -0.539324349686299) km/s
Particles: (-11.623422700314535, -18.440572802490234, -0.3291927731417527) km/sWithout mass-weighting:
- gas: volume or :no weighting
- particles: no weighting
println( "Gas: ", bulk_velocity(gas, :km_s, weighting=:volume) , " km/s" )
println( "Particles: ", bulk_velocity(particles, :km_s, weighting=:no) , " km/s" )Gas: (1.5248458901822848, -8.770913864354457, -0.5037635305158428) km/s
Particles: (-11.594477384589647, -18.388591187193725, -0.30977462952679713) km/sprintln( "Gas: ", average_velocity(gas, :km_s, weighting=:volume) , " km/s" )
println( "Particles: ", average_velocity(particles, :km_s, weighting=:no) , " km/s" )Gas: (1.5248458901822848, -8.770913864354457, -0.5037635305158428) km/s
Particles: (-11.594477384589647, -18.388591187193725, -0.30977462952679713) km/smethods(bulk_velocity)2 methods for generic function <b>bulkvelocity</b> from Mera:<ul><li> bulkvelocity(dataobject::<b>ContainMassDataSetType</b>; <i>unit, weighting, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L199" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:199</a></li> <li> bulkvelocity(dataobject::<b>ContainMassDataSetType</b>, unit::<b>Symbol</b>; <i>weighting, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L194" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basic_calc.jl:194</a></li> </ul>
methods(average_velocity)2 methods for generic function <b>averagevelocity</b> from Mera:<ul><li> averagevelocity(dataobject::<b>ContainMassDataSetType</b>; <i>unit, weighting, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L240" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basiccalc.jl:240</a></li> <li> averagevelocity(dataobject::<b>ContainMassDataSetType</b>, unit::<b>Symbol</b>; <i>weighting, mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L236" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basic_calc.jl:236</a></li> </ul>
Mass Weighted Average
The functions center_of_mass and bulk_velocity use the function average_mweighted (average_mass-weighted) in the backend which can be feeded with any kind of variable that is pre-defined for the getvar() function or exists in the datatable. See the defined method and at getvar() below:
methods( average_mweighted )1 method for generic function <b>averagemweighted</b> from Mera:<ul><li> averagemweighted(dataobject::<b>ContainMassDataSetType</b>, var::<b>Symbol</b>; <i>mask</i>) in Mera at <a href="https://github.com/ManuelBehrendt/Mera.jl/tree/292a22d2938a49234887580dc7b36bab9e28350d//src/functions/basiccalc.jl#L171" target="blank">/Users/mabe/Documents/codes/github/Mera.jl/src/functions/basic_calc.jl:171</a></li> </ul>
<a id="Statistics"></a>
Get Predefined Quantities
Here, we only show the examples with the hydro-data:
info = getinfo(1, "/Volumes/FASTStorage/Simulations/Mera-Tests//manu_stable_2019", verbose=false);
gas = gethydro(info, [:rho, :vx, :vy, :vz], verbose=false); Reading data...
Progress: 100%|█████████████████████████████████████████| Time: 0:00:20Use getvar to extract variables or derive predefined quantities from the database, dependent on the data type. See the possible variables:
getvar()Predefined vars that can be calculated for each cell/particle:
----------------------------------------------------------------
=============================[gas]:=============================
-all the non derived hydro vars-
:cpu, :level, :rho, :cx, :cy, :cz, :vx, :vy, :vz, :p, var6,...
-derived hydro vars-
:x, :y, :z
:mass, :cellsize, :volume, :freefall_time
:cs, :mach, :machx, :machy, :machz, :jeanslength, :jeansnumber
:T, :Temp, :Temperature with p/rho
:h, :hx, :hy, :hz (specific angular momentum)
==========================[particles]:==========================
-all the non derived particle vars-
:cpu, :level, :id, :family, :tag
:x, :y, :z, :vx, :vy, :vz, :mass, :birth, :metal....
-derived particle vars-
:age
===========================[gravity]:===========================
-all the non derived gravity vars-
:cpu, :level, cx, cy, cz, :epot, :ax, :ay, :az
-derived gravity vars-
:x, :y, :z
:cellsize, :volume
===========================[clumps]:===========================
:peak_x or :x, :peak_y or :y, :peak_z or :z
:v, :ekin,...
=====================[gas or particles]:=======================
:v, :ekin
related to a given center:
---------------------------
:r_cylinder, :r_sphere (radial components)
:ϕ
:vr_cylinder
:vϕ_cylinder
----------------------------------------------------------------Get a Single Quantity
In the following example, we calculate the mass for each cell of the hydro data.
- The output is a 1dim array in code units by default (mass1).
- Each element/cell can be scaled to Msol units by the elementwise multiplikation gas.scale.Msol (mass2).
- The
getvarfunction supports intrinsic scaling to a selected unit (mass3). - The selected unit does not need a keyword argument if the following order is maintained: dataobject, variable, unit
mass1 = getvar(gas, :mass) # [code units]
mass2 = getvar(gas, :mass) * gas.scale.Msol # scale the result (1dim array) from code units to solar masses
mass3 = getvar(gas, :mass, unit=:Msol) # unit calculation, provided by a keyword argument [Msol]
mass4 = getvar(gas, :mass, :Msol) # unit calculation provided by an argument [Msol]
# construct a three dimensional array to compare the three created arrays column wise:
mass_overview = [mass1 mass2 mass3 mass4] 37898393×4 Matrix{Float64}:
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
8.9407e-7 894.07 894.07 894.07
⋮
1.02889e-7 102.889 102.889 102.889
1.02889e-7 102.889 102.889 102.889
1.94423e-7 194.423 194.423 194.423
1.94423e-7 194.423 194.423 194.423
8.90454e-8 89.0454 89.0454 89.0454
8.90454e-8 89.0454 89.0454 89.0454
2.27641e-8 22.7641 22.7641 22.7641
2.27641e-8 22.7641 22.7641 22.7641
8.42157e-9 8.42157 8.42157 8.42157
8.42157e-9 8.42157 8.42157 8.42157
3.65085e-8 36.5085 36.5085 36.5085
3.65085e-8 36.5085 36.5085 36.5085Furthermore, we provide a simple function to get the mass of each cell in code units:
getmass(gas)37898393-element Vector{Float64}:
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
8.940696716308594e-7
⋮
1.0288910576564388e-7
1.0288910576564388e-7
1.9442336261293343e-7
1.9442336261293343e-7
8.90453891574347e-8
8.90453891574347e-8
2.276412192306883e-8
2.276412192306883e-8
8.421571563820485e-9
8.421571563820485e-9
3.650851622718898e-8
3.650851622718898e-8Get Multiple Quantities
Get several quantities with one function call by passing an array containing the selected variables. getvar returns a dictionary containing 1dim arrays for each quantity in code units:
quantities = getvar(gas, [:mass, :ekin])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…
:ekin => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 2.28274e-7, 2.…The units for each quantity can by passed as an array to the keyword argument "units" (plural, compare with single quantitiy call above) by preserving the order of the vars argument:
quantities = getvar(gas, [:mass, :ekin], units=[:Msol, :erg])Dict{Any, Any} with 2 entries:
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894…
:ekin => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 1.95354e49, 1.…The function can be called without any keywords by preserving the following order: dataobject, variables, units
quantities = getvar(gas, [:mass, :ekin], [:Msol, :erg])Dict{Any, Any} with 2 entries:
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894…
:ekin => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 1.95354e49, 1.…The arrays of the single quantities can be accessed from the dictionary:
quantities[:mass]37898393-element Vector{Float64}:
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
894.0696716308591
⋮
102.88910576564386
102.88910576564386
194.42336261293337
194.42336261293337
89.04538915743468
89.04538915743468
22.764121923068824
22.764121923068824
8.421571563820482
8.421571563820482
36.50851622718897
36.50851622718897If all selected variables should be of the same unit use the following arguments: dataobject, array of quantities, unit (no array needed):
quantities = getvar(gas, [:vx, :vy, :vz], :km_s)Dict{Any, Any} with 3 entries:
:vy => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … -97.5301, -97.53…
:vz => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0…
:vx => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … -24.307, -24.307…Get Quantities related to a center
Some quantities are related to a given center, e.g. radius in cylindrical coordinates, see the overview :
getvar()Predefined vars that can be calculated for each cell/particle:
----------------------------------------------------------------
=============================[gas]:=============================
-all the non derived hydro vars-
:cpu, :level, :rho, :cx, :cy, :cz, :vx, :vy, :vz, :p, var6,...
-derived hydro vars-
:x, :y, :z
:mass, :cellsize, :volume, :freefall_time
:cs, :mach, :machx, :machy, :machz, :jeanslength, :jeansnumber
:T, :Temp, :Temperature with p/rho
:h, :hx, :hy, :hz (specific angular momentum)
==========================[particles]:==========================
-all the non derived particle vars-
:cpu, :level, :id, :family, :tag
:x, :y, :z, :vx, :vy, :vz, :mass, :birth, :metal....
-derived particle vars-
:age
===========================[gravity]:===========================
-all the non derived gravity vars-
:cpu, :level, cx, cy, cz, :epot, :ax, :ay, :az
-derived gravity vars-
:x, :y, :z
:cellsize, :volume
===========================[clumps]:===========================
:peak_x or :x, :peak_y or :y, :peak_z or :z
:v, :ekin,...
=====================[gas or particles]:=======================
:v, :ekin
related to a given center:
---------------------------
:r_cylinder, :r_sphere (radial components)
:ϕ
:vr_cylinder
:vϕ_cylinder
----------------------------------------------------------------The unit of the provided center-array (in cartesian coordinates: x,y.z) is given by the keyword argument center_unit (default: code units). The function returns the quantitites in code units:
cv = (gas.boxlen / 2.) * gas.scale.kpc # provide the box-center in kpc
# e.g. for :mass the center keyword is ignored
quantities = getvar(gas, [:mass, :r_cylinder], center=[cv, cv, cv], center_unit=:kpc) Dict{Any, Any} with 2 entries:
:r_cylinder => [70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583…
:mass => [8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407e-7, 8.9407…Here, the function returns the result in the units that are provided. Note: E.g. the quantities :mass and :v (velocity) are not affected by the given center.
quantities = getvar(gas, [:mass, :r_cylinder, :v], units=[:Msol, :kpc, :km_s], center=[cv, cv, cv], center_unit=:kpc)Dict{Any, Any} with 3 entries:
:r_cylinder => [70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583…
:v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 100.513,…
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.0…Use the short notation for the box center :bc or :boxcenter for all dimensions (x,y,z). In this case the keyword center_unit is ignored:
quantities = getvar(gas, [:mass, :r_cylinder, :v], units=[:Msol, :kpc, :km_s], center=[:boxcenter])Dict{Any, Any} with 3 entries:
:r_cylinder => [70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583…
:v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 100.513,…
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.0…quantities = getvar(gas, [:mass, :r_cylinder, :v], units=[:Msol, :kpc, :km_s], center=[:bc])Dict{Any, Any} with 3 entries:
:r_cylinder => [70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583, 70.1583…
:v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 100.513,…
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.0…Use the box center notation for individual dimensions, here x,z. The keyword center_unit is needed for the y-coordinates:
quantities = getvar(gas, [:mass, :r_cylinder, :v], units=[:Msol, :kpc, :km_s], center=[:bc, 24., :bc], center_unit=:kpc)Dict{Any, Any} with 3 entries:
:r_cylinder => [54.9408, 54.9408, 54.9408, 54.9408, 54.9408, 54.9408, 54.9408…
:v => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 100.513,…
:mass => [894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.07, 894.0…Create Costum Quantities
Example1: Represent the positions of the data as the radius for a disk, centred in the simulation box (cylindrical coordinates):
boxlen = info.boxlen
cv = boxlen / 2. # box-center
levels = getvar(gas, :level) # get the level of each cell
cellsize = boxlen ./ 2 .^levels # calculate the cellsize for each cell (code units)
# or use the predefined quantity
cellsize = getvar(gas, :cellsize)
# convert the cell-number (related to the levels) into positions (code units), relative to the box center
x = getvar(gas, :cx) .* cellsize .- cv # (code units)
y = getvar(gas, :cy) .* cellsize .- cv # (code units)
# or use the predefined quantity
x = getvar(gas, :x, center=[:bc])
y = getvar(gas, :y, center=[:bc])
# calculate the cylindrical radius and scale from code units to kpc
radius = sqrt.(x.^2 .+ y.^2) .* info.scale.kpc37898393-element Vector{Float64}:
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
70.15825094589823
⋮
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808
20.08587520654808Use IndexedTables Functions
using Mera.IndexedTablesExample: Get the mass for each gas cell: mi = ρi * cellvolumei = ρ_i * (boxlen / 2^level)^3
Version 1
Use the select function and calculate the mass for each cell:
boxlen = gas.boxlen
level = select(gas.data, :level ) # get level information from each cell
cellvol = (boxlen ./ 2 .^level).^3 # calculate volume for each cell
mass1 = select(gas.data, :rho) .* cellvol .* info.scale.Msol; # calculate the mass for each cell in Msol unitsVersion 2
Use a single time the select function to do the calculations from above :
mass2 = select( gas.data, (:rho, :level)=>p->p.rho * (boxlen / 2^p.level)^3 ) .* info.scale.Msol;Version 3
Use the map function to do the calculations from above :
mass3 = map(p->p.rho * (boxlen / 2^p.level)^3, gas.data) .* info.scale.Msol;Comparison of the results:
[mass1 mass2 mass3]37898393×3 Matrix{Float64}:
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
894.07 894.07 894.07
⋮
102.889 102.889 102.889
102.889 102.889 102.889
194.423 194.423 194.423
194.423 194.423 194.423
89.0454 89.0454 89.0454
89.0454 89.0454 89.0454
22.7641 22.7641 22.7641
22.7641 22.7641 22.7641
8.42157 8.42157 8.42157
8.42157 8.42157 8.42157
36.5085 36.5085 36.5085
36.5085 36.5085 36.5085Statistical Quantities
info = getinfo(400, "/Volumes/FASTStorage/Simulations/Mera-Tests/manu_sim_sf_L14", verbose=false);
gas = gethydro(info, [:rho, :vx, :vy, :vz], lmax=8, smallr=1e-5, verbose=false);
particles = getparticles(info, [:mass, :vx, :vy, :vz], verbose=false)
clumps = getclumps(info, verbose=false);Reading data...
Progress: 100%|█████████████████████████████████████████| Time: 0:02:47
Progress: 100%|█████████████████████████████████████████| Time: 0:00:19Pass any kind of Array{<:Real,1} (Float, Integer,...) to the wstat function to get several unweighted statistical quantities at once:
stats_gas = wstat( getvar(gas, :vx, :km_s) )
stats_particles = wstat( getvar(particles, :vx, :km_s) )
stats_clumps = wstat( getvar(clumps, :rho_av, :Msol_pc3) );The result is an object that contains several fields with the statistical quantities:
println( typeof(stats_gas) )
println( typeof(stats_particles) )
println( typeof(stats_clumps) )
propertynames(stats_gas)Mera.WStatType
Mera.WStatType
Mera.WStatType
(:mean, :median, :std, :skewness, :kurtosis, :min, :max)println( "Gas <vx>_allcells : ", stats_gas.mean, " km/s" )
println( "Particles <vx>_allparticles : ", stats_particles.mean, " km/s" )
println( "Clumps <rho_av>_allclumps : ", stats_clumps.mean, " Msol/pc^3" )Gas <vx>_allcells : -2.9318774650713717 km/s
Particles <vx>_allparticles : -11.594477384589647 km/s
Clumps <rho_av>_allclumps : 594.7315900915926 Msol/pc^3println( "Gas min/max_allcells : ", stats_gas.min, "/", stats_gas.max, " km/s" )
println( "Particles min/max_allparticles : ", stats_particles.min,"/", stats_particles.max, " km/s" )
println( "Clumps min/max_allclumps : ", stats_clumps.min, "/", stats_clumps.max, " Msol/pc^3" )Gas min/max_allcells : -676.5464963488397/894.9181733956399 km/s
Particles min/max_allparticles : -874.6440509326601/670.7956741234592 km/s
Clumps min/max_allclumps : 125.4809686796669/5357.370234867635 Msol/pc^3Weighted Statistics
Pass any kind of Array{<:Real,1} (Float, Integer,...) for the given variables and one for the weighting with the same length. The weighting goes cell by cell, particle by particle, clump by clump, etc...:
stats_gas = wstat( getvar(gas, :vx, :km_s), weight=getvar(gas, :mass ));
stats_particles = wstat( getvar(particles, :vx, :km_s), weight=getvar(particles, :mass ));
stats_clumps = wstat( getvar(clumps, :peak_x, :kpc ), weight=getvar(clumps, :mass_cl)) ;Without the keyword weight the following order for the given arrays has to be maintained: values, weight
stats_gas = wstat( getvar(gas, :vx, :km_s), getvar(gas, :mass ));
stats_particles = wstat( getvar(particles, :vx, :km_s), getvar(particles, :mass ));
stats_clumps = wstat( getvar(clumps, :peak_x, :kpc ), getvar(clumps, :mass_cl)) ;propertynames(stats_gas)(:mean, :median, :std, :skewness, :kurtosis, :min, :max)println( "Gas <vx>_allcells : ", stats_gas.mean, " km/s (mass weighted)" )
println( "Particles <vx>_allparticles : ", stats_particles.mean, " km/s (mass weighted)" )
println( "Clumps <peak_x>_allclumps : ", stats_clumps.mean, " kpc (mass weighted)" )Gas <vx>_allcells : -1.1999253584798584 km/s (mass weighted)
Particles <vx>_allparticles : -11.623422700314533 km/s (mass weighted)
Clumps <peak_x>_allclumps : 23.135765457064576 kpc (mass weighted)println( "Gas min/max_allcells : ", stats_gas.min, "/", stats_gas.max, " km/s" )
println( "Particles min/max_allparticles : ", stats_particles.min,"/", stats_particles.max, " km/s" )
println( "Clumps min/max_allclumps : ", stats_clumps.min, "/", stats_clumps.max, " Msol/pc^3" )Gas min/max_allcells : -676.5464963488397/894.9181733956399 km/s
Particles min/max_allparticles : -874.6440509326601/670.7956741234592 km/s
Clumps min/max_allclumps : 10.29199219000667/38.17382813002474 Msol/pc^3For the average of the gas-density use volume weighting:
stats_gas = wstat( getvar(gas, :rho, :g_cm3), weight=getvar(gas, :volume) );println( "Gas <rho>_allcells : ", stats_gas.mean, " g/cm^3 (volume weighted)" )Gas <rho>_allcells : 1.8958545012297553e-26 g/cm^3 (volume weighted)Helpful Functions
Get the x,y,z positions of every cell relative to a given center:
x,y,z = getpositions(gas, :kpc, center=[24.,24.,24.], center_unit=:kpc); # returns a Tuple of 3 arraysThe box-center can be calculated automatically:
x,y,z = getpositions(gas, :kpc, center=[:boxcenter]);[x y z] # preview of the output849332×3 Matrix{Float64}:
-23.25 -23.25 -23.25
-23.25 -23.25 -22.5
-23.25 -23.25 -21.75
-23.25 -23.25 -21.0
-23.25 -23.25 -20.25
-23.25 -23.25 -19.5
-23.25 -23.25 -18.75
-23.25 -23.25 -18.0
-23.25 -23.25 -17.25
-23.25 -23.25 -16.5
-23.25 -23.25 -15.75
-23.25 -23.25 -15.0
-23.25 -23.25 -14.25
⋮
16.125 3.9375 0.1875
16.125 3.9375 0.375
16.125 3.9375 0.5625
16.125 3.9375 0.75
16.125 4.125 -0.5625
16.125 4.125 -0.375
16.125 4.125 -0.1875
16.125 4.125 0.0
16.125 4.125 0.1875
16.125 4.125 0.375
16.125 4.125 0.5625
16.125 4.125 0.75Get the extent of the dataset-domain:
getextent(gas) # returns Tuple of (xmin, xmax), (ymin ,ymax ), (zmin ,zmax )((0.0, 48.0), (0.0, 48.0), (0.0, 48.0))Get the extent relative to a given center:
getextent(gas, center=[:boxcenter])((-24.0, 24.0), (-24.0, 24.0), (-24.0, 24.0))Get simulation time in code unit oder physical unit
gettime(info)39.9019537349027gettime(info, :Myr)594.9774920106152gettime(gas, :Myr)594.9774920106152