Units

In synthesizer all quantities a user interacts with (that are not dimensionless) have units associated with them, implemented with the unyt package.

Synthesizer objects and methods should always be provided with quantites and associated units. This can be easily achieved with the unyt package.

from unyt import Mpc

# Define a variable with units
x = 1 * Mpc

print(x)
print("x is now a unyt_quantity: type(x)=", type(x))

1 Mpc
x is now a unyt_quantity: type(x)= <class 'unyt.array.unyt_quantity'>

All unit functionality in synthesizer is contained in the units module. In this module we have defined a dictionary containing the default units of all attributes throughout synthesizer.

from synthesizer.units import default_units

print(default_units)

{'lam': Å, 'obslam': Å, 'wavelength': Å, 'vacuum_wavelength': Å, 'original_lam': Å, 'lam_min': Å, 'lam_max': Å, 'lam_eff': Å, 'lam_fwhm': Å, 'mean_lams': Å, 'pivot_lams': Å, 'nu': Hz, 'obsnu': Hz, 'nuz': Hz, 'original_nu': Hz, 'luminosity': erg/s, 'luminosities': erg/s, 'bolometric_luminosity': erg/s, 'bolometric_luminosities': erg/s, 'lnu': erg/(Hz*s), 'llam': erg/(s*Å), 'continuum': erg/(Hz*s), 'flux': erg/(cm**2*s), 'fnu': nJy, 'flam': erg/(cm**2*s*Å), 'equivalent_width': Å, 'coordinates': Mpc, 'radii': Mpc, 'smoothing_lengths': Mpc, 'softening_length': Mpc, 'velocities': km/s, 'mass': unyt_quantity(1., 'Msun'), 'masses': unyt_quantity(1., 'Msun'), 'initial_masses': unyt_quantity(1., 'Msun'), 'initial_mass': unyt_quantity(1., 'Msun'), 'current_masses': unyt_quantity(1., 'Msun'), 'dust_masses': unyt_quantity(1., 'Msun'), 'ages': yr, 'accretion_rate': unyt_quantity(1., 'Msun/yr'), 'accretion_rates': unyt_quantity(1., 'Msun/yr'), 'bb_temperature': K, 'bb_temperatures': K, 'inclination': degree, 'inclinations': degree, 'resolution': Mpc, 'fov': Mpc, 'orig_resolution': Mpc, 'centre': Mpc, 'photo_lnu': erg/(Hz*s), 'photo_fnu': erg/(Hz*cm**2*s), 'softening_lengths': Mpc}

As you can see, most units are defined symbollically, however, you may notice units such as Msun (the default used for masses) is defined as a compound unit in terms of kgs. Shortly we will cover working with the quantities returned by synthesizer and how compound units work in practice.

The Units object

The unit system is defined by the Units object. This object contains a collection of attributes defining the units associated to each quantity throughout synthesizer which is not dimensionless. Importantly, Units is a Singleton object. This means there can only ever be one instance of Units; if a second is instantiated then the first is returned. This ensures that the unit system remains consistent when running synthesizer.

from synthesizer.units import Units

# Define multiple Units instances
units1 = Units()
units2 = Units()

print("Both units instances are the same object:", units1 is units2)

Both units instances are the same object: True

You can take a look at the unit system by printing the instance of Units.

print(units1)

Unit System:
lam:                Å
obslam:             Å
wavelength:         Å
vacuum_wavelength:  Å
original_lam:       Å
lam_min:            Å
lam_max:            Å
lam_eff:            Å
lam_fwhm:           Å
mean_lams:          Å
pivot_lams:         Å
nu:                 Hz
obsnu:              Hz
nuz:                Hz
original_nu:        Hz
luminosity:         erg/s
luminosities:       erg/s
bolometric_luminosity: erg/s
bolometric_luminosities: erg/s
lnu:                erg/(Hz*s)
llam:               erg/(s*Å)
continuum:          erg/(Hz*s)
flux:               erg/(cm**2*s)
fnu:                nJy
flam:               erg/(cm**2*s*Å)
equivalent_width:   Å
coordinates:        Mpc
radii:              Mpc
smoothing_lengths:  Mpc
softening_length:   Mpc
velocities:         km/s
mass:               0.9999999999999999 Msun
masses:             0.9999999999999999 Msun
initial_masses:     0.9999999999999999 Msun
initial_mass:       0.9999999999999999 Msun
current_masses:     0.9999999999999999 Msun
dust_masses:        0.9999999999999999 Msun
ages:               yr
accretion_rate:     0.9999999999999999 Msun/yr
accretion_rates:    0.9999999999999999 Msun/yr
bb_temperature:     K
bb_temperatures:    K
inclination:        degree
inclinations:       degree
resolution:         Mpc
fov:                Mpc
orig_resolution:    Mpc
centre:             Mpc
photo_lnu:          erg/(Hz*s)
photo_fnu:          erg/(Hz*cm**2*s)
softening_lengths:  Mpc

Modifying the default unit system

If the default unit system works for your needs then you don’t need to do anything. You will never interact with the Units object and all quantites will have the default units associated to them automatically. However, if you need to change one or more of the units used you can import Units and instantiate it with a dictionary of the modified quantities.

from unyt import Msun, Myr, kpc

# Make the dictionary containing the units we want to change
new_units = {
    "coordinates": kpc,
    "smoothing_lengths": kpc,
    "softening_length": kpc,
    "ages": Myr,
}

# Set up the modified unit system
units = Units(new_units)

print()
print(units)

Unit System:
lam:                Å
obslam:             Å
wavelength:         Å
vacuum_wavelength:  Å
original_lam:       Å
lam_min:            Å
lam_max:            Å
lam_eff:            Å
lam_fwhm:           Å
mean_lams:          Å
pivot_lams:         Å
nu:                 Hz
obsnu:              Hz
nuz:                Hz
original_nu:        Hz
luminosity:         erg/s
luminosities:       erg/s
bolometric_luminosity: erg/s
bolometric_luminosities: erg/s
lnu:                erg/(Hz*s)
llam:               erg/(s*Å)
continuum:          erg/(Hz*s)
flux:               erg/(cm**2*s)
fnu:                nJy
flam:               erg/(cm**2*s*Å)
equivalent_width:   Å
coordinates:        Mpc
radii:              Mpc
smoothing_lengths:  Mpc
softening_length:   Mpc
velocities:         km/s
mass:               0.9999999999999999 Msun
masses:             0.9999999999999999 Msun
initial_masses:     0.9999999999999999 Msun
initial_mass:       0.9999999999999999 Msun
current_masses:     0.9999999999999999 Msun
dust_masses:        0.9999999999999999 Msun
ages:               yr
accretion_rate:     0.9999999999999999 Msun/yr
accretion_rates:    0.9999999999999999 Msun/yr
bb_temperature:     K
bb_temperatures:    K
inclination:        degree
inclinations:       degree
resolution:         Mpc
fov:                Mpc
orig_resolution:    Mpc
centre:             Mpc
photo_lnu:          erg/(Hz*s)
photo_fnu:          erg/(Hz*cm**2*s)
softening_lengths:  Mpc

Something has gone wrong… but recall that the unit system will return the original if one exists, so actually this should be completely expected.

This issue highlights the need to set up Units before doing anything else. If any computations have been done the Units instance will exist and will not be modifiable after the fact. However, should you fall in this trap the code will warn you as above - no hidden gotchas here!

Now, lets go against the advice above, and use the highly inadvisable force argument to get a new Unit system. But please note, in a real use case, forcing a modified unit system WILL NOT convert existing quantities to the new unit system.

# Set up the modified unit system
units = Units(new_units, force=True)

print()
print(units)

Redefining unit system:
coordinates: kpc
smoothing_lengths: kpc
softening_length: kpc
ages: Myr

Unit System:
lam:                Å
obslam:             Å
wavelength:         Å
vacuum_wavelength:  Å
original_lam:       Å
lam_min:            Å
lam_max:            Å
lam_eff:            Å
lam_fwhm:           Å
mean_lams:          Å
pivot_lams:         Å
nu:                 Hz
obsnu:              Hz
nuz:                Hz
original_nu:        Hz
luminosity:         erg/s
luminosities:       erg/s
bolometric_luminosity: erg/s
bolometric_luminosities: erg/s
lnu:                erg/(Hz*s)
llam:               erg/(s*Å)
continuum:          erg/(Hz*s)
flux:               erg/(cm**2*s)
fnu:                nJy
flam:               erg/(cm**2*s*Å)
equivalent_width:   Å
coordinates:        kpc
radii:              Mpc
smoothing_lengths:  kpc
softening_length:   kpc
velocities:         km/s
mass:               0.9999999999999999 Msun
masses:             0.9999999999999999 Msun
initial_masses:     0.9999999999999999 Msun
initial_mass:       0.9999999999999999 Msun
current_masses:     0.9999999999999999 Msun
dust_masses:        0.9999999999999999 Msun
ages:               Myr
accretion_rate:     0.9999999999999999 Msun/yr
accretion_rates:    0.9999999999999999 Msun/yr
bb_temperature:     K
bb_temperatures:    K
inclination:        degree
inclinations:       degree
resolution:         Mpc
fov:                Mpc
orig_resolution:    Mpc
centre:             Mpc
photo_lnu:          erg/(Hz*s)
photo_fnu:          erg/(Hz*cm**2*s)
softening_lengths:  Mpc

Working with Quantity objects

There is no need to work with the Units object itself beyond initially defining a modified unit system. Beyond this, all unit operations are handled “behind the scenes”. This hidden functionality is enabled by the Quantity object.

All attributes on synthesizer objects which carry units are in fact Quantity objects. Quantitiy objects carry a reference to the global unit system, and extract the appropriate units depending on the name of the variable storing the Quantity. As such, a user will never instantiate a quantity themselves, but their usage is important.

One simple thing to keep in mind is how to return the value with or without units. This is achieved by the application or omission of a leading underscore to a variable name.

Lets create an Sed object, which has a wavelength array stored under lam.

import numpy as np
from unyt import Hz, angstrom, erg, s

from synthesizer.sed import Sed

# Make an sed with arbitrary arguments
sed = Sed(
    lam=np.linspace(10, 1000, 10) * angstrom, lnu=np.ones(10) * erg / s / Hz
)

We can access this attribute with units as you would expect to access any attribute.

print(sed.lam)

[  10.  120.  230.  340.  450.  560.  670.  780.  890. 1000.] Å

Or we can append a leading underscore and return it without units.

print(sed._lam)

[  10.  120.  230.  340.  450.  560.  670.  780.  890. 1000.]

In the case of compound units this is somewhat less elegant. Lets demonstrate with a Stars object.

from synthesizer.particle.stars import Stars

# Create a dummy Stars object
stars = Stars(
    initial_masses=np.random.rand(10) * Msun,
    ages=np.ones(10) * Myr,
    metallicities=np.ones(10),
)

If we print the initial_masses with units we get the compound version in kg.

print(stars.initial_masses)

[0.96109166 0.61018815 0.83004024 0.10923202 0.88023096 0.02437266
 0.09106734 0.48094055 0.59908529 0.39421454] Msun

However, if we extract the values alone we get the values we expect in \(M_\odot\).

print(stars._initial_masses)

[0.96109166 0.61018815 0.83004024 0.10923202 0.88023096 0.02437266
 0.09106734 0.48094055 0.59908529 0.39421454]

Its worth keeping this in mind whenever extracting masses from synthesizer objects.

Automatic unit conversion

Finally, let’s utilise some automatic unit conversion. If we input a mixture of properties to a synthesizer object, all with different units to the global unit system, we don’t have to convert them all before inputting them. As long as we pass them to synthesizer with unyt units attached, the conversion will be handled automatically. Here we use a Stars object again to demonstrate.

from unyt import Mpc, g, m

# Create a dummy Stars object
stars = Stars(
    initial_masses=np.random.rand(10) * 10**34.0 * g,
    ages=np.ones(10) * Myr,
    metallicities=np.ones(10),
    coordinates=np.random.rand(10, 3) * Mpc,
    smoothing_lengths=np.random.rand(10) * 10**22.0 * m,
)

print(
    "stars.initial_masses[0]=",
    stars.initial_masses[0],
    "=",
    stars._initial_masses[0],
    "Msun",
)
print("stars.ages[0]=", stars.ages[0])
print("stars.coordinates[0]=", stars.coordinates[0])
print("stars.smoothing_lengths[0]=", stars.smoothing_lengths[0])

stars.initial_masses[0]= 1.7033893800504492 Msun = 1.7033893800504494 Msun
stars.ages[0]= 1000000.0 yr
stars.coordinates[0]= [0.20752606 0.19493777 0.79362484] Mpc
stars.smoothing_lengths[0]= 0.09713545106505901 Mpc