omegas

CHANGELOG.md

@@ -1,3 +1,3 @@
 # 0.1.0
 
-Initial stable release.
+Initial release.

README.md

@@ -18,8 +18,3 @@ This project requires [`Rust`](https://www.rust-lang.org/tools/install). Once in
 cargo build
 cargo test
 ```
-
-## Guidelines
-
-* Newtypes should be used for all units to provide compile-time checking of 
-dimensions.

src/constants.rs

@@ -2,7 +2,10 @@ use once_cell::sync::Lazy;
 
 use crate::{
     eV,
-    units::{Meters3PerKgPerSecond2, MetersPerSecond, PositiveFloat},
+    units::{
+        JoulePerKelvin, JoulePerMeter3Kelvin4, JouleSeconds, Meters3PerKgPerSecond2,
+        MetersPerSecond, PositiveFloat, WattsPerMeters2Kelvin4,
+    },
     DimensionlessPositiveFloat,
 };
 
@@ -23,7 +26,35 @@ pub static T_NU_TO_T_GAMMA_RATIO: Lazy<DimensionlessPositiveFloat> =
     Lazy::new(|| PositiveFloat((4. / 11_f64).powf(1. / 3.)));
 
 // SI units: https://en.wikipedia.org/wiki/2019_redefinition_of_the_SI_base_units
+
+/// Speed of light
 pub static C_M_PER_S: MetersPerSecond = 299792458.;
 
 /// Gravitational constant
 pub static G: Meters3PerKgPerSecond2 = 6.6743e-11;
+
+/// Boltzmann constant
+pub static BOLTZMANN: JoulePerKelvin = 1.380649e-23;
+
+/// Stefan-Boltzmann constant
+pub static STEFAN_BOLTZMANN: WattsPerMeters2Kelvin4 = 5.670374e-8;
+
+/// Reduced Planck constant
+pub static H_BAR: JouleSeconds = 1.05457e-34;
+
+/// $\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.)));
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn alpha() {
+        assert!(*ALPHA > 7.0e-16);
+        assert!(*ALPHA < 8.0e-16);
+    }
+}

src/cosmology.rs

@@ -28,7 +28,7 @@ pub struct FLRWCosmology {
     pub omega: OmegaFactors,
 
     /// Temperature of the CMB at `z=0`.
-    pub T_CMB0: Kelvin,
+    pub T_CMB0: Option<Kelvin>,
     /// Number of effective neutrino species.
     pub N_eff: DimensionlessPositiveFloat,
     /// Mass of neutrino species in eV.
@@ -42,7 +42,7 @@ impl FLRWCosmology {
         reference: Option<String>,
         H_0: f64,
         omega: OmegaFactors,
-        T_CMB0: Option<Kelvin>,
+        T_CMB0: Option<f64>,
         N_eff: Option<DimensionlessPositiveFloat>,
         m_nu: Option<Vec<eV>>,
     ) -> Result<Self, anyhow::Error> {
@@ -62,7 +62,7 @@ impl FLRWCosmology {
             reference,
             H_0,
             omega,
-            T_CMB0: T_CMB0.unwrap_or(*ZERO),
+            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()),
         })
@@ -72,8 +72,9 @@ impl FLRWCosmology {
         PositiveFloat(
             (self.omega_m0().0 * (1. + z).powf(3.)
                 + self.omega_k0().0 * (1. + z).powf(2.)
-                + self.omega_de0().0)
-                .sqrt(),
+                + self.omega_de0().0
+                + (self.omega_gamma0().0 + self.omega_nu0().0) * (1. + z).powf(4.))
+            .sqrt(),
         )
     }
 
@@ -113,45 +114,103 @@ impl FLRWCosmology {
 
     /// Critical mass density at redshift z.
     pub fn critical_density(&self, z: Redshift) -> units::KilogramsPerMeter3 {
-        PositiveFloat(
-            3. * self.H(z).powf(2.)
-                / (8. * constants::PI * constants::G * units::MPC_TO_KILOMETERS.powf(2.)),
-        )
+        if z == 0.0 {
+            PositiveFloat(
+                3. * self.H_0.powf(2.)
+                    / (8. * constants::PI * constants::G * units::MPC_TO_KILOMETERS.powf(2.)),
+            )
+        } else {
+            PositiveFloat(
+                3. * self.H(z).powf(2.)
+                    / (8. * constants::PI * constants::G * units::MPC_TO_KILOMETERS.powf(2.)),
+            )
+        }
     }
 
     /// Dimensionless photon density (density/critical density) at `z=0`.
+    ///
+    /// Eqn. 2.28 from Ryden divided by the critical density at `z=0`
     pub fn omega_gamma0(&self) -> DimensionlessPositiveFloat {
-        // TODO
-        PositiveFloat::new(0.0).unwrap()
+        match self.T_CMB0 {
+            Some(T_CMB0) => PositiveFloat(
+                *constants::ALPHA * T_CMB0.powf(4.)
+                    / (self.critical_density(0.).0 * C_M_PER_S.powf(2.)),
+            ),
+            None => *ZERO,
+        }
+    }
+
+    /// Dimensionless photon density (density/critical density) at `z>0`
+    pub fn omega_gamma(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_gamma0().0 * (1.0 + z).powf(4.) * 1. / self.E(z).0.powf(2.))
     }
 
     /// Dimensionless neutrino density (density/critical density) at `z=0`
     pub fn omega_nu0(&self) -> DimensionlessPositiveFloat {
-        // TODO
-        PositiveFloat::new(0.0).unwrap()
+        match self.T_CMB0 {
+            Some(_) => PositiveFloat(
+                7. / 8. * (4.0f64 / 11.).powf(4. / 3.) * self.N_eff.0 * self.omega_gamma0().0,
+            ),
+            None => *ZERO,
+        }
+    }
+
+    /// Dimensionless neutrino density (density/critical density) at `z>0`
+    pub fn omega_nu(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_nu0().0 * (1.0 + z).powf(4.) * 1. / self.E(z).0.powf(2.))
     }
 
     /// Dimensionless dark matter density (density/critical density) at `z=0`
     pub fn omega_dm0(&self) -> DimensionlessPositiveFloat {
         self.omega.omega_dark_matter_density_0()
     }
 
+    /// Dimensionless dark matter density (density/critical density) at `z>0`
+    pub fn omega_dm(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_dm0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+    }
+
     /// Dimensionless effective curvature density (density/critical density) at `z=0`
     pub fn omega_k0(&self) -> DimensionlessPositiveFloat {
         self.omega
             .curvature_density_0(self.omega_nu0(), self.omega_gamma0())
     }
 
+    /// Dimensionless effective curvature density (density/critical density) at `z>0`
+    pub fn omega_k(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_k0().0 * (1.0 + z).powf(2.) * 1. / self.E(z).0.powf(2.))
+    }
+
     /// Dimensionless matter density (density/critical density) at `z=0`
     pub fn omega_m0(&self) -> DimensionlessPositiveFloat {
         self.omega.Omega_M0
     }
 
+    /// Dimensionless matter density (density/critical density) at `z>0`
+    pub fn omega_m(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_m0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+    }
+
+    /// Dimensionless baryon density (density/critical density) at `z=0`
+    pub fn omega_b0(&self) -> DimensionlessPositiveFloat {
+        self.omega.Omega_b0
+    }
+
+    /// Dimensionless baryon density (density/critical density) at `z>0`
+    pub fn omega_b(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_b0().0 * (1.0 + z).powf(3.) * 1. / self.E(z).0.powf(2.))
+    }
+
     /// Dimensionless dark energy density (density/critical density) at `z=0`
     pub fn omega_de0(&self) -> DimensionlessPositiveFloat {
         self.omega.Omega_DE0
     }
 
+    /// Dimensionless dark energy density (density/critical density) at `z>0`.
+    pub fn omega_de(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        PositiveFloat(self.omega_de0().0 / self.E(z).0.powf(2.))
+    }
+
     /// Dimensionless total density (density/critical density) at `z=0`.
     pub fn omega_tot0(&self) -> DimensionlessPositiveFloat {
         self.omega_m0()
@@ -161,13 +220,25 @@ impl FLRWCosmology {
             + self.omega_k0()
     }
 
+    /// Dimensionless total density (density/critical density) at `z>0`.
+    pub fn omega_tot(&self, z: Redshift) -> DimensionlessPositiveFloat {
+        self.omega_m(z)
+            + self.omega_gamma(z)
+            + self.omega_nu(z)
+            + self.omega_de(z)
+            + self.omega_k(z)
+    }
+
     /// Whether this cosmology is spatially flat
     pub fn is_flat(&self) -> bool {
         self.omega_k0() == *ZERO && self.omega_tot0() == *ONE
     }
 
     /// Neutrino temperature at `z=0`.
     pub fn neutrino_temperature0(&self) -> Kelvin {
-        PositiveFloat(self.T_CMB0.0 * (*constants::T_NU_TO_T_GAMMA_RATIO).0)
+        match self.T_CMB0 {
+            Some(T_cmb) => PositiveFloat(T_cmb.0 * (*constants::T_NU_TO_T_GAMMA_RATIO).0),
+            None => *ZERO,
+        }
     }
 }

src/cosmology/omega_factors.rs

@@ -34,6 +34,9 @@ impl OmegaFactors {
         omega_nu0: DimensionlessPositiveFloat,
         omega_gamma0: DimensionlessPositiveFloat,
     ) -> DimensionlessPositiveFloat {
+        if self.Omega_M0 + self.Omega_DE0 + omega_nu0 + omega_gamma0 > PositiveFloat(1.0) {
+            unimplemented!();
+        }
         PositiveFloat(1.0) - self.Omega_M0 - self.Omega_DE0 - omega_nu0 - omega_gamma0
     }
 }

src/dark_energy.rs

@@ -4,4 +4,4 @@ pub trait DarkEnergyEquationOfState {
     fn w(&self, z: Vec<Redshift>) -> Vec<f64>;
 }
 
-// TODO: impl DarkEnergyEquationOfState for Cosmology
+// TODO: impl DarkEnergyEquationOfState for FLRWCosmology

src/distances.rs

@@ -31,23 +31,24 @@ impl Distances for FLRWCosmology {
         todo!()
     }
 
-    fn angular_diameter_distance(&self, _z: Redshift) -> Mpc {
-        todo!()
+    fn angular_diameter_distance(&self, z: Redshift) -> Mpc {
+        PositiveFloat(self.transverse_comoving_distance(z).0 / (1. + z))
     }
 
-    fn luminosity_distance(&self, _z: Redshift) -> Mpc {
-        todo!()
+    fn luminosity_distance(&self, z: Redshift) -> Mpc {
+        // TODO: K-CORRECTIONS
+        PositiveFloat(self.transverse_comoving_distance(z).0 * (1. + z))
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use crate::cosmology::OmegaFactors;
+    use crate::{constants::ZERO, cosmology::OmegaFactors};
 
     use super::*;
 
     #[test]
-    fn flat_universe_distances() {
+    fn flat_universe_distances_no_relativistic_contribution() {
         let omegas = OmegaFactors::new(0.286, 0.714, 0.05).unwrap();
         let cosmology = FLRWCosmology::new(None, None, 69.6, omegas, None, None, None).unwrap();
 
@@ -56,4 +57,42 @@ mod tests {
         assert!(cosmology.radial_comoving_distance(3.0).0 > 6482.5);
         assert!(cosmology.radial_comoving_distance(3.0).0 < 6483.0);
     }
+
+    #[test]
+    fn flat_universe_distances_with_radiation_but_no_neutrinos() {
+        let omegas = OmegaFactors::new(0.25, 0.7, 0.05).unwrap();
+        let cosmology = FLRWCosmology::new(
+            None,
+            None,
+            69.6,
+            omegas,
+            Some(2.7255),
+            Some(PositiveFloat(0.)),
+            Some(vec![]),
+        )
+        .unwrap();
+
+        // Megaparsecs
+        assert!(cosmology.radial_comoving_distance(3.0).0 > 6598.5);
+        assert!(cosmology.radial_comoving_distance(3.0).0 < 6599.0);
+    }
+
+    #[test]
+    fn flat_universe_distances_with_radiation_and_neutrinos() {
+        let omegas = OmegaFactors::new(0.25, 0.7, 0.05).unwrap();
+        let cosmology = FLRWCosmology::new(
+            None,
+            None,
+            69.6,
+            omegas,
+            Some(2.7255),
+            Some(PositiveFloat(3.04)),
+            Some(vec![*ZERO, *ZERO, *ZERO]),
+        )
+        .unwrap();
+
+        // Megaparsecs
+        assert!(cosmology.radial_comoving_distance(3.0).0 > 6598.);
+        assert!(cosmology.radial_comoving_distance(3.0).0 < 6598.5);
+    }
 }

src/lib.rs

@@ -5,9 +5,13 @@ pub mod dark_energy;
 pub mod distances;
 pub mod redshift;
 pub mod units;
+pub mod utils;
 
 pub use cosmology::FLRWCosmology;
 
 // Common units are re-exported from the crate root for convenience.
 pub use redshift::Redshift;
-pub use units::{eV, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Kelvin, KmPerSecPerMpc, Mpc};
+pub use units::{
+    eV, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Joule, Kelvin, Kilogram, KmPerSecPerMpc,
+    Mpc,
+};

src/units.rs

@@ -10,6 +10,11 @@ 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;
@@ -30,6 +35,11 @@ pub type KilogramsPerMeter3 = PositiveFloat;
 // TODO: Work out a better way to handle types for composite unit information
 pub type MetersPerSecond = f64;
 pub type Meters3PerKgPerSecond2 = f64;
+pub type Meters2KgPerSecond2Kelvin = f64;
+pub type JouleSeconds = f64;
+pub type JoulePerMeter3Kelvin4 = f64;
+pub type WattsPerMeters2Kelvin4 = f64;
+pub type JoulePerKelvin = f64;
 
 // Conversions
 pub const KILOMETER_TO_METER: f64 = 1000.;

src/utils.rs

@@ -0,0 +1,26 @@
+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);
+    }
+}