0.1.0

Cargo.toml

@@ -8,6 +8,7 @@ keywords = ["cosmology", "physics", "astronomy"]
 categories = ["physics"]
 include = ["Cargo.toml", "src", "README.md"]
 edition = "2021"
+license = "MIT OR Apache-2.0"
 
 [dependencies]
 anyhow = "1"

README.md

@@ -1,12 +1,24 @@
 # cosmocalc.rs
 
+[![Crates.io][crates-badge]][crates-url]
+[![Documentation][docs-badge]][docs-url]
+
+[crates-badge]: https://img.shields.io/crates/v/cosmocalc.svg
+[crates-url]: https://crates.io/crates/cosmocalc
+[docs-badge]: https://docs.rs/cosmocalc/badge.svg
+[docs-url]: https://docs.rs/cosmocalc
+
 A library for computing quantities in cosmology in the Rust programming language
 
 # Features
 
 ## Cosmological distance calculations
 
-TODO
+```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.); 
+```
 
 # Developers
 

benchmarks/README.md

@@ -0,0 +1,12 @@
+# Benchmarks
+
+We generate test vectors (and eventually benchmarks) against other projects,
+currently just astro.py.
+
+Install Python requirements:
+
+```
+python3 -m venv .venv
+source .venv/bin/activate
+pip install -r 
+```

benchmarks/astropy_vectors.py

@@ -0,0 +1,42 @@
+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()

benchmarks/requirements.txt

@@ -0,0 +1,2 @@
+astropy
+scipy

src/constants.rs

@@ -11,13 +11,14 @@ use crate::{
 
 pub const PI: f64 = std::f64::consts::PI;
 
-pub static ZERO: Lazy<DimensionlessPositiveFloat> =
-    Lazy::new(|| DimensionlessPositiveFloat::new(0.0).unwrap());
-pub static ONE: Lazy<DimensionlessPositiveFloat> =
-    Lazy::new(|| DimensionlessPositiveFloat::new(1.0).unwrap());
-
 // Neutrinos
-pub static DEFAULT_NEUTRINO_MASSES: Lazy<[eV; 3]> = Lazy::new(|| [*ZERO, *ZERO, *ZERO]);
+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)
 

src/cosmology.rs

@@ -5,16 +5,27 @@ mod omega_factors;
 pub use omega_factors::OmegaFactors;
 
 use crate::{
-    constants::{self, C_M_PER_S, DEFAULT_NEUTRINO_MASSES, DEFAULT_N_EFF, ONE, ZERO},
+    constants::{self, C_M_PER_S, DEFAULT_NEUTRINO_MASSES, DEFAULT_N_EFF},
     eV, units,
     units::{PositiveFloat, Seconds},
-    DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Kelvin, KmPerSecPerMpc, Mpc, Redshift,
+    DimensionlessFloat, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Kelvin, KmPerSecPerMpc,
+    Mpc, Redshift,
 };
 
 /// Represents an FLRW cosmology.
 ///
 /// This represents an homogenous and isotropic cosmology based
 /// on the FLRW (Friedmann-Lemaitre-Robertson-Walker) metric.
+///
+/// # Examples
+///
+/// ```
+/// use cosmocalc::{Distances, FLRWCosmology};
+///
+/// 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.);
+/// ```
 pub struct FLRWCosmology {
     /// A descriptive name.
     pub name: Option<String>,
@@ -36,6 +47,20 @@ pub struct FLRWCosmology {
 }
 
 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,
+            H_0,
+            omega,
+            T_CMB0: Some(PositiveFloat::zero()),
+            N_eff: PositiveFloat::zero(),
+            m_nu: vec![],
+        }
+    }
+
     /// Instantiate a new FLRW cosmology.
     pub fn new(
         name: Option<String>,
@@ -130,89 +155,89 @@ impl FLRWCosmology {
     /// 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 {
+    pub fn omega_gamma0(&self) -> DimensionlessFloat {
         match self.T_CMB0 {
-            Some(T_CMB0) => PositiveFloat(
+            Some(T_CMB0) => DimensionlessFloat(
                 *constants::ALPHA * T_CMB0.powf(4.)
                     / (self.critical_density(0.).0 * C_M_PER_S.powf(2.)),
             ),
-            None => *ZERO,
+            None => DimensionlessFloat::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.))
+    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.))
     }
 
     /// Dimensionless neutrino density (density/critical density) at `z=0`
-    pub fn omega_nu0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_nu0(&self) -> DimensionlessFloat {
         match self.T_CMB0 {
-            Some(_) => PositiveFloat(
+            Some(_) => DimensionlessFloat(
                 7. / 8. * (4.0f64 / 11.).powf(4. / 3.) * self.N_eff.0 * self.omega_gamma0().0,
             ),
-            None => *ZERO,
+            None => DimensionlessFloat::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.))
+    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.))
     }
 
     /// Dimensionless dark matter density (density/critical density) at `z=0`
-    pub fn omega_dm0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_dm0(&self) -> DimensionlessFloat {
         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.))
+    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.))
     }
 
     /// Dimensionless effective curvature density (density/critical density) at `z=0`
-    pub fn omega_k0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_k0(&self) -> DimensionlessFloat {
         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.))
+    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.))
     }
 
     /// Dimensionless matter density (density/critical density) at `z=0`
-    pub fn omega_m0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_m0(&self) -> DimensionlessFloat {
         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.))
+    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.))
     }
 
     /// Dimensionless baryon density (density/critical density) at `z=0`
-    pub fn omega_b0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_b0(&self) -> DimensionlessFloat {
         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.))
+    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.))
     }
 
     /// Dimensionless dark energy density (density/critical density) at `z=0`
-    pub fn omega_de0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_de0(&self) -> DimensionlessFloat {
         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.))
+    pub fn omega_de(&self, z: Redshift) -> DimensionlessFloat {
+        DimensionlessFloat(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 {
+    pub fn omega_tot0(&self) -> DimensionlessFloat {
         self.omega_m0()
             + self.omega_gamma0()
             + self.omega_nu0()
@@ -221,7 +246,7 @@ impl FLRWCosmology {
     }
 
     /// Dimensionless total density (density/critical density) at `z>0`.
-    pub fn omega_tot(&self, z: Redshift) -> DimensionlessPositiveFloat {
+    pub fn omega_tot(&self, z: Redshift) -> DimensionlessFloat {
         self.omega_m(z)
             + self.omega_gamma(z)
             + self.omega_nu(z)
@@ -231,14 +256,15 @@ impl FLRWCosmology {
 
     /// Whether this cosmology is spatially flat
     pub fn is_flat(&self) -> bool {
-        self.omega_k0() == *ZERO && self.omega_tot0() == *ONE
+        self.omega_k0() == DimensionlessFloat::zero()
+            && 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 => *ZERO,
+            None => DimensionlessPositiveFloat::zero(),
         }
     }
 }

src/cosmology/omega_factors.rs

@@ -1,13 +1,13 @@
-use crate::{units::PositiveFloat, DimensionlessPositiveFloat};
+use crate::DimensionlessFloat;
 
 /// Represents a collection of dimensionless density parameters.
 pub struct OmegaFactors {
     /// Ratio of non-relativistic matter to critical density at `z=0`.
-    pub Omega_M0: DimensionlessPositiveFloat,
+    pub Omega_M0: DimensionlessFloat,
     /// Ratio of dark energy density to critical density at `z=0`.
-    pub Omega_DE0: DimensionlessPositiveFloat,
+    pub Omega_DE0: DimensionlessFloat,
     /// Ratio of baryon density to critical density at `z=0`.
-    pub Omega_b0: DimensionlessPositiveFloat,
+    pub Omega_b0: DimensionlessFloat,
 }
 
 impl OmegaFactors {
@@ -17,26 +17,23 @@ impl OmegaFactors {
         }
 
         Ok(OmegaFactors {
-            Omega_M0: DimensionlessPositiveFloat::new(Omega_M0)?,
-            Omega_DE0: DimensionlessPositiveFloat::new(Omega_DE0)?,
-            Omega_b0: DimensionlessPositiveFloat::new(Omega_b0)?,
+            Omega_M0: DimensionlessFloat::new(Omega_M0)?,
+            Omega_DE0: DimensionlessFloat::new(Omega_DE0)?,
+            Omega_b0: DimensionlessFloat::new(Omega_b0)?,
         })
     }
 
     /// Dark matter density at `z=0` is matter at `z=0` minus baryons at `z=0`.
-    pub fn omega_dark_matter_density_0(&self) -> DimensionlessPositiveFloat {
+    pub fn omega_dark_matter_density_0(&self) -> DimensionlessFloat {
         self.Omega_M0 - self.Omega_b0
     }
 
     /// Curvature density at `z=0` given densities of relativistic particles at `z=0`.
     pub fn curvature_density_0(
         &self,
-        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
+        omega_nu0: DimensionlessFloat,
+        omega_gamma0: DimensionlessFloat,
+    ) -> DimensionlessFloat {
+        DimensionlessFloat(1.0) - self.Omega_M0 - self.Omega_DE0 - omega_nu0 - omega_gamma0
     }
 }