diff --git a/CHANGELOG.md b/CHANGELOG.md
index 69afd52..587db06 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,12 @@
+# 0.2.0
+
+* Adds lookback time and distance calculations.
+* Adds comoving volume calculations.
+* Updates values for constants based on CODATA 2018.
+* Refactors handling of base units.
+* Add initial benchmarks.
+* Performance improvements.
+
 # 0.1.0
 
 Initial release.
diff --git a/Cargo.toml b/Cargo.toml
index b4624ca..2f57587 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -13,3 +13,16 @@ license = "MIT OR Apache-2.0"
 [dependencies]
 anyhow = "1"
 once_cell = "1.8"
+
+[dev-dependencies]
+criterion = { version = "0.3", features=["html_reports"] }
+
+[[bench]]
+name = "distances"
+harness = false
+
+[profile.release]
+debug = 1
+
+[rust]
+debuginfo-level = 1
diff --git a/README.md b/README.md
index 2e21b71..46b50af 100644
--- a/README.md
+++ b/README.md
@@ -16,8 +16,7 @@ A library for computing quantities in cosmology in the Rust programming language
 
 ```rust
 let cosmology = FLRWCosmology::two_component(0.286, 0.714, 69.6);
-assert!(cosmology.radial_comoving_distance(2.0).0 > 5273.);
-assert!(cosmology.radial_comoving_distance(2.0).0 < 5274.); 
+assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) > Mpc::new(5273.));
 ```
 
 # Developers
@@ -30,3 +29,13 @@ This project requires [`Rust`](https://www.rust-lang.org/tools/install). Once in
 cargo build
 cargo test
 ```
+
+## Benchmarks
+
+Run `criterion` benchmarks using:
+
+```
+cargo bench
+```
+
+This will generate a report at `target/criterion/report/index.html`.
diff --git a/benches/distances.rs b/benches/distances.rs
new file mode 100644
index 0000000..1263f70
--- /dev/null
+++ b/benches/distances.rs
@@ -0,0 +1,25 @@
+use cosmocalc::{cosmology::OmegaFactors, Distances, FLRWCosmology, FloatingPointUnit, Redshift};
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
+
+pub fn bench_luminosity_distance(c: &mut Criterion) {
+    let mut group = c.benchmark_group("d_L");
+    let omegas = OmegaFactors::new(0.27, 0.73, 0.044).unwrap();
+    let cosmology = FLRWCosmology::new(None, None, 70.0, omegas, None, None, None).unwrap();
+
+    let z_1 = Redshift::new(1.);
+    group.bench_with_input(
+        BenchmarkId::new("d_L(z=1)", format!("{:?}", z_1)),
+        &z_1,
+        |b, z| b.iter(|| cosmology.luminosity_distance(*z)),
+    );
+    let z_2 = Redshift::new(2.);
+    group.bench_with_input(
+        BenchmarkId::new("d_L(z=2)", format!("{:?}", z_2)),
+        &z_2,
+        |b, z| b.iter(|| cosmology.luminosity_distance(*z)),
+    );
+    group.finish();
+}
+
+criterion_group!(benches, bench_luminosity_distance);
+criterion_main!(benches);
diff --git a/benchmarks/astropy_vectors.py b/benchmarks/astropy_vectors.py
deleted file mode 100644
index 4d0a5eb..0000000
--- a/benchmarks/astropy_vectors.py
+++ /dev/null
@@ -1,42 +0,0 @@
-from astropy.cosmology import FLRW, FlatLambdaCDM, LambdaCDM, Planck18
-
-def flat_universe_distances_no_relativistic_contribution():
-    H0 = 69.6
-    Om0 = 0.286
-    Ode0 = 0.714
-    Ob0 = 0.05
-    cosmo = FlatLambdaCDM(H0, Om0, Ode0, 0, Ob0=Ob0)
-
-    print('comoving transverse distance to z=3')
-    print(cosmo.comoving_transverse_distance(3))
-    # 6482.549296743232 Mpc
-    print('angular diameter distance to z=3')
-    print(cosmo.angular_diameter_distance(3))
-    # 1620.637324185808 Mpc
-    print('luminosity distance to z=3')
-    print(cosmo.luminosity_distance(3))
-    # 25930.197186972928 Mpc
-
-
-def flat_universe_distances_with_radiation_but_no_neutrinos():
-    H0 = 69.6
-    Om0 = 0.299
-    Ode0 = 0.7
-    Ob0 = 0.05
-    Tcmb0 = 2.7255
-    cosmo = FlatLambdaCDM(H0, Om0, Ode0, 2.7255, 0, Ob0=Ob0)
-
-    print('comoving transverse distance to z=3')
-    print(cosmo.comoving_transverse_distance(3))
-    # 6398.504909397802 Mpc
-    print('angular diameter distance to z=3')
-    print(cosmo.angular_diameter_distance(3))
-    # 1599.6262273494506 Mpc
-    print('luminosity distance to z=3')
-    print(cosmo.luminosity_distance(3))
-    # 25594.01963759121 Mpc
-
-
-if __name__=="__main__":
-    flat_universe_distances_no_relativistic_contribution()
-    flat_universe_distances_with_radiation_but_no_neutrinos()
\ No newline at end of file
diff --git a/benchmarks/README.md b/python_benchmarks/README.md
similarity index 100%
rename from benchmarks/README.md
rename to python_benchmarks/README.md
diff --git a/python_benchmarks/astropy_vectors.py b/python_benchmarks/astropy_vectors.py
new file mode 100644
index 0000000..2e756ce
--- /dev/null
+++ b/python_benchmarks/astropy_vectors.py
@@ -0,0 +1,86 @@
+from astropy.cosmology import FLRW, FlatLambdaCDM, LambdaCDM, Planck18
+import timeit
+
+def flat_universe_distances_no_relativistic_contribution():
+    H0 = 69.6
+    Om0 = 0.286
+    Ode0 = 0.714
+    Ob0 = 0.05
+    cosmo = FlatLambdaCDM(H0, Om0, Ode0, 0, Ob0=Ob0)
+
+    print('comoving transverse distance to z=3')
+    print(cosmo.comoving_transverse_distance(3))
+    # 6482.549296743232 Mpc
+    print('angular diameter distance to z=3')
+    print(cosmo.angular_diameter_distance(3))
+    # 1620.637324185808 Mpc
+    print('luminosity distance to z=3')
+    print(cosmo.luminosity_distance(3))
+    # 25930.197186972928 Mpc
+
+
+def flat_universe_distances_with_radiation_but_no_neutrinos():
+    H0 = 69.6
+    Om0 = 0.299
+    Ode0 = 0.7
+    Ob0 = 0.05
+    Tcmb0 = 2.7255
+    cosmo = FlatLambdaCDM(H0, Om0, Ode0, 2.7255, 0, Ob0=Ob0)
+
+    print('comoving transverse distance to z=3')
+    print(cosmo.comoving_transverse_distance(3))
+    # 6398.504909397802 Mpc
+    print('angular diameter distance to z=3')
+    print(cosmo.angular_diameter_distance(3))
+    # 1599.6262273494506 Mpc
+    print('luminosity distance to z=3')
+    print(cosmo.luminosity_distance(3))
+    # 25594.01963759121 Mpc
+
+
+def lookback_time():
+   Om0 = 0.27
+   Ode0 = 0.73
+   Ob0 = 0.044
+   H0 = 70.0
+   cosmo = FlatLambdaCDM(H0, Om0, Ode0, 2.7255, 0, Ob0=Ob0)
+   print('lookback time to z=3')
+   print(cosmo.lookback_time(3))
+   # 11.641584565145552 Gyr
+
+
+def comoving_volume():
+   Om0 = 0.27
+   Ode0 = 0.73
+   Ob0 = 0.044
+   H0 = 70.0
+   cosmo = FlatLambdaCDM(H0, Om0, Ode0, 2.7255, 0, Ob0=Ob0)
+   print('hubble distance')
+   print(cosmo.hubble_distance)
+   # 4282.749400000001 Mpc
+   print('comoving volume to z=3')
+   print(cosmo.comoving_volume(3))
+   # 1179374442443.6943 Mpc3
+
+
+def bench_luminosity_distance():
+    Om0 = 0.27
+    Ode0 = 0.73
+    Ob0 = 0.044
+    H0 = 70.0
+    cosmo = FlatLambdaCDM(H0, Om0, Ode0, 2.7255, 0, Ob0=Ob0)
+    n = 10000
+    print('d_L(z=1): ', timeit.timeit(lambda: cosmo.luminosity_distance(1), number=n)/n, 's')
+    # 15us
+    # Compared to 11us Rust
+    print('d_L(z=2): ', timeit.timeit(lambda: cosmo.luminosity_distance(2), number=n)/n, 's')
+    # 15us
+    # Compared to 22us Rust
+
+
+if __name__=="__main__":
+    # flat_universe_distances_no_relativistic_contribution()
+    # flat_universe_distances_with_radiation_but_no_neutrinos()
+    # lookback_time()
+    # comoving_volume()
+    bench_luminosity_distance()
\ No newline at end of file
diff --git a/benchmarks/requirements.txt b/python_benchmarks/requirements.txt
similarity index 100%
rename from benchmarks/requirements.txt
rename to python_benchmarks/requirements.txt
diff --git a/src/constants.rs b/src/constants.rs
index d3c52bb..2e98477 100644
--- a/src/constants.rs
+++ b/src/constants.rs
@@ -1,3 +1,9 @@
+/// This module contains constants.
+///
+/// Fundamental constants use the [NIST] [CODATA].
+///
+/// [NIST]: <https://physics.nist.gov/cgi-bin/cuu/Value?u|search_for=abbr_in!>
+/// [CODATA 2018]: <https://www.nist.gov/publications/codata-recommended-values-fundamental-physical-constants-2018>
 use once_cell::sync::Lazy;
 
 use crate::{
@@ -6,48 +12,44 @@ use crate::{
         JoulePerKelvin, JoulePerMeter3Kelvin4, JouleSeconds, Meters3PerKgPerSecond2,
         MetersPerSecond, PositiveFloat, WattsPerMeters2Kelvin4,
     },
-    DimensionlessPositiveFloat,
+    DimensionlessPositiveFloat, FloatingPointUnit,
 };
 
+/// PI
 pub const PI: f64 = std::f64::consts::PI;
 
-// Neutrinos
-pub static DEFAULT_NEUTRINO_MASSES: Lazy<[eV; 3]> = Lazy::new(|| {
-    [
-        DimensionlessPositiveFloat::zero(),
-        DimensionlessPositiveFloat::zero(),
-        DimensionlessPositiveFloat::zero(),
-    ]
-});
-pub static DEFAULT_N_EFF: Lazy<DimensionlessPositiveFloat> =
-    Lazy::new(|| DimensionlessPositiveFloat::new(3.04).unwrap()); // WMAP (ApJ Spergel et al 2007)
+/// Speed of light [CODATA 2018]
+pub static C_M_PER_S: MetersPerSecond = 299792458.;
 
-/// Ratio of neutrino to photon temperature
-pub static T_NU_TO_T_GAMMA_RATIO: Lazy<DimensionlessPositiveFloat> =
-    Lazy::new(|| PositiveFloat((4. / 11_f64).powf(1. / 3.)));
+/// Gravitational constant [CODATA 2018]
+pub static G: Meters3PerKgPerSecond2 = 6.67430e-11;
 
-// SI units: https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units
+/// Boltzmann constant [CODATA 2018]
+pub static BOLTZMANN: JoulePerKelvin = 1.380649e-23;
 
-/// Speed of light
-pub static C_M_PER_S: MetersPerSecond = 299792458.;
+/// Stefan-Boltzmann constant [CODATA 2018]
+pub static STEFAN_BOLTZMANN: WattsPerMeters2Kelvin4 = 5.6703744194e-8;
 
-/// Gravitational constant
-pub static G: Meters3PerKgPerSecond2 = 6.6743e-11;
+/// Reduced Planck constant [CODATA 2018]
+pub static H_BAR: JouleSeconds = 1.054571817e-34;
 
-/// Boltzmann constant
-pub static BOLTZMANN: JoulePerKelvin = 1.380649e-23;
+/// Vector of neutrino masses (defaults to 3 massless neutrinos)
+pub static DEFAULT_NEUTRINO_MASSES: Lazy<[eV; 3]> =
+    Lazy::new(|| [eV::zero(), eV::zero(), eV::zero()]);
 
-/// Stefan-Boltzmann constant
-pub static STEFAN_BOLTZMANN: WattsPerMeters2Kelvin4 = 5.670374e-8;
+/// Effective number of neutrinos
+pub static DEFAULT_N_EFF: Lazy<DimensionlessPositiveFloat> =
+    Lazy::new(|| DimensionlessPositiveFloat::new(3.04).unwrap()); // WMAP (ApJ Spergel et al 2007)
 
-/// Reduced Planck constant
-pub static H_BAR: JouleSeconds = 1.05457e-34;
+/// Ratio of neutrino to photon temperature
+pub static T_NU_TO_T_GAMMA_RATIO: Lazy<DimensionlessPositiveFloat> =
+    Lazy::new(|| PositiveFloat((4. / 11_f64).powf(1. / 3.)));
 
 /// $\alpha = \frac{\pi^2 k^4}{15 \hbar^3 c^3}$
 ///
 /// Eqn 2.29 from Ryden
 pub static ALPHA: Lazy<JoulePerMeter3Kelvin4> =
-    Lazy::new(|| PI.powf(2.) * BOLTZMANN.powf(4.) / (15. * H_BAR.powf(3.) * C_M_PER_S.powf(3.)));
+    Lazy::new(|| PI.powi(2) * BOLTZMANN.powi(4) / (15. * H_BAR.powi(3) * C_M_PER_S.powi(3)));
 
 #[cfg(test)]
 mod tests {
diff --git a/src/cosmology.rs b/src/cosmology.rs
index b89b37a..fcdf557 100644
--- a/src/cosmology.rs
+++ b/src/cosmology.rs
@@ -6,10 +6,11 @@ pub use omega_factors::OmegaFactors;
 
 use crate::{
     constants::{self, C_M_PER_S, DEFAULT_NEUTRINO_MASSES, DEFAULT_N_EFF},
-    eV, units,
-    units::{PositiveFloat, Seconds},
-    DimensionlessFloat, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Kelvin, KmPerSecPerMpc,
-    Mpc, Redshift,
+    eV,
+    units::length::{KILOMETER_TO_METER, MPC_TO_KILOMETERS},
+    units::{HInvMpc, PositiveFloat},
+    DimensionlessFloat, DimensionlessPositiveFloat, FloatingPointUnit, Gyr, Kelvin,
+    KilogramsPerMeter3, KmPerSecPerMpc, Meter, Mpc, Redshift, Seconds,
 };
 
 /// Represents an FLRW cosmology.
@@ -20,11 +21,10 @@ use crate::{
 /// # Examples
 ///
 /// ```
-/// use cosmocalc::{Distances, FLRWCosmology};
+/// use cosmocalc::{Distances, Redshift, Mpc, FLRWCosmology, FloatingPointUnit};
 ///
 /// let cosmology = FLRWCosmology::two_component(0.286, 0.714, 69.6);
-/// assert!(cosmology.radial_comoving_distance(2.0).0 > 5273.);
-/// assert!(cosmology.radial_comoving_distance(2.0).0 < 5274.);
+/// assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) > Mpc::new(5273.));
 /// ```
 pub struct FLRWCosmology {
     /// A descriptive name.
@@ -35,8 +35,16 @@ pub struct FLRWCosmology {
     /// Hubble constant at `z=0` (km/(s/Mpc)).
     pub H_0: KmPerSecPerMpc,
 
-    /// Omega factors for this cosmology.
+    /// Base omega factors for this cosmology.
     pub omega: OmegaFactors,
+    /// Omega curvature at z=0
+    pub omega_k0: DimensionlessFloat,
+    /// Omega gamma at z=0
+    pub omega_gamma0: DimensionlessFloat,
+    /// Omega nu at z=0
+    pub omega_nu0: DimensionlessFloat,
+    /// Total omega at z=0
+    pub omega_tot0: DimensionlessFloat,
 
     /// Temperature of the CMB at `z=0`.
     pub T_CMB0: Option<Kelvin>,
@@ -50,15 +58,16 @@ impl FLRWCosmology {
     /// Instantiate a simple two component cosmology.
     pub fn two_component(Omega_M0: f64, Omega_DE0: f64, H_0: f64) -> Self {
         let omega = OmegaFactors::new(Omega_M0, Omega_DE0, Omega_M0).unwrap();
-        Self {
-            name: None,
-            reference: None,
+        Self::new(
+            None,
+            None,
             H_0,
             omega,
-            T_CMB0: Some(PositiveFloat::zero()),
-            N_eff: PositiveFloat::zero(),
-            m_nu: vec![],
-        }
+            Some(0.),
+            Some(PositiveFloat::zero()),
+            Some(vec![]),
+        )
+        .unwrap()
     }
 
     /// Instantiate a new FLRW cosmology.
@@ -71,34 +80,54 @@ impl FLRWCosmology {
         N_eff: Option<DimensionlessPositiveFloat>,
         m_nu: Option<Vec<eV>>,
     ) -> Result<Self, anyhow::Error> {
-        if N_eff.unwrap_or(*DEFAULT_N_EFF).floor()
-            != m_nu
-                .clone()
-                .unwrap_or_else(|| DEFAULT_NEUTRINO_MASSES.to_vec())
-                .len() as f64
-        {
+        let N_eff = N_eff.unwrap_or(*DEFAULT_N_EFF);
+        let m_nu = m_nu.unwrap_or_else(|| DEFAULT_NEUTRINO_MASSES.to_vec());
+
+        if N_eff.floor() != m_nu.len() as f64 {
             return Err(anyhow!(
                 "number of neutrino masses must match the number of effective neutrino species"
             ));
         }
 
+        let critical_density0 = PositiveFloat(
+            3. * H_0.powi(2) / (8. * constants::PI * constants::G * MPC_TO_KILOMETERS.powi(2)),
+        );
+        let omega_gamma0 = match T_CMB0 {
+            Some(T_CMB0) => DimensionlessFloat(
+                *constants::ALPHA * T_CMB0.powi(4) / ((critical_density0).0 * C_M_PER_S.powi(2)),
+            ),
+            None => DimensionlessFloat::zero(),
+        };
+        let omega_nu0 = match T_CMB0 {
+            Some(_) => DimensionlessFloat(
+                7. / 8. * (4.0f64 / 11.).powf(4. / 3.) * N_eff.0 * omega_gamma0.0,
+            ),
+            None => DimensionlessFloat::zero(),
+        };
+        let omega_k0 = omega.curvature_density_0(omega_nu0, omega_gamma0);
+        let omega_tot0 = omega.Omega_M0 + omega_gamma0 + omega_nu0 + omega.Omega_DE0 + omega_k0;
+
         Ok(Self {
             name,
             reference,
             H_0,
             omega,
-            T_CMB0: T_CMB0.map(PositiveFloat),
-            N_eff: N_eff.unwrap_or(*DEFAULT_N_EFF),
-            m_nu: m_nu.unwrap_or_else(|| DEFAULT_NEUTRINO_MASSES.to_vec()),
+            omega_k0,
+            omega_gamma0,
+            omega_nu0,
+            omega_tot0,
+            T_CMB0: T_CMB0.map(Kelvin),
+            N_eff,
+            m_nu,
         })
     }
 
     pub fn E(&self, z: Redshift) -> Mpc {
-        PositiveFloat(
-            (self.omega_m0().0 * (1. + z).powf(3.)
-                + self.omega_k0().0 * (1. + z).powf(2.)
-                + self.omega_de0().0
-                + (self.omega_gamma0().0 + self.omega_nu0().0) * (1. + z).powf(4.))
+        Mpc::new(
+            (self.omega.Omega_M0.0 * (1. + z.0).powi(3)
+                + self.omega_k0.0 * (1. + z.0).powi(2)
+                + self.omega.Omega_DE0.0
+                + (self.omega_gamma0.0 + self.omega_nu0.0) * (1. + z.0).powi(4))
             .sqrt(),
         )
     }
@@ -110,7 +139,7 @@ impl FLRWCosmology {
 
     /// Scale factor at redshift z.
     pub fn scale_factor(&self, z: Redshift) -> DimensionlessPositiveFloat {
-        PositiveFloat(1.0 / (z + 1.0))
+        PositiveFloat(1.0 / (z.0 + 1.0))
     }
 
     /// Dimensionless hubble parameter h where 100 km/s/Mpc * h = H0
@@ -122,32 +151,55 @@ impl FLRWCosmology {
     pub fn hubble_time(&self) -> Seconds {
         // H_0 units are km/s/Mpc so we need to convert Mpc to km
         // such that the distance units cancel
-        PositiveFloat(1. / self.H_0 * units::MPC_TO_KILOMETERS)
+        Seconds::new(1. / self.H_0 * MPC_TO_KILOMETERS)
     }
 
     /// Hubble distance in Mpc: $D_H = c / H_0$.
-    pub fn hubble_distance(&self) -> KmPerSecPerMpc {
+    pub fn hubble_distance(&self) -> Mpc {
         // Factor of 1000 to convert c in m/s to c in km/s so that
         // the units cancel.
-        C_M_PER_S / (self.H_0 * units::KILOMETER_TO_METER)
+        Mpc::new(C_M_PER_S / (self.H_0 * KILOMETER_TO_METER))
     }
 
     /// Hubble distance in h^{-1} Mpc.
-    pub fn hubble_distance_little_h(&self) -> HInvKmPerSecPerMpc {
+    pub fn hubble_distance_little_h(&self) -> HInvMpc {
         C_M_PER_S / (1.0e5)
     }
 
+    /// CMB temperature at redshift z.
+    pub fn T_CMB(&self, z: Redshift) -> Kelvin {
+        if z == Redshift::zero() {
+            self.T_CMB0.unwrap_or_default()
+        } else {
+            Kelvin::new(self.T_CMB0.unwrap_or_default().0 * (z.0 + 1.))
+        }
+    }
+
+    /// Neutrino temperature at redshift z.
+    pub fn T_nu(&self, z: Redshift) -> Kelvin {
+        let T_nu = match self.T_CMB0 {
+            Some(T_cmb) => Kelvin(T_cmb.0 * (*constants::T_NU_TO_T_GAMMA_RATIO).0),
+            None => Kelvin::zero(),
+        };
+
+        if z == Redshift::zero() {
+            T_nu
+        } else {
+            Kelvin::new(T_nu.0 * (z.0 + 1.))
+        }
+    }
+
     /// Critical mass density at redshift z.
-    pub fn critical_density(&self, z: Redshift) -> units::KilogramsPerMeter3 {
-        if z == 0.0 {
+    pub fn critical_density(&self, z: Redshift) -> KilogramsPerMeter3 {
+        if z == Redshift::zero() {
             PositiveFloat(
-                3. * self.H_0.powf(2.)
-                    / (8. * constants::PI * constants::G * units::MPC_TO_KILOMETERS.powf(2.)),
+                3. * self.H_0.powi(2)
+                    / (8. * constants::PI * constants::G * MPC_TO_KILOMETERS.powi(2)),
             )
         } else {
             PositiveFloat(
-                3. * self.H(z).powf(2.)
-                    / (8. * constants::PI * constants::G * units::MPC_TO_KILOMETERS.powf(2.)),
+                3. * self.H(z).powi(2)
+                    / (8. * constants::PI * constants::G * MPC_TO_KILOMETERS.powi(2)),
             )
         }
     }
@@ -156,33 +208,22 @@ impl FLRWCosmology {
     ///
     /// Eqn. 2.28 from Ryden divided by the critical density at `z=0`
     pub fn omega_gamma0(&self) -> DimensionlessFloat {
-        match self.T_CMB0 {
-            Some(T_CMB0) => DimensionlessFloat(
-                *constants::ALPHA * T_CMB0.powf(4.)
-                    / (self.critical_density(0.).0 * C_M_PER_S.powf(2.)),
-            ),
-            None => DimensionlessFloat::zero(),
-        }
+        self.omega_gamma0
     }
 
     /// Dimensionless photon density (density/critical density) at `z>0`
     pub fn omega_gamma(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_gamma0().0 * (1.0 + z).powf(4.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_gamma0.0 * (1.0 + z.0).powi(4) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless neutrino density (density/critical density) at `z=0`
     pub fn omega_nu0(&self) -> DimensionlessFloat {
-        match self.T_CMB0 {
-            Some(_) => DimensionlessFloat(
-                7. / 8. * (4.0f64 / 11.).powf(4. / 3.) * self.N_eff.0 * self.omega_gamma0().0,
-            ),
-            None => DimensionlessFloat::zero(),
-        }
+        self.omega_nu0
     }
 
     /// Dimensionless neutrino density (density/critical density) at `z>0`
     pub fn omega_nu(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_nu0().0 * (1.0 + z).powf(4.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_nu0().0 * (1.0 + z.0).powi(4) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless dark matter density (density/critical density) at `z=0`
@@ -192,18 +233,17 @@ impl FLRWCosmology {
 
     /// Dimensionless dark matter density (density/critical density) at `z>0`
     pub fn omega_dm(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_dm0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_dm0().0 * (1.0 + z.0).powi(3) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless effective curvature density (density/critical density) at `z=0`
     pub fn omega_k0(&self) -> DimensionlessFloat {
-        self.omega
-            .curvature_density_0(self.omega_nu0(), self.omega_gamma0())
+        self.omega_k0
     }
 
     /// Dimensionless effective curvature density (density/critical density) at `z>0`
     pub fn omega_k(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_k0().0 * (1.0 + z).powf(2.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_k0().0 * (1.0 + z.0).powi(2) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless matter density (density/critical density) at `z=0`
@@ -213,7 +253,7 @@ impl FLRWCosmology {
 
     /// Dimensionless matter density (density/critical density) at `z>0`
     pub fn omega_m(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_m0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_m0().0 * (1.0 + z.0).powi(3) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless baryon density (density/critical density) at `z=0`
@@ -223,7 +263,7 @@ impl FLRWCosmology {
 
     /// Dimensionless baryon density (density/critical density) at `z>0`
     pub fn omega_b(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_b0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_b0().0 * (1.0 + z.0).powi(3) * 1. / self.E(z).0.powi(2))
     }
 
     /// Dimensionless dark energy density (density/critical density) at `z=0`
@@ -233,16 +273,12 @@ impl FLRWCosmology {
 
     /// Dimensionless dark energy density (density/critical density) at `z>0`.
     pub fn omega_de(&self, z: Redshift) -> DimensionlessFloat {
-        DimensionlessFloat(self.omega_de0().0 / self.E(z).0.powf(2.))
+        DimensionlessFloat(self.omega_de0().0 / self.E(z).0.powi(2))
     }
 
     /// Dimensionless total density (density/critical density) at `z=0`.
     pub fn omega_tot0(&self) -> DimensionlessFloat {
-        self.omega_m0()
-            + self.omega_gamma0()
-            + self.omega_nu0()
-            + self.omega_de0()
-            + self.omega_k0()
+        self.omega_tot0
     }
 
     /// Dimensionless total density (density/critical density) at `z>0`.
@@ -260,11 +296,26 @@ impl FLRWCosmology {
             && self.omega_tot0() == DimensionlessFloat::one()
     }
 
-    /// Neutrino temperature at `z=0`.
-    pub fn neutrino_temperature0(&self) -> Kelvin {
-        match self.T_CMB0 {
-            Some(T_cmb) => PositiveFloat(T_cmb.0 * (*constants::T_NU_TO_T_GAMMA_RATIO).0),
-            None => DimensionlessPositiveFloat::zero(),
+    /// Lookback time
+    ///
+    /// The difference in ages of the universe from now to when the light
+    /// was emitted from the object at `z`.
+    pub fn lookback_time(&self, z: Redshift) -> Gyr {
+        let mut integrand: f64 = 0.0;
+        let mut z_prime = 0.0;
+        let DZ = 0.0001;
+        while z_prime < z.0 {
+            z_prime += DZ / 2.;
+            integrand += (DZ / 2.) / ((1. + z_prime) * self.E(Redshift::new(z_prime)).0);
         }
+        Seconds::new(self.hubble_time().0 * integrand).into()
+    }
+
+    /// Lookback distance
+    ///
+    /// Proper distance between now and redshift z
+    pub fn lookback_distance(&self, z: Redshift) -> Mpc {
+        let lookback_time_seconds: Seconds = self.lookback_time(z).into();
+        Meter::new(lookback_time_seconds.0 * constants::C_M_PER_S).into()
     }
 }
diff --git a/src/cosmology/omega_factors.rs b/src/cosmology/omega_factors.rs
index 79cc44e..67c6a82 100644
--- a/src/cosmology/omega_factors.rs
+++ b/src/cosmology/omega_factors.rs
@@ -1,4 +1,4 @@
-use crate::DimensionlessFloat;
+use crate::{units::FloatingPointUnit, DimensionlessFloat};
 
 /// Represents a collection of dimensionless density parameters.
 pub struct OmegaFactors {
@@ -17,9 +17,9 @@ impl OmegaFactors {
         }
 
         Ok(OmegaFactors {
-            Omega_M0: DimensionlessFloat::new(Omega_M0)?,
-            Omega_DE0: DimensionlessFloat::new(Omega_DE0)?,
-            Omega_b0: DimensionlessFloat::new(Omega_b0)?,
+            Omega_M0: DimensionlessFloat::new(Omega_M0),
+            Omega_DE0: DimensionlessFloat::new(Omega_DE0),
+            Omega_b0: DimensionlessFloat::new(Omega_b0),
         })
     }
 
diff --git a/src/distances.rs b/src/distances.rs
index c246083..379698d 100644
--- a/src/distances.rs
+++ b/src/distances.rs
@@ -1,10 +1,14 @@
-use crate::{units::PositiveFloat, DimensionlessFloat, FLRWCosmology, Mpc, Redshift};
+use crate::{
+    constants,
+    units::{FloatingPointUnit, Mpc3},
+    DimensionlessFloat, FLRWCosmology, Mpc, Redshift,
+};
 
 /// Bin width in redshift integrals.
 const DZ: f64 = 0.0001;
 
-/// Cosmological distances following Hogg 2000
-/// https://arxiv.org/pdf/astro-ph/9905116.pdf
+/// Cosmological distances following [Hogg 2000]
+/// [Hogg 2000]: <https://arxiv.org/pdf/astro-ph/9905116.pdf>
 pub trait Distances {
     /// Line of sight (radial) comoving distance in Megaparsecs.
     fn radial_comoving_distance(&self, z: Redshift) -> Mpc;
@@ -13,18 +17,33 @@ pub trait Distances {
     /// Angular diameter distance in Megaparsecs.
     fn angular_diameter_distance(&self, z: Redshift) -> Mpc;
     /// Luminosity distance in Megaparsecs.
+    ///
+    /// This should be used with bolometric quantities, i.e.
+    /// it does not include K-corrections.
     fn luminosity_distance(&self, z: Redshift) -> Mpc;
+    /// Comoving volume.
+    fn comoving_volume(&self, z: Redshift) -> Mpc3;
 }
 
 impl Distances for FLRWCosmology {
     fn radial_comoving_distance(&self, z: Redshift) -> Mpc {
-        let mut integrand: f64 = 0.0;
-        let mut z_prime = 0.0;
-        while z_prime < z {
-            z_prime += DZ / 2.;
-            integrand += (DZ / 2.) / self.E(z_prime).0;
-        }
-        PositiveFloat(self.hubble_distance() * integrand)
+        // TODO: To speed up further pick x equal size bins
+        // We operate over 1e4
+        let max_range = (z.0 as i64) * 10000;
+        let step = DZ;
+        // Quantity inside the sqrt is E()
+        let integrand: f64 = (0..max_range)
+            .map(|z_prime_e4| {
+                step / (self.omega.Omega_M0.0 * (1. + z_prime_e4 as f64 / 10000.).powi(3)
+                    + self.omega_k0.0 * (1. + z_prime_e4 as f64 / 10000.).powi(2)
+                    + self.omega.Omega_DE0.0
+                    + (self.omega_gamma0.0 + self.omega_nu0.0)
+                        * (1. + z_prime_e4 as f64 / 10000.).powi(4))
+                .sqrt()
+            })
+            .sum();
+
+        Mpc::new(self.hubble_distance().0 * integrand)
     }
 
     fn transverse_comoving_distance(&self, z: Redshift) -> Mpc {
@@ -33,9 +52,9 @@ impl Distances for FLRWCosmology {
         if omega_k > DimensionlessFloat::zero() {
             // Negative curvature (open)
             let sqrt_omega_k = (omega_k.0).sqrt();
-            PositiveFloat(
-                self.hubble_distance() * 1. / sqrt_omega_k
-                    * f64::sinh(sqrt_omega_k * radial_comoving.0 / self.hubble_distance()),
+            Mpc::new(
+                self.hubble_distance().0 * 1. / sqrt_omega_k
+                    * f64::sinh(sqrt_omega_k * radial_comoving.0 / self.hubble_distance().0),
             )
         } else if omega_k == DimensionlessFloat::zero() {
             // Flat
@@ -43,26 +62,60 @@ impl Distances for FLRWCosmology {
         } else {
             // Positive curvature (closed)
             let abs_sqrt_omega_k = (-1. * omega_k.0).sqrt();
-            PositiveFloat(
-                self.hubble_distance() * 1. / abs_sqrt_omega_k
-                    * f64::sin(abs_sqrt_omega_k * radial_comoving.0 / self.hubble_distance()),
+            Mpc::new(
+                self.hubble_distance().0 * 1. / abs_sqrt_omega_k
+                    * f64::sin(abs_sqrt_omega_k * radial_comoving.0 / self.hubble_distance().0),
             )
         }
     }
 
     fn angular_diameter_distance(&self, z: Redshift) -> Mpc {
-        PositiveFloat(self.transverse_comoving_distance(z).0 / (1. + z))
+        Mpc::new(self.transverse_comoving_distance(z).0 / (1. + z.0))
     }
 
     fn luminosity_distance(&self, z: Redshift) -> Mpc {
         // TODO: K-CORRECTIONS
-        PositiveFloat(self.transverse_comoving_distance(z).0 * (1. + z))
+        Mpc::new(self.transverse_comoving_distance(z).0 * (1. + z.0))
+    }
+
+    /// Comoving volume
+    fn comoving_volume(&self, z: Redshift) -> Mpc3 {
+        // KmPerSecPerMpc
+        let d_H = self.hubble_distance().0;
+        let d_H_cubed = d_H.powi(3);
+
+        let transverse_comoving = self.transverse_comoving_distance(z);
+        let omega_k = self.omega_k(z);
+        if omega_k > DimensionlessFloat::zero() {
+            // Negative curvature (open)
+            let sqrt_omega_k = (omega_k.0).sqrt();
+            let coefficient = 4. * constants::PI * d_H_cubed / (2. * omega_k.0);
+            let term_1_in_parens = transverse_comoving.0 / d_H
+                * (1. + omega_k.0 * transverse_comoving.powi(2) / d_H.powi(2)).sqrt();
+            let term_2_in_parens =
+                1. / sqrt_omega_k * f64::asinh(sqrt_omega_k * transverse_comoving.0 / d_H);
+
+            coefficient * (term_1_in_parens - term_2_in_parens)
+        } else if omega_k == DimensionlessFloat::zero() {
+            // Flat
+            4. * constants::PI * transverse_comoving.powi(3) / 3.
+        } else {
+            // Positive curvature (closed)
+            let sqrt_omega_k = (-1. * omega_k.0).sqrt();
+            let coefficient = 4. * constants::PI * d_H_cubed / (2. * omega_k.0);
+            let term_1_in_parens = transverse_comoving.0 / d_H
+                * (1. + omega_k.0 * transverse_comoving.powi(2) / d_H.powi(2)).sqrt();
+            let term_2_in_parens =
+                1. / sqrt_omega_k * f64::asin(sqrt_omega_k * transverse_comoving.0 / d_H);
+
+            coefficient * (term_1_in_parens - term_2_in_parens)
+        }
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use crate::{cosmology::OmegaFactors, DimensionlessPositiveFloat};
+    use crate::{cosmology::OmegaFactors, eV, units::PositiveFloat};
 
     use super::*;
 
@@ -72,13 +125,12 @@ mod tests {
         let omegas = OmegaFactors::new(0.286, 0.714, 0.05).unwrap();
         let cosmology = FLRWCosmology::new(None, None, 69.6, omegas, None, None, None).unwrap();
 
-        // Megaparsecs
-        assert!(cosmology.radial_comoving_distance(3.0).0 > 6482.5);
-        assert!(cosmology.radial_comoving_distance(3.0).0 < 6482.6);
-        assert!(cosmology.angular_diameter_distance(3.0).0 > 1620.6);
-        assert!(cosmology.angular_diameter_distance(3.0).0 < 1620.7);
-        assert!(cosmology.luminosity_distance(3.0).0 > 25930.0);
-        assert!(cosmology.luminosity_distance(3.0).0 < 25930.2);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) > Mpc::new(6482.5));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) < Mpc::new(6482.8));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) > Mpc::new(1620.6));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) < Mpc::new(1620.7));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) > Mpc::new(25930.0));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) < Mpc::new(25931.0));
     }
 
     #[test]
@@ -96,13 +148,12 @@ mod tests {
         )
         .unwrap();
 
-        // Megaparsecs
-        assert!(cosmology.radial_comoving_distance(3.0).0 > 6395.0);
-        assert!(cosmology.radial_comoving_distance(3.0).0 < 6399.0);
-        assert!(cosmology.angular_diameter_distance(3.0).0 > 1599.0);
-        assert!(cosmology.angular_diameter_distance(3.0).0 < 1600.0);
-        assert!(cosmology.luminosity_distance(3.0).0 > 25588.);
-        assert!(cosmology.luminosity_distance(3.0).0 < 25589.);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) > Mpc::new(6395.0));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) < Mpc::new(6399.0));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) > Mpc::new(1599.0));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) < Mpc::new(1600.0));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) > Mpc::new(25589.));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) < Mpc::new(25594.));
     }
 
     #[test]
@@ -115,21 +166,16 @@ mod tests {
             omegas,
             Some(2.7255),
             Some(PositiveFloat(3.04)),
-            Some(vec![
-                DimensionlessPositiveFloat::zero(),
-                DimensionlessPositiveFloat::zero(),
-                DimensionlessPositiveFloat::zero(),
-            ]),
+            Some(vec![eV::zero(), eV::zero(), eV::zero()]),
         )
         .unwrap();
 
-        // Megaparsecs
-        assert!(cosmology.radial_comoving_distance(3.0).0 > 6598.);
-        assert!(cosmology.radial_comoving_distance(3.0).0 < 6598.5);
-        assert!(cosmology.angular_diameter_distance(3.0).0 > 1600.5);
-        assert!(cosmology.angular_diameter_distance(3.0).0 < 1700.0);
-        assert!(cosmology.luminosity_distance(3.0).0 > 25000.);
-        assert!(cosmology.luminosity_distance(3.0).0 < 27000.);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) > Mpc::new(6598.));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) < Mpc::new(6598.5));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) > Mpc::new(1600.5));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) < Mpc::new(1700.0));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) > Mpc::new(25000.));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) < Mpc::new(27000.));
     }
 
     #[test]
@@ -138,13 +184,13 @@ mod tests {
         let cosmology = FLRWCosmology::new(None, None, 69.6, omegas, None, None, None).unwrap();
 
         // Megaparsecs
-        assert!(cosmology.radial_comoving_distance(3.0).0 > 5200.);
-        assert!(cosmology.radial_comoving_distance(3.0).0 < 5300.);
-        assert!(cosmology.angular_diameter_distance(3.0).0 > 1250.);
-        assert!(cosmology.angular_diameter_distance(3.0).0 < 1600.);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) > Mpc::new(5200.));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(3.0)) < Mpc::new(5300.));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) > Mpc::new(1250.));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(3.0)) < Mpc::new(1600.));
         // No k-corrections here
-        assert!(cosmology.luminosity_distance(3.0).0 > 22000.);
-        assert!(cosmology.luminosity_distance(3.0).0 < 24000.);
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) > Mpc::new(22000.));
+        assert!(cosmology.luminosity_distance(Redshift::new(3.0)) < Mpc::new(24000.));
     }
 
     #[test]
@@ -153,19 +199,40 @@ mod tests {
         let cosmology = FLRWCosmology::new(None, None, 69.6, omegas, None, None, None).unwrap();
 
         // Megaparsecs
-        assert!(cosmology.radial_comoving_distance(2.0).0 > 5000.);
-        assert!(cosmology.radial_comoving_distance(2.0).0 < 6000.);
-        assert!(cosmology.angular_diameter_distance(2.0).0 > 1500.);
-        assert!(cosmology.angular_diameter_distance(2.0).0 < 2000.);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) > Mpc::new(5000.));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) < Mpc::new(6000.));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(2.0)) > Mpc::new(1500.));
+        assert!(cosmology.angular_diameter_distance(Redshift::new(2.0)) < Mpc::new(2000.));
         // No k-corrections here
-        assert!(cosmology.luminosity_distance(2.0).0 > 14000.);
-        assert!(cosmology.luminosity_distance(2.0).0 < 16000.);
+        assert!(cosmology.luminosity_distance(Redshift::new(2.0)) > Mpc::new(14000.));
+        assert!(cosmology.luminosity_distance(Redshift::new(2.0)) < Mpc::new(16000.));
     }
 
     #[test]
     fn simple_two_component() {
         let cosmology = FLRWCosmology::two_component(0.286, 0.714, 69.6);
-        assert!(cosmology.radial_comoving_distance(2.0).0 > 5273.);
-        assert!(cosmology.radial_comoving_distance(2.0).0 < 5274.);
+        assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) > Mpc::new(5273.));
+        assert!(cosmology.radial_comoving_distance(Redshift::new(2.0)) < Mpc::new(5274.));
+    }
+
+    /// This test is here for profilng purposes. It runs the luminosity distance computation
+    /// for many iterations such that profiling data can be collected.
+    #[ignore]
+    #[test]
+    fn simple_luminosity() {
+        let omegas = OmegaFactors::new(0.27, 0.73, 0.044).unwrap();
+        let cosmology = FLRWCosmology::new(None, None, 70.0, omegas, None, None, None).unwrap();
+        for _ in 0..10000000 {
+            cosmology.luminosity_distance(Redshift::new(2.0));
+        }
+    }
+
+    #[test]
+    fn comoving_volume() {
+        // TESTED vs: astro.py 5.1 FlatLambdaCDM. Within 10e8 Mpc3.
+        let omegas = OmegaFactors::new(0.27, 0.73, 0.044).unwrap();
+        let cosmology = FLRWCosmology::new(None, None, 70.0, omegas, None, None, None).unwrap();
+        assert!(cosmology.comoving_volume(Redshift::new(3.0)) > 1179361698730.);
+        assert!(cosmology.comoving_volume(Redshift::new(3.0)) < 1179470000000.);
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index c22cf3e..c608edc 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,11 +1,11 @@
 #![allow(non_snake_case)]
+#![allow(non_camel_case_types)]
 pub mod constants;
 pub mod cosmology;
 pub mod dark_energy;
 pub mod distances;
 pub mod redshift;
 pub mod units;
-pub mod utils;
 
 pub use cosmology::FLRWCosmology;
 pub use distances::Distances;
@@ -13,6 +13,14 @@ pub use distances::Distances;
 // Common units are re-exported from the crate root for convenience.
 pub use redshift::Redshift;
 pub use units::{
-    eV, DimensionlessFloat, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Joule, Kelvin,
-    Kilogram, KmPerSecPerMpc, Mpc,
+    energy::{eV, Joule},
+    length::{Kilometer, Meter, Mpc},
+    mass::{Gram, Kilogram},
+    temperature::Kelvin,
+    time::{Gyr, Seconds},
+    DimensionlessFloat, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, KilogramsPerMeter3,
+    KmPerSecPerMpc,
 };
+
+// Common traits are re-exported from the crate root also.
+pub use crate::units::traits::FloatingPointUnit;
diff --git a/src/redshift.rs b/src/redshift.rs
index 3be221c..9bd42f3 100644
--- a/src/redshift.rs
+++ b/src/redshift.rs
@@ -1 +1,5 @@
-pub type Redshift = f64;
+use std::ops::{Add, Sub};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { Redshift }
diff --git a/src/units.rs b/src/units.rs
index dab463d..993523d 100644
--- a/src/units.rs
+++ b/src/units.rs
@@ -1,33 +1,24 @@
 use anyhow::anyhow;
 
+pub mod dimensionless;
+pub mod energy;
+pub mod length;
+pub(crate) mod macros;
+pub mod mass;
+pub mod temperature;
+pub mod time;
+pub mod traits;
+
+pub use dimensionless::DimensionlessFloat;
+pub use traits::FloatingPointUnit;
+
 // Continuous positive quantities that are dimensionless (e.g. ratios like the omegas)
 pub type DimensionlessPositiveFloat = PositiveFloat;
 
-// Time
-pub type Gyr = PositiveFloat;
-pub type Seconds = PositiveFloat;
-
-// Energy
-#[allow(non_camel_case_types)]
-pub type eV = PositiveFloat;
-pub type Joule = PositiveFloat;
-
-// Mass
-pub type Kilogram = PositiveFloat;
-pub type Gram = PositiveFloat;
-
-// Temperatures
-pub type Kelvin = PositiveFloat;
-
 // Hubble parameter units
 pub type KmPerSecPerMpc = f64;
 pub type HInvKmPerSecPerMpc = f64;
 
-// Distances
-pub type Meter = PositiveFloat;
-pub type Km = PositiveFloat;
-pub type Mpc = PositiveFloat;
-
 // Densities
 pub type KilogramsPerMeter3 = PositiveFloat;
 
@@ -40,54 +31,8 @@ pub type JouleSeconds = f64;
 pub type JoulePerMeter3Kelvin4 = f64;
 pub type WattsPerMeters2Kelvin4 = f64;
 pub type JoulePerKelvin = f64;
-
-// Conversions
-pub const KILOMETER_TO_METER: f64 = 1000.;
-pub const MPC_TO_METERS: f64 = 3.086e+22;
-pub const MPC_TO_KILOMETERS: f64 = 3.086e+19;
-
-/// Represents continuous dimensionless physical quantities.
-#[derive(Clone, Copy, Debug, PartialOrd, PartialEq)]
-pub struct DimensionlessFloat(pub f64);
-
-impl DimensionlessFloat {
-    pub fn new(x: f64) -> Result<Self, anyhow::Error> {
-        Ok(Self(x))
-    }
-
-    pub fn zero() -> Self {
-        Self(0.)
-    }
-
-    pub fn one() -> Self {
-        Self(1.)
-    }
-
-    // Passthrough methods for convenience
-    pub fn floor(&self) -> f64 {
-        self.0.floor()
-    }
-
-    pub fn powf(&self, exp: f64) -> f64 {
-        self.0.powf(exp)
-    }
-}
-
-impl std::ops::Sub for DimensionlessFloat {
-    type Output = DimensionlessFloat;
-
-    fn sub(self, rhs: Self) -> Self::Output {
-        DimensionlessFloat(self.0 - rhs.0)
-    }
-}
-
-impl std::ops::Add for DimensionlessFloat {
-    type Output = DimensionlessFloat;
-
-    fn add(self, rhs: Self) -> Self::Output {
-        DimensionlessFloat(self.0 + rhs.0)
-    }
-}
+pub type Mpc3 = f64;
+pub type HInvMpc = f64;
 
 /// Represents continuous physical quantities that _cannot_ be negative.
 #[derive(Clone, Copy, Debug, PartialOrd, PartialEq)]
diff --git a/src/units/dimensionless.rs b/src/units/dimensionless.rs
new file mode 100644
index 0000000..677e6dc
--- /dev/null
+++ b/src/units/dimensionless.rs
@@ -0,0 +1,5 @@
+use std::ops::{Add, Sub};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { DimensionlessFloat }
diff --git a/src/units/energy.rs b/src/units/energy.rs
new file mode 100644
index 0000000..6e35b92
--- /dev/null
+++ b/src/units/energy.rs
@@ -0,0 +1,6 @@
+use std::ops::{Add, Sub};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { eV }
+floating_point_unit_impl! { Joule }
diff --git a/src/units/length.rs b/src/units/length.rs
new file mode 100644
index 0000000..5ba2e4b
--- /dev/null
+++ b/src/units/length.rs
@@ -0,0 +1,36 @@
+use std::ops::{Add, Sub};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { Meter }
+floating_point_unit_impl! { Kilometer }
+floating_point_unit_impl! { Mpc }
+
+// Conversions
+pub const KILOMETER_TO_METER: f64 = 1000.;
+pub const MPC_TO_METERS: f64 = 3.086e+22;
+pub const MPC_TO_KILOMETERS: f64 = 3.086e+19;
+
+impl From<Kilometer> for Meter {
+    fn from(km: Kilometer) -> Meter {
+        Meter(1000. * km.0)
+    }
+}
+
+impl From<Mpc> for Meter {
+    fn from(mpc: Mpc) -> Self {
+        Meter(mpc.0 * MPC_TO_METERS)
+    }
+}
+
+impl From<Mpc> for Kilometer {
+    fn from(mpc: Mpc) -> Self {
+        Kilometer(mpc.0 * MPC_TO_KILOMETERS)
+    }
+}
+
+impl From<Meter> for Mpc {
+    fn from(meter: Meter) -> Self {
+        Mpc(meter.0 / MPC_TO_METERS)
+    }
+}
diff --git a/src/units/macros.rs b/src/units/macros.rs
new file mode 100644
index 0000000..99a4fa8
--- /dev/null
+++ b/src/units/macros.rs
@@ -0,0 +1,38 @@
+/// This macro exists to implement the basic methods that should
+/// exist on unit/dimension structs that behave like floating point units.
+///
+/// You can add and subtract two values of the same unit.
+macro_rules! floating_point_unit_impl {
+    ($outer : ident) => {
+        #[derive(Clone, Copy, Debug, PartialOrd, PartialEq)]
+        pub struct $outer(pub f64);
+
+        impl FloatingPointUnit for $outer {
+            fn new(x: f64) -> Self {
+                Self(x)
+            }
+
+            fn inner(&self) -> f64 {
+                self.0
+            }
+        }
+
+        impl Add<$outer> for $outer {
+            type Output = $outer;
+
+            fn add(self, b: $outer) -> $outer {
+                $outer((self.0.add(&b.0)))
+            }
+        }
+
+        impl Sub<$outer> for $outer {
+            type Output = $outer;
+
+            fn sub(self, b: $outer) -> $outer {
+                $outer((self.0.sub(&b.0)))
+            }
+        }
+    };
+}
+
+pub(crate) use floating_point_unit_impl;
diff --git a/src/units/mass.rs b/src/units/mass.rs
new file mode 100644
index 0000000..a991f16
--- /dev/null
+++ b/src/units/mass.rs
@@ -0,0 +1,41 @@
+use std::ops::{Add, Sub};
+
+use crate::{
+    constants,
+    units::{macros::floating_point_unit_impl, traits::FloatingPointUnit},
+    Joule,
+};
+
+floating_point_unit_impl! { Kilogram }
+floating_point_unit_impl! { Gram }
+
+impl From<Joule> for Kilogram {
+    fn from(energy: Joule) -> Self {
+        // Convert between energy and mass using $E=mc^2$
+        Kilogram::new(energy.0 / (constants::C_M_PER_S).powi(2))
+    }
+}
+
+impl From<Kilogram> for Joule {
+    fn from(mass: Kilogram) -> Self {
+        // Convert between energy and mass using $E=mc^2$
+        Joule::new(mass.0 * (constants::C_M_PER_S).powi(2))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn mass_energy_equivalance() {
+        // ~1kg
+        let equivalent_mass: Kilogram = Joule::new(8.9e16).into();
+        assert!(equivalent_mass > Kilogram::new(0.99));
+        assert!(equivalent_mass < Kilogram::new(1.01));
+
+        let equivalent_energy: Joule = Kilogram::new(1.).into();
+        assert!(equivalent_energy > Joule::new(8.9e16));
+        assert!(equivalent_energy < Joule::new(9.05e16));
+    }
+}
diff --git a/src/units/temperature.rs b/src/units/temperature.rs
new file mode 100644
index 0000000..590487a
--- /dev/null
+++ b/src/units/temperature.rs
@@ -0,0 +1,14 @@
+use std::{
+    default::Default,
+    ops::{Add, Sub},
+};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { Kelvin }
+
+impl Default for Kelvin {
+    fn default() -> Self {
+        Self(0.)
+    }
+}
diff --git a/src/units/time.rs b/src/units/time.rs
new file mode 100644
index 0000000..2e559d8
--- /dev/null
+++ b/src/units/time.rs
@@ -0,0 +1,21 @@
+use std::ops::{Add, Sub};
+
+use crate::units::{macros::floating_point_unit_impl, traits::FloatingPointUnit};
+
+floating_point_unit_impl! { Gyr }
+floating_point_unit_impl! { Seconds }
+
+pub const SECONDS_PER_YR: f64 = 3.154e+7;
+pub const SECONDS_PER_GYR: f64 = 3.154e+16;
+
+impl From<Seconds> for Gyr {
+    fn from(seconds: Seconds) -> Self {
+        Gyr::new(seconds.0 / SECONDS_PER_GYR)
+    }
+}
+
+impl From<Gyr> for Seconds {
+    fn from(gyrs: Gyr) -> Self {
+        Seconds::new(gyrs.0 * SECONDS_PER_GYR)
+    }
+}
diff --git a/src/units/traits.rs b/src/units/traits.rs
new file mode 100644
index 0000000..cdd820d
--- /dev/null
+++ b/src/units/traits.rs
@@ -0,0 +1,40 @@
+pub trait FloatingPointUnit {
+    /// Create the value.
+    fn new(inner: f64) -> Self;
+
+    /// Get the inner unit.
+    fn inner(&self) -> f64;
+
+    /// Default implementations
+
+    /// Get the zero value for this unit.
+    fn zero() -> Self
+    where
+        Self: std::marker::Sized,
+    {
+        Self::new(0.)
+    }
+
+    /// Get the one value for this unit.
+    fn one() -> Self
+    where
+        Self: std::marker::Sized,
+    {
+        Self::new(1.)
+    }
+
+    /// Round the value down.
+    fn floor(&self) -> f64 {
+        self.inner().floor()
+    }
+
+    /// The inner part of this value raised to a fp power.
+    fn powf(&self, exp: f64) -> f64 {
+        self.inner().powf(exp)
+    }
+
+    /// The inner part of this value raised to an integer.
+    fn powi(&self, exp: i32) -> f64 {
+        self.inner().powi(exp)
+    }
+}
diff --git a/src/utils.rs b/src/utils.rs
deleted file mode 100644
index 81037c8..0000000
--- a/src/utils.rs
+++ /dev/null
@@ -1,26 +0,0 @@
-use crate::{constants, units::PositiveFloat, Joule, Kilogram};
-
-/// Convert energy to mass using $E=mc^2$
-pub fn energy_to_mass(energy: Joule) -> Kilogram {
-    PositiveFloat(energy.0 / (constants::C_M_PER_S).powf(2.))
-}
-
-/// Convert mass to energy using $E=mc^2$
-pub fn mass_to_energy(mass: Kilogram) -> Joule {
-    PositiveFloat(mass.0 * (constants::C_M_PER_S).powf(2.))
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn mass_energy_equivalance() {
-        // ~1kg
-        assert!(energy_to_mass(PositiveFloat(8.9e16)).0 > 0.99);
-        assert!(energy_to_mass(PositiveFloat(8.9e16)).0 < 1.01);
-
-        assert!(mass_to_energy(PositiveFloat(1.)).0 > 8.9e16);
-        assert!(mass_to_energy(PositiveFloat(1.)).0 < 9.05e16);
-    }
-}
diff --git a/tests/test_cosmology.rs b/tests/test_cosmology.rs
index c8e701b..f3ca2fa 100644
--- a/tests/test_cosmology.rs
+++ b/tests/test_cosmology.rs
@@ -1,6 +1,6 @@
 use cosmocalc::{
     cosmology::{FLRWCosmology, OmegaFactors},
-    DimensionlessPositiveFloat,
+    DimensionlessPositiveFloat, FloatingPointUnit, Gyr, Mpc, Redshift,
 };
 
 #[test]
@@ -21,13 +21,13 @@ fn hubble_distance_and_time() {
     assert!(cosmology.hubble_distance_little_h() < 3000.0);
 
     // For H_0 = 70, should be 4285.7 Mpc
-    assert!(cosmology.hubble_distance() > 4000.0);
-    assert!(cosmology.hubble_distance() < 4500.0);
+    assert!(cosmology.hubble_distance() > Mpc::new(4000.0));
+    assert!(cosmology.hubble_distance() < Mpc::new(4500.0));
 
     // D_H in units of h^{-1} Mpc should be equal to D_H in units of Mpc
     assert!(
         cosmology.hubble_distance_little_h()
-            - (cosmology.hubble_distance() * cosmology.little_h().0)
+            - (cosmology.hubble_distance().0 * cosmology.little_h().0)
             < 0.01
     );
 
@@ -41,6 +41,17 @@ fn densities() {
     let omegas = OmegaFactors::new(0.27, 0.73, 0.044).unwrap();
     let cosmology = FLRWCosmology::new(None, None, 70.0, omegas, None, None, None).unwrap();
 
-    assert!(cosmology.critical_density(0.).0 > 8.7e-27);
-    assert!(cosmology.critical_density(0.).0 < 9.5e-27);
+    assert!(cosmology.critical_density(Redshift::zero()).0 > 8.7e-27);
+    assert!(cosmology.critical_density(Redshift::zero()).0 < 9.5e-27);
+}
+
+#[test]
+fn lookback_time() {
+    // TESTED vs: astro.py 5.1 FlatLambdaCDM
+    // Accurate to 0.1 Gyr
+    let omegas = OmegaFactors::new(0.27, 0.73, 0.044).unwrap();
+    let cosmology = FLRWCosmology::new(None, None, 70.0, omegas, None, None, None).unwrap();
+    assert!(cosmology.lookback_time(Redshift::zero()) == Gyr::zero());
+    assert!(cosmology.lookback_time(Redshift::new(3.0)) > Gyr::new(11.64));
+    assert!(cosmology.lookback_time(Redshift::new(3.0)) < Gyr::new(11.65));
 }