# Comoving units and code units¶

Notebook

In yt 3.0, we want to make it easier to access "raw" simulation data that a code writes directly to disk. The new unit system makes it much simpler to convert back and forth between physical coordinates and the unscaled "raw" coordinate system used internally in the simulation code. In some cases, this conversion involves transforming to comoving coordinates, so that is also covered here.

### Code units¶

Let's take a look at a cosmological enzo dataset to play with converting between physical units and code units:

In [1]:
import yt


The conversion factors between Enzo's internal unit system and the physical CGS system are stored in the dataset's unit_registry object. Code units have names like code_length and code_time. Let's take a look at the names of all of the code units, along with their CGS conversion factors for this cosmological enzo dataset:

In [2]:
reg = ds.unit_registry

for un in reg.keys():
if un.startswith('code_'):
fmt_tup = (un, reg.lut[un][0], str(reg.lut[un][1]))
print ("Unit name:      {:<15}\nCGS conversion: {:<15}\nDimensions:     {:<15}\n".format(*fmt_tup))

Unit name:      code_magnetic
CGS conversion: 0.005133186931548846
Dimensions:     sqrt((mass))/(sqrt((length))*(time))

Unit name:      code_temperature
CGS conversion: 1.0
Dimensions:     (temperature)

Unit name:      code_mass
CGS conversion: 4.891307900838622e+50
Dimensions:     (mass)

Unit name:      code_pressure
CGS conversion: 2.0968351867734223e-06
Dimensions:     (mass)/((length)*(time)**2)

Unit name:      code_time
CGS conversion: 648009786344140.2
Dimensions:     (time)

Unit name:      code_length
CGS conversion: 5.5551728502644205e+26
Dimensions:     (length)

Unit name:      code_metallicity
CGS conversion: 1.0
Dimensions:     1

Unit name:      code_density
CGS conversion: 2.853200371985856e-30
Dimensions:     (mass)/(length)**3

Unit name:      code_velocity
CGS conversion: 8588289701.681004
Dimensions:     (length)/(time)

Unit name:      code_specific_energy
CGS conversion: 7.375871999999998e+19
Dimensions:     (length)**2/(time)**2


In [3]:
fmt_tup

Out[3]:
('code_specific_energy', 7.375871999999998e+19, '(length)**2/(time)**2')

Most of the time you will not have to deal with the unit registry. For example, the conversion factors to code units are stored as attributes of the dataset object:

In [4]:
print ("Length unit: ", ds.length_unit)
print ("Time unit: ", ds.time_unit)
print ("Mass unit: ", ds.mass_unit)
print ("Velocity unit: ", ds.velocity_unit)

Length unit:  128.0 Mpccm/h
Time unit:  648009786344140.2 s
Mass unit:  4.891307900838622e+50 g
Velocity unit:  8588289701.681004 cm/s


Conversion factors will be supplied in CGS by default. We can also ask what the conversion factors are in code units.

In [5]:
print ("Length unit: ", ds.length_unit.in_units('code_length'))
print ("Time unit: ", ds.time_unit.in_units('code_time'))
print ("Mass unit: ", ds.mass_unit.in_units('code_mass'))
print ("Velocity unit: ", ds.velocity_unit.in_units('code_velocity'))

Length unit:  1.0 code_length
Time unit:  1.0 code_time
Mass unit:  1.0 code_mass
Velocity unit:  1.0 code_velocity


as expected, all the conversion factors are unity in code units.

We can also play with unit conversions on ds.domain_width. First, we see for enzo how code length units are defined relative to the domain width:

In [6]:
ds.domain_width

Out[6]:
YTArray([ 1.,  1.,  1.]) code_length
In [7]:
ds.domain_width.in_cgs()

Out[7]:
YTArray([  5.55517285e+26,   5.55517285e+26,   5.55517285e+26]) cm
In [8]:
ds.domain_width.in_units('Mpccm/h')

Out[8]:
YTArray([ 128.,  128.,  128.]) Mpccm/h

### Comoving units¶

This last example uses a cosmological unit. In english, I asked for the domain width in comoving megaparsecs, scaled as if the hubble constant were 100 km/s/Mpc. Although $h$ isn't really a unit, yt treats it as one for the purposes of the unit system.

As an aside, Darren Croton's research note on the history, use, and interpretation of $h$ as it appears in the astronomical literature is pretty much required reading for anyone who has to deal with factors of $h$ every now and then.

In yt, comoving length unit symbols are named following the pattern (length symbol)cm, i.e. pccm for comoving parsec or mcm for a comoving meter. A comoving length unit is different from the normal length unit by a factor of $(1+z)$:

In [9]:
z = ds.current_redshift

print (ds.quan(1, 'Mpc')/ds.quan(1, 'Mpccm'))
print (1+z)

1.0013930880640793 dimensionless
1.0013930880640796


As we saw before, $h$ is treated like any other unit symbol. It has dimensionless units, just like a scalar:

In [10]:
print (ds.quan(1, 'Mpc')/ds.quan(1, 'Mpc/h'))
print (ds.hubble_constant)

0.71 dimensionless
0.71


These units can be used in readily used in plots and anywhere a length unit is appropriate in yt.

In [11]:
slc = yt.SlicePlot(ds, 0, 'density', width=(128, 'Mpccm/h'))
slc.set_figure_size(6)

Out[11]:

### The unit registry¶

When you create a YTArray without referring to a unit registry, yt uses the default unit registry, which does not include code units or comoving units.

In [12]:
from yt import YTQuantity

a = YTQuantity(3, 'cm')

print (a.units.registry.keys())

dict_keys(['Myr', 'g', 'fs', 'lbf', 'ns', 'au', 'kpc', 'gauss', 'solMass', 'degree', 'kg', 'keV', 't_geom', 'dimensionless', 'angstrom', 'Mjup', 'esu', 'R_jup', 'day', 'sr', 'rad', 'mile', 'Zsun', 'Tsun', 'lon', 'km', 'q_pl', 'steradian', 's', 'yr', 'r_earth', 'GeV', 'lat', 'C', 'degC', 'lbm', 'photons', 'solRad', 'me', 'deg', 'G', 'statA', 'counts', 'Mm', 'Gyr', 'mol', 'h', 'eV', 'cm', 'msun', 'min', 'dyn', 'Lsun', 'c', 'Mpc', 'Jy', 'ms', 'mg', 'kyr', 'Msun', 'A', 'm_geom', 'Angstrom', 'hourangle', 'statV', 'V', 'nm', 'R_earth', 'erg', 'mp', 'm', 'um', 'Mearth', 'rsun', 'amu', 'ps', 'degF', 'Pa', 'uG', 'T', 'W', 'mas', 'ohm', 'R', 'Rsun', 'radian', 'K', 'J', 'fm', 'r_sun', 'ft', 'hr', 'pc', 'm_pl', 'MeV', 'l_pl', 'Sv', 'Gpc', 'mm', 'E_pl', 'Fr', 'statohm', 'rayleigh', 'Hz', 'arcsec', 'dyne', 'pm', 'l_geom', 'solLum', 'arcmin', 'r_jup', 'T_pl', 'N', 't_pl', 'atm', 'statC', 'd', 'AU', 'ly', 'R_sun'])


When a dataset is loaded, yt infers conversion factors from the internal simulation unit system to the CGS unit system. These conversion factors are stored in a unit_registry along with conversion factors to the other known unit symbols. For the cosmological Enzo dataset we loaded earlier, we can see there are a number of additional unit symbols not defined in the default unit lookup table:

In [13]:
print (sorted([k for k in ds.unit_registry.keys() if k not in a.units.registry.keys()]))

['AUcm', 'Mpccm', 'Ppc', 'Tpc', 'a', 'aucm', 'code_density', 'code_length', 'code_magnetic', 'code_mass', 'code_metallicity', 'code_pressure', 'code_specific_energy', 'code_temperature', 'code_time', 'code_velocity', 'mcm', 'pccm', 'unitary']


Since code units do not appear in the default unit symbol lookup table, one must explicitly refer to a unit registry when creating a YTArray to be able to convert to the unit system of a simulation.

To make this as clean as possible, there are array and quantity-creating convenience functions attached to the Dataset object:

• ds.arr()
• ds.quan()

These functions make it straightforward to create arrays and quantities that can be converted to code units or comoving units. For example:

In [14]:
a = ds.quan(3, 'code_length')

print (a)
print (a.in_cgs())
print (a.in_units('Mpccm/h'))

3.0 code_length
1.666551855079326e+27 cm
384.0 Mpccm/h

In [15]:
b = ds.arr([3, 4, 5], 'Mpccm/h')
print (b)
print (b.in_cgs())

[ 3.  4.  5.] Mpccm/h
[  1.30199364e+25   1.73599152e+25   2.16998939e+25] cm


### Overriding Code Unit Definitions¶

On occasion, you might have a dataset for a supported frontend that does not have the conversions to code units accessible (for example, Athena data) or you may want to change them outright. yt provides a mechanism so that one may provide their own code unit definitions to load, which override the default rules for a given frontend for defining code units. This is provided through the units_override dictionary. We'll use an example of an Athena dataset. First, a call to load without units_override:

In [16]:
ds1 = yt.load("MHDSloshing/virgo_low_res.0054.vtk")
print (ds1.length_unit)
print (ds1.mass_unit)
print (ds1.time_unit)
sp1 = ds1.sphere("c",(0.1,"unitary"))
print (sp1["density"])

1.0 cm
1.0 g
1.0 s
[ 0.05134981  0.05134912  0.05109047 ...,  0.14608461  0.14489453
0.14385277] g/cm**3


This is a galaxy cluster dataset, so it is not likely that the units of density are correct. We happen to know that the unit definitions are different, so we can override the units:

In [17]:
units_override = {"length_unit":(1.0,"Mpc"),
"time_unit":(1.0,"Myr"),
"mass_unit":(1.0e14,"Msun")}


units_override can take the following keys:

• length_unit
• time_unit
• mass_unit
• magnetic_unit
• temperature_unit

and the associated values can be (value, unit) tuples, YTQuantities, or floats (in the latter case they are assumed to have the corresponding cgs unit).

In [18]:
ds2 = yt.load("MHDSloshing/virgo_low_res.0054.vtk", units_override=units_override)
print (ds2.length_unit)
print (ds2.mass_unit)
print (ds2.time_unit)
sp2 = ds2.sphere("c",(0.1,"unitary"))
print (sp2["density"])

1.0 Mpc
100000000000000.0 Msun
1.0 Myr
[  3.47531683e-28   3.47527018e-28   3.45776515e-28 ...,   9.88689766e-28
9.80635384e-28   9.73584863e-28] g/cm**3


This option should be used very carefully, and only if you know that the dataset does not provide units or that the unit definitions generated are incorrect for some reason.

## Units for Cosmological Datasets¶

yt has additional capabilities to handle the comoving coordinate system used internally in cosmological simulations. Simulations that use comoving coordinates, all length units have three other counterparts corresponding to comoving units, scaled comoving units, and scaled proper units. In all cases ‘scaled’ units refer to scaling by the reduced Hubble parameter - i.e. the length unit is what it would be in a universe where Hubble’s parameter is 100 km/s/Mpc.

To access these different units, yt has a common naming system. Scaled units are denoted by dividing by the scaled Hubble parameter h (which is in itself a unit). Comoving units are denoted by appending cm to the end of the unit name.

Using the parsec as an example,

pc
Proper parsecs, $\rm{pc}$.
pccm
Comoving parsecs, $\rm{pc}/(1+z)$.
pccm/h
Comoving parsecs normalized by the scaled hubble constant, $\rm{pc}/h/(1+z)$.
pc/h
Proper parsecs, normalized by the scaled hubble constant, $\rm{pc}/h$.