From 999e1036734840f978e6c3e9b723b03b91b5da64 Mon Sep 17 00:00:00 2001 From: redshiftzero Date: Thu, 8 Sep 2022 18:43:50 -0400 Subject: [PATCH] 0.1.0 --- Cargo.toml | 1 + README.md | 14 ++++- benchmarks/README.md | 12 +++++ benchmarks/astropy_vectors.py | 42 +++++++++++++++ benchmarks/requirements.txt | 2 + src/constants.rs | 13 ++--- src/cosmology.rs | 88 ++++++++++++++++++++------------ src/cosmology/omega_factors.rs | 27 +++++----- src/distances.rs | 93 ++++++++++++++++++++++++++++++---- src/lib.rs | 5 +- src/units.rs | 51 +++++++++++++++++++ 11 files changed, 283 insertions(+), 65 deletions(-) create mode 100644 benchmarks/README.md create mode 100644 benchmarks/astropy_vectors.py create mode 100644 benchmarks/requirements.txt diff --git a/Cargo.toml b/Cargo.toml index d946a23..b4624ca 100644 --- a/Cargo.toml +++ b/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" diff --git a/README.md b/README.md index f61dd69..2e21b71 100644 --- a/README.md +++ b/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 diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 0000000..45bae7b --- /dev/null +++ b/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 +``` diff --git a/benchmarks/astropy_vectors.py b/benchmarks/astropy_vectors.py new file mode 100644 index 0000000..4d0a5eb --- /dev/null +++ b/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() \ No newline at end of file diff --git a/benchmarks/requirements.txt b/benchmarks/requirements.txt new file mode 100644 index 0000000..4460a4e --- /dev/null +++ b/benchmarks/requirements.txt @@ -0,0 +1,2 @@ +astropy +scipy diff --git a/src/constants.rs b/src/constants.rs index 3df583e..d3c52bb 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -11,13 +11,14 @@ use crate::{ pub const PI: f64 = std::f64::consts::PI; -pub static ZERO: Lazy = - Lazy::new(|| DimensionlessPositiveFloat::new(0.0).unwrap()); -pub static ONE: Lazy = - 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 = Lazy::new(|| DimensionlessPositiveFloat::new(3.04).unwrap()); // WMAP (ApJ Spergel et al 2007) diff --git a/src/cosmology.rs b/src/cosmology.rs index 2296583..b89b37a 100644 --- a/src/cosmology.rs +++ b/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, @@ -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, @@ -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(), } } } diff --git a/src/cosmology/omega_factors.rs b/src/cosmology/omega_factors.rs index f3e9aed..79cc44e 100644 --- a/src/cosmology/omega_factors.rs +++ b/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 } } diff --git a/src/distances.rs b/src/distances.rs index d5faa6c..c246083 100644 --- a/src/distances.rs +++ b/src/distances.rs @@ -1,4 +1,4 @@ -use crate::{units::PositiveFloat, FLRWCosmology, Mpc, Redshift}; +use crate::{units::PositiveFloat, DimensionlessFloat, FLRWCosmology, Mpc, Redshift}; /// Bin width in redshift integrals. const DZ: f64 = 0.0001; @@ -27,8 +27,27 @@ impl Distances for FLRWCosmology { PositiveFloat(self.hubble_distance() * integrand) } - fn transverse_comoving_distance(&self, _z: Redshift) -> Mpc { - todo!() + fn transverse_comoving_distance(&self, z: Redshift) -> Mpc { + let radial_comoving = self.radial_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(); + PositiveFloat( + self.hubble_distance() * 1. / sqrt_omega_k + * f64::sinh(sqrt_omega_k * radial_comoving.0 / self.hubble_distance()), + ) + } else if omega_k == DimensionlessFloat::zero() { + // Flat + radial_comoving + } 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()), + ) + } } fn angular_diameter_distance(&self, z: Redshift) -> Mpc { @@ -43,24 +62,29 @@ impl Distances for FLRWCosmology { #[cfg(test)] mod tests { - use crate::{constants::ZERO, cosmology::OmegaFactors}; + use crate::{cosmology::OmegaFactors, DimensionlessPositiveFloat}; use super::*; #[test] fn flat_universe_distances_no_relativistic_contribution() { + // TESTED vs: astro.py 5.1 FlatLambdaCDM 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 - // Ned Wright calculator yields 6481.1 assert!(cosmology.radial_comoving_distance(3.0).0 > 6482.5); - assert!(cosmology.radial_comoving_distance(3.0).0 < 6483.0); + 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); } #[test] fn flat_universe_distances_with_radiation_but_no_neutrinos() { - let omegas = OmegaFactors::new(0.25, 0.7, 0.05).unwrap(); + // TESTED vs: astro.py 5.1 FlatLambdaCDM + let omegas = OmegaFactors::new(0.299, 0.7, 0.05).unwrap(); let cosmology = FLRWCosmology::new( None, None, @@ -73,8 +97,12 @@ mod tests { .unwrap(); // Megaparsecs - assert!(cosmology.radial_comoving_distance(3.0).0 > 6598.5); - assert!(cosmology.radial_comoving_distance(3.0).0 < 6599.0); + 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.); } #[test] @@ -87,12 +115,57 @@ mod tests { omegas, Some(2.7255), Some(PositiveFloat(3.04)), - Some(vec![*ZERO, *ZERO, *ZERO]), + Some(vec![ + DimensionlessPositiveFloat::zero(), + DimensionlessPositiveFloat::zero(), + DimensionlessPositiveFloat::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.); + } + + #[test] + fn open_universe_distances_no_relativistic_contribution() { + let omegas = OmegaFactors::new(0.286, 0.0, 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 > 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.); + // No k-corrections here + assert!(cosmology.luminosity_distance(3.0).0 > 22000.); + assert!(cosmology.luminosity_distance(3.0).0 < 24000.); + } + + #[test] + fn closed_universe_distances_no_relativistic_contribution() { + let omegas = OmegaFactors::new(0.286, 0.8, 0.05).unwrap(); + 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.); + // No k-corrections here + assert!(cosmology.luminosity_distance(2.0).0 > 14000.); + assert!(cosmology.luminosity_distance(2.0).0 < 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.); } } diff --git a/src/lib.rs b/src/lib.rs index ff83d09..c22cf3e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -8,10 +8,11 @@ pub mod units; pub mod utils; pub use cosmology::FLRWCosmology; +pub use distances::Distances; // Common units are re-exported from the crate root for convenience. pub use redshift::Redshift; pub use units::{ - eV, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Joule, Kelvin, Kilogram, KmPerSecPerMpc, - Mpc, + eV, DimensionlessFloat, DimensionlessPositiveFloat, HInvKmPerSecPerMpc, Joule, Kelvin, + Kilogram, KmPerSecPerMpc, Mpc, }; diff --git a/src/units.rs b/src/units.rs index 487d0ec..dab463d 100644 --- a/src/units.rs +++ b/src/units.rs @@ -46,6 +46,49 @@ 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 { + 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) + } +} + /// Represents continuous physical quantities that _cannot_ be negative. #[derive(Clone, Copy, Debug, PartialOrd, PartialEq)] pub struct PositiveFloat(pub f64); @@ -58,6 +101,14 @@ impl PositiveFloat { 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()