lookback time and distance, comoving volume (#1)

* refactoring units, add my cool macro lol

* more refactoring

* structs for all base units

* use CODATA 2018 values for constants

* lookback time

* lookback distance

* comoving volume

* add Rust benchmarks

* initial profiling improvements

* lots of perf improvements

* fixins

* vectorizing distances

* final bits of optimization

* changelog

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.

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

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`.

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);

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()

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 {