use core::f64::consts::PI;
use metrology::{C, G, HBAR};
pub fn graviton_count(h0: f64, f: f64) -> f64 {
let w = 2.0 * PI * f;
h0 * h0 * C.powi(5) / (32.0 * PI * G * HBAR * w * w)
}
pub fn occupation_per_cell(h: f64, f: f64) -> f64 {
let w = 2.0 * PI * f;
C * C * h * h * w * w / (32.0 * PI * G) / (HBAR * w) * (C / f).powi(3)
}
pub fn stimulated_rate(vs: f64, m: f64, h: f64, l: f64) -> f64 {
vs * vs * m * h * h / (4.0 * l * l * PI.powi(3) * HBAR)
}
pub fn spontaneous_rate(rho: f64, vs: f64, r: f64, l: f64) -> f64 {
8.0 * PI * G * rho * vs.powi(4) * r * r / (l * C.powi(5))
}
pub fn optimal_bar_mass(vs: f64, h0: f64, chirp: f64, f: f64) -> f64 {
let w = 2.0 * PI * f;
24.0 * PI * PI / 5.0 * HBAR / (h0 * h0 * vs * vs) * (G * chirp / (2.0 * C.powi(3))).powf(5.0 / 3.0) * w.powf(8.0 / 3.0)
}
pub fn bar_cross_section(m: f64, vs: f64, l: f64, r: f64, poisson: f64, n: f64, theta: f64) -> f64 {
let correction = 1.0 + poisson * (1.0 - 2.0 * poisson) * (n * PI * r / l).powi(2) / 2.0;
8.0 / PI * m / (n * n) * (G / C) * (vs / C).powi(2) * theta.sin().powi(4) * correction
}
#[cfg(test)]
mod tests {
use super::*;
use crate::inspiral::energy_spectrum;
use core::f64::consts::FRAC_PI_2;
use metrology::{MPC, SOLAR_MASS};
fn close(a: f64, b: f64, rel: f64) -> bool {
(a / b - 1.0).abs() < rel
}
#[test]
fn counting_gravitons() {
assert!(close(graviton_count(1e-21, 150.0), 4e36, 0.05));
assert!(close(occupation_per_cell(1e-22, 1e3), 2.15031e35, 1e-5));
assert!(close(occupation_per_cell(7e-28, 1e18), 1.05365e-5, 1e-5));
assert!(close(occupation_per_cell(1e-30, 4.63715e12), 1.0, 1e-5));
}
#[test]
fn tobar_rates_and_masses() {
assert!(close(stimulated_rate(5400.0, 1800.0, 5e-22, 1.0), 1.00326, 1e-5)); assert!(close(spontaneous_rate(8570.0, 5e3, 0.5, 1.0), 9.27562e-34, 1e-5)); let be = optimal_bar_mass(13e3, 2e-22, 1.19 * SOLAR_MASS, 100.0);
let al = optimal_bar_mass(5.4e3, 2e-22, 1.19 * SOLAR_MASS, 150.0);
assert!(close(be, 12.8458, 1e-5) && close(al, 219.501, 1e-5));
assert!(close(be, 15.0, 0.25) && close(al, 250.0, 0.25)); }
#[test]
fn gw170817_in_a_bar() {
let sigma = bar_cross_section(2300.0, 5400.0, 3.0, 0.6, 0.345, 1.0, FRAC_PI_2);
assert!(close(sigma, 4.31988e-25, 1e-5)); let fluence = energy_spectrum(1.188 * SOLAR_MASS, 900.0) / (4.0 * PI * (40.0 * MPC).powi(2));
let e = sigma * fluence;
assert!(close(e, 1.15207e-31, 1e-5));
assert!(close(e / (HBAR * 2.0 * PI * 900.0), 0.193189, 1e-5));
assert!(close(e / (C * 2300.0), 1.67083e-43, 1e-5));
}
}