diff --git a/Numerics/NumericTest/FlrwModelTests.cs b/Numerics/NumericTest/FlrwModelTests.cs new file mode 100644 index 0000000..3ffbc8a --- /dev/null +++ b/Numerics/NumericTest/FlrwModelTests.cs @@ -0,0 +1,116 @@ +using System; +using CSharpNumerics.Physics.Constants; +using CSharpNumerics.Physics.Cosmology; + +namespace NumericsTests; + +[TestClass] +public class FlrwModelTests +{ + private const double Mpc = PhysicsConstants.Megaparsec; + private const double SecondsPerYear = 3.15576e7; // Julian year + + private static FlrwModel Planck() => FlrwModel.Planck2018(); + + [TestMethod] + public void Planck2018_HasFlatGeometry_AndExpectedParameters() + { + var m = Planck(); + Assert.AreEqual(67.4, m.HubbleConstantKmSMpc, 1e-9); + Assert.AreEqual(0.315, m.OmegaMatter, 1e-9); + Assert.AreEqual(0.685, m.OmegaLambda, 1e-9); + Assert.AreEqual(0.0, m.OmegaCurvature, 1e-9); // flat + } + + [TestMethod] + public void DimensionlessHubble_AtPresent_IsOne() + { + Assert.AreEqual(1.0, Planck().DimensionlessHubble(0.0), 1e-12); + } + + [TestMethod] + public void HubbleParameter_AtPresent_EqualsHubbleConstant() + { + var m = Planck(); + Assert.AreEqual(m.HubbleConstant, m.HubbleParameter(0.0), m.HubbleConstant * 1e-12); + Assert.IsTrue(m.HubbleParameter(1.0) > m.HubbleParameter(0.0)); // expansion was faster + } + + [TestMethod] + public void CriticalDensity_Today_IsAbout8Point5e_minus27() + { + // ρ_c,0 ≈ 8.5×10⁻²⁷ kg/m³ for H₀ = 67.4 + Assert.AreEqual(8.5e-27, Planck().CriticalDensity(), 0.3e-27); + } + + [TestMethod] + public void AgeOfUniverse_IsAbout13Point8Gyr() + { + double ageYears = Planck().AgeOfUniverse() / SecondsPerYear; + Assert.AreEqual(13.8e9, ageYears, 0.3e9); + } + + [TestMethod] + public void ComovingDistance_AtZ1_IsAbout3400Mpc() + { + double dcMpc = Planck().ComovingDistance(1.0) / Mpc; + Assert.AreEqual(3402.0, dcMpc, 30.0); + } + + [TestMethod] + public void ComovingDistance_AtLowRedshift_FollowsHubbleLaw() + { + // D_C(z) → D_H·z as z → 0 + var m = Planck(); + double dc = m.ComovingDistance(0.001); + Assert.AreEqual(m.HubbleDistance * 0.001, dc, m.HubbleDistance * 0.001 * 0.005); + } + + [TestMethod] + public void DistanceDuality_LuminosityEqualsOnePlusZSquaredTimesAngular() + { + // Etherington relation: D_L = (1+z)²·D_A + var m = Planck(); + double z = 2.0; + double dl = m.LuminosityDistance(z); + double da = m.AngularDiameterDistance(z); + Assert.AreEqual((1.0 + z) * (1.0 + z) * da, dl, dl * 1e-9); + } + + [TestMethod] + public void DistanceModulus_AtZ1_IsAbout44() + { + Assert.AreEqual(44.16, Planck().DistanceModulus(1.0), 0.1); + } + + [TestMethod] + public void ScaleFactor_IsInverseOfOnePlusZ() + { + Assert.AreEqual(0.5, Planck().ScaleFactor(1.0), 1e-12); + } + + [TestMethod] + public void LookbackTimePlusAge_EqualsAgeOfUniverse() + { + var m = Planck(); + double total = m.AgeOfUniverse(); + double sum = m.LookbackTime(1.0) + m.Age(1.0); + Assert.AreEqual(total, sum, total * 1e-3); + } + + [TestMethod] + public void OpenUniverse_TransverseDistanceExceedsComoving() + { + // Ω_k > 0 (open) → sinh form makes D_M > D_C + var open = new FlrwModel(70.0, 0.3, 0.0); // Ω_k = 0.7 + Assert.IsTrue(open.OmegaCurvature > 0.0); + Assert.IsTrue(open.TransverseComovingDistance(2.0) > open.ComovingDistance(2.0)); + } + + [TestMethod] + [ExpectedException(typeof(ArgumentOutOfRangeException))] + public void ComovingDistance_NegativeRedshift_Throws() + { + Planck().ComovingDistance(-0.5); + } +} diff --git a/Numerics/Numerics/Physics/Constants/PhysicsConstants.cs b/Numerics/Numerics/Physics/Constants/PhysicsConstants.cs index 2a69f7c..9080d4a 100644 --- a/Numerics/Numerics/Physics/Constants/PhysicsConstants.cs +++ b/Numerics/Numerics/Physics/Constants/PhysicsConstants.cs @@ -55,8 +55,15 @@ public static class PhysicsConstants public const double AstronomicalUnit = 1.495978707e11; // meters public const double LightYear = 9.4607e15; // meters public const double Parsec = 3.0857e16; // meters + public const double Megaparsec = 3.0857e22; // meters (10^6 parsecs) public const double EarthRadius = 6.371e6; // meters public const double MoonRadius = 1.7371e6; // meters public const double SolarLuminosity = 3.828e26; // Watts public const double HubbleConstant = 2.2e-18; // 1/s + public const double CmbTemperature = 2.72548; // Kelvin (cosmic microwave background) + + // Planck units + public const double PlanckMass = 2.176434e-8; // kg + public const double PlanckLength = 1.616255e-35; // meters + public const double PlanckTime = 5.391247e-44; // seconds } diff --git a/Numerics/Numerics/Physics/Cosmology/FlrwModel.cs b/Numerics/Numerics/Physics/Cosmology/FlrwModel.cs new file mode 100644 index 0000000..2b40242 --- /dev/null +++ b/Numerics/Numerics/Physics/Cosmology/FlrwModel.cs @@ -0,0 +1,211 @@ +using System; +using CSharpNumerics.Physics.Constants; + +namespace CSharpNumerics.Physics.Cosmology; + +/// +/// A homogeneous, isotropic Friedmann–Lemaître–Robertson–Walker (FLRW) cosmology. +/// Holds the present-day density parameters and Hubble constant and derives the +/// expansion history H(z), the cosmological distance measures, lookback time and +/// the age of the universe by numerical integration of the Friedmann equation. +/// +/// Conventions: H₀ is supplied in km·s⁻¹·Mpc⁻¹; distances are returned in metres, +/// times in seconds. The curvature density is fixed by flatness: +/// Ω_k = 1 − Ω_m − Ω_Λ − Ω_r. +/// +public class FlrwModel +{ + private const double C = PhysicsConstants.SpeedOfLight; + private const double G = PhysicsConstants.GravitationalConstant; + private const double Mpc = PhysicsConstants.Megaparsec; + + private const int IntegrationSteps = 2000; + private const double ScaleFactorFloor = 1e-8; // lower bound for age integrals + + /// Hubble constant H₀ in s⁻¹. + public double HubbleConstant { get; } + + /// Hubble constant H₀ in km·s⁻¹·Mpc⁻¹. + public double HubbleConstantKmSMpc { get; } + + /// Present-day matter density parameter Ω_m. + public double OmegaMatter { get; } + + /// Present-day dark-energy density parameter Ω_Λ. + public double OmegaLambda { get; } + + /// Present-day radiation density parameter Ω_r. + public double OmegaRadiation { get; } + + /// Curvature density parameter Ω_k = 1 − Ω_m − Ω_Λ − Ω_r. + public double OmegaCurvature { get; } + + /// Hubble distance D_H = c/H₀ (metres). + public double HubbleDistance => C / HubbleConstant; + + /// Hubble time 1/H₀ (seconds). + public double HubbleTime => 1.0 / HubbleConstant; + + /// + /// Creates an FLRW cosmology. + /// + /// H₀ in km·s⁻¹·Mpc⁻¹. + /// Present-day matter density parameter Ω_m. + /// Present-day dark-energy density parameter Ω_Λ. + /// Present-day radiation density parameter Ω_r (default 0). + public FlrwModel(double hubbleConstantKmSMpc, double omegaMatter, double omegaLambda, double omegaRadiation = 0.0) + { + if (hubbleConstantKmSMpc <= 0) + throw new ArgumentOutOfRangeException(nameof(hubbleConstantKmSMpc), "H₀ must be positive."); + + HubbleConstantKmSMpc = hubbleConstantKmSMpc; + HubbleConstant = hubbleConstantKmSMpc * 1000.0 / Mpc; // (km/s)/Mpc → 1/s + OmegaMatter = omegaMatter; + OmegaLambda = omegaLambda; + OmegaRadiation = omegaRadiation; + OmegaCurvature = 1.0 - omegaMatter - omegaLambda - omegaRadiation; + } + + /// + /// The Planck 2018 flat ΛCDM concordance cosmology + /// (H₀ = 67.4, Ω_m = 0.315, Ω_Λ = 0.685). + /// + public static FlrwModel Planck2018() => new FlrwModel(67.4, 0.315, 0.685); + + /// + /// Dimensionless Hubble function E(z) = H(z)/H₀ + /// = √(Ω_m(1+z)³ + Ω_r(1+z)⁴ + Ω_k(1+z)² + Ω_Λ). + /// + public double DimensionlessHubble(double z) + { + double zp1 = 1.0 + z; + double zp1_2 = zp1 * zp1; + double zp1_3 = zp1_2 * zp1; + double zp1_4 = zp1_3 * zp1; + + return Math.Sqrt( + OmegaMatter * zp1_3 + + OmegaRadiation * zp1_4 + + OmegaCurvature * zp1_2 + + OmegaLambda); + } + + /// Hubble parameter H(z) = H₀·E(z) in s⁻¹. + public double HubbleParameter(double z) => HubbleConstant * DimensionlessHubble(z); + + /// Scale factor a = 1/(1+z) (a = 1 today). + public double ScaleFactor(double z) => 1.0 / (1.0 + z); + + /// + /// Critical density ρ_c(z) = 3H(z)²/(8πG) in kg·m⁻³. + /// With no argument, the present-day value ρ_c,0. + /// + public double CriticalDensity(double z = 0.0) + { + double h = HubbleParameter(z); + return 3.0 * h * h / (8.0 * Math.PI * G); + } + + /// + /// Line-of-sight comoving distance D_C = D_H·∫₀ᶻ dz'/E(z') in metres. + /// + public double ComovingDistance(double z) + { + if (z < 0) throw new ArgumentOutOfRangeException(nameof(z), "Redshift must be non-negative."); + if (z == 0) return 0.0; + + double integral = Simpson(zp => 1.0 / DimensionlessHubble(zp), 0.0, z, IntegrationSteps); + return HubbleDistance * integral; + } + + /// + /// Transverse comoving distance D_M (metres), accounting for spatial curvature. + /// Equals D_C for a flat universe; uses the sinh/sin forms otherwise. + /// + public double TransverseComovingDistance(double z) + { + double dc = ComovingDistance(z); + double dh = HubbleDistance; + + if (Math.Abs(OmegaCurvature) < 1e-12) + return dc; + + double sqrtOk = Math.Sqrt(Math.Abs(OmegaCurvature)); + if (OmegaCurvature > 0) // open + return dh / sqrtOk * Math.Sinh(sqrtOk * dc / dh); + + return dh / sqrtOk * Math.Sin(sqrtOk * dc / dh); // closed + } + + /// + /// Luminosity distance D_L = (1+z)·D_M in metres. + /// + public double LuminosityDistance(double z) => (1.0 + z) * TransverseComovingDistance(z); + + /// + /// Angular-diameter distance D_A = D_M/(1+z) in metres. + /// + public double AngularDiameterDistance(double z) => TransverseComovingDistance(z) / (1.0 + z); + + /// + /// Distance modulus μ = 5·log₁₀(D_L / 10 pc) in magnitudes. + /// + public double DistanceModulus(double z) + { + double dLparsecs = LuminosityDistance(z) / PhysicsConstants.Parsec; + return 5.0 * Math.Log10(dLparsecs) - 5.0; + } + + /// + /// Age of the universe at redshift in seconds: + /// t(a) = (1/H₀)·∫₀ᵃ da'/(a'·E(a')), with a = 1/(1+z). + /// + public double Age(double z) + { + if (z < 0) throw new ArgumentOutOfRangeException(nameof(z), "Redshift must be non-negative."); + return AgeIntegral(ScaleFactorFloor, ScaleFactor(z)); + } + + /// Present-day age of the universe t₀ in seconds. + public double AgeOfUniverse() => AgeIntegral(ScaleFactorFloor, 1.0); + + /// + /// Lookback time to redshift in seconds: + /// the difference between the present age and the age at that redshift. + /// + public double LookbackTime(double z) + { + if (z < 0) throw new ArgumentOutOfRangeException(nameof(z), "Redshift must be non-negative."); + return AgeIntegral(ScaleFactor(z), 1.0); + } + + // t = (1/H₀)·∫ da/(a·E(a)) with E(a) = E(z) at a = 1/(1+z). + private double AgeIntegral(double aLow, double aHigh) + { + if (aHigh <= aLow) return 0.0; + + double integral = Simpson(a => + { + double zp = 1.0 / a - 1.0; + return 1.0 / (a * DimensionlessHubble(zp)); + }, aLow, aHigh, IntegrationSteps); + + return integral / HubbleConstant; + } + + // Composite Simpson's rule over [a, b] with an even number of intervals. + private static double Simpson(Func f, double a, double b, int n) + { + if (n % 2 != 0) n++; + double h = (b - a) / n; + double sum = f(a) + f(b); + + for (int i = 1; i < n; i++) + { + double x = a + i * h; + sum += (i % 2 == 0 ? 2.0 : 4.0) * f(x); + } + + return sum * h / 3.0; + } +} diff --git a/Numerics/Numerics/Physics/README.md b/Numerics/Numerics/Physics/README.md index 5df54fe..03ed216 100644 --- a/Numerics/Numerics/Physics/README.md +++ b/Numerics/Numerics/Physics/README.md @@ -692,6 +692,73 @@ double esiVenus = AstronomyExtensions.CalculateEsi(0.95, 0.95, 0.93, 737); // --- +## 🌌 Cosmology + +The `Physics.Cosmology` namespace provides the `FlrwModel` — a homogeneous, isotropic Friedmann–Lemaître–Robertson–Walker cosmology. It derives the expansion history H(z), the cosmological distance measures, lookback time and the age of the universe by numerically integrating the Friedmann equation. H₀ is given in km·s⁻¹·Mpc⁻¹; distances are returned in metres, times in seconds. + +### Building a model + +```csharp +using CSharpNumerics.Physics.Cosmology; + +// Planck 2018 flat ΛCDM (H₀ = 67.4, Ω_m = 0.315, Ω_Λ = 0.685) +var cosmos = FlrwModel.Planck2018(); + +// Or specify your own — curvature is fixed by flatness: +// Ω_k = 1 − Ω_m − Ω_Λ − Ω_r +var open = new FlrwModel(hubbleConstantKmSMpc: 70.0, omegaMatter: 0.3, omegaLambda: 0.0); // Ω_k = 0.7 + +double dh = cosmos.HubbleDistance; // c/H₀ (m) +double th = cosmos.HubbleTime; // 1/H₀ (s) +double ok = cosmos.OmegaCurvature; // ≈ 0 for the flat default +``` + +### Expansion history & density + +```csharp +// Dimensionless Hubble E(z) = H(z)/H₀ = √(Ω_m(1+z)³ + Ω_r(1+z)⁴ + Ω_k(1+z)² + Ω_Λ) +double E = cosmos.DimensionlessHubble(1.0); +double Hz = cosmos.HubbleParameter(1.0); // H(z) in s⁻¹ + +double a = cosmos.ScaleFactor(1.0); // 1/(1+z) = 0.5 + +// Critical density ρ_c = 3H²/(8πG) — today (≈ 8.5×10⁻²⁷ kg/m³) or at redshift z +double rhoC = cosmos.CriticalDensity(); +double rhoCz = cosmos.CriticalDensity(2.0); +``` + +### Cosmological distances + +```csharp +double dC = cosmos.ComovingDistance(1.0); // line-of-sight comoving (m) +double dM = cosmos.TransverseComovingDistance(1.0); // accounts for curvature +double dL = cosmos.LuminosityDistance(1.0); // (1+z)·D_M +double dA = cosmos.AngularDiameterDistance(1.0); // D_M/(1+z) +double mu = cosmos.DistanceModulus(1.0); // 5·log₁₀(D_L/10pc) ≈ 44.2 at z = 1 + +// At low redshift D_C → (c/H₀)·z (Hubble's law); the Etherington relation +// D_L = (1+z)²·D_A holds exactly. +``` + +### Cosmic time + +```csharp +double t0 = cosmos.AgeOfUniverse(); // ≈ 13.8 Gyr (in seconds) +double tz = cosmos.Age(2.0); // age of the universe at z = 2 +double look = cosmos.LookbackTime(1.0); // t₀ − t(z), light-travel time since z = 1 +``` + +| Member | Returns | +|--------|---------| +| `DimensionlessHubble(z)` / `HubbleParameter(z)` | E(z) / H(z) | +| `CriticalDensity(z = 0)` | ρ_c = 3H²/(8πG) | +| `ComovingDistance` / `TransverseComovingDistance` | D_C / D_M | +| `LuminosityDistance` / `AngularDiameterDistance` | D_L / D_A | +| `DistanceModulus` | μ = 5·log₁₀(D_L/10pc) | +| `AgeOfUniverse` / `Age(z)` / `LookbackTime(z)` | cosmic time measures | + +--- + ## 🔄 Oscillations The `Physics.Oscillations` namespace provides one-dimensional oscillator models with both **analytic** and **numerical** solutions. All oscillators implement `IOscillator` for a consistent API.