Feat/relativity

.gitignore

@@ -31,3 +31,4 @@ _ReSharper*/
 .vs/
 #Nuget packages folder
 packages/
+/.claude/settings.json

Numerics/NumericTest/GeneralRelativityTests.cs

@@ -0,0 +1,72 @@
+using System;
+using CSharpNumerics.Physics.Constants;
+using CSharpNumerics.Physics.Relativity;
+
+namespace NumericsTests;
+
+[TestClass]
+public class GeneralRelativityTests
+{
+    private const double SolarMass = PhysicsConstants.SolarMass;
+    private const double RadToArcsec = 180.0 / Math.PI * 3600.0;
+
+    // Mercury's orbit
+    private const double MercurySemiMajorAxis = 5.79e10;   // m
+    private const double MercuryEccentricity = 0.2056;
+    private const double MercuryPeriodDays = 87.969;
+
+    [TestMethod]
+    public void PerihelionPrecession_Mercury_PerOrbit()
+    {
+        // ≈ 5.0×10⁻⁷ rad per orbit
+        double perOrbit = GeneralRelativity.PerihelionPrecessionPerOrbit(
+            MercurySemiMajorAxis, MercuryEccentricity, SolarMass);
+        Assert.AreEqual(5.02e-7, perOrbit, 0.1e-7);
+    }
+
+    [TestMethod]
+    public void PerihelionPrecession_Mercury_IsAbout43ArcsecPerCentury()
+    {
+        double arcsec = GeneralRelativity.PerihelionPrecessionArcsecPerCentury(
+            MercurySemiMajorAxis, MercuryEccentricity, SolarMass, MercuryPeriodDays);
+        Assert.AreEqual(43.0, arcsec, 1.0);
+    }
+
+    [TestMethod]
+    public void LightDeflection_AtSunLimb_IsAbout1Point75Arcsec()
+    {
+        double sunRadius = 6.957e8; // m
+        double radians = GeneralRelativity.LightDeflection(sunRadius, SolarMass);
+        Assert.AreEqual(1.75, radians * RadToArcsec, 0.02);
+    }
+
+    [TestMethod]
+    public void LightDeflection_ScalesInverselyWithImpactParameter()
+    {
+        double a = GeneralRelativity.LightDeflection(1.0e9, SolarMass);
+        double b = GeneralRelativity.LightDeflection(2.0e9, SolarMass);
+        Assert.AreEqual(a / 2.0, b, a * 1e-12);
+    }
+
+    [TestMethod]
+    public void ShapiroDelay_IsPositive_AndGrowsWithMass()
+    {
+        double au = PhysicsConstants.AstronomicalUnit;
+        double sunRadius = 6.957e8;
+        double delay = GeneralRelativity.ShapiroDelay(au, au, sunRadius, SolarMass);
+
+        Assert.IsTrue(delay > 0.0);
+        // Earth–graze–Earth one-way excess delay is on the order of 0.1 ms.
+        Assert.AreEqual(1.2e-4, delay, 0.3e-4);
+
+        double heavier = GeneralRelativity.ShapiroDelay(au, au, sunRadius, 2.0 * SolarMass);
+        Assert.IsTrue(heavier > delay);
+    }
+
+    [TestMethod]
+    [ExpectedException(typeof(ArgumentOutOfRangeException))]
+    public void ShapiroDelay_RadiusBelowImpactParameter_Throws()
+    {
+        GeneralRelativity.ShapiroDelay(1.0e8, 1.0e12, 5.0e8, SolarMass);
+    }
+}

Numerics/NumericTest/SchwarzschildTests.cs

@@ -0,0 +1,66 @@
+using System;
+using CSharpNumerics.Physics.Constants;
+using CSharpNumerics.Physics.Relativity;
+
+namespace NumericsTests;
+
+[TestClass]
+public class SchwarzschildTests
+{
+    private const double C = PhysicsConstants.SpeedOfLight;
+    private const double SolarMass = PhysicsConstants.SolarMass;
+
+    [TestMethod]
+    public void Radius_OfSun_IsAboutTwoPointNineKilometres()
+    {
+        // r_s(Sun) ≈ 2953 m
+        Assert.AreEqual(2953.0, Schwarzschild.Radius(SolarMass), 5.0);
+    }
+
+    [TestMethod]
+    public void PhotonSphere_And_Isco_AreMultiplesOfSchwarzschildRadius()
+    {
+        double rs = Schwarzschild.Radius(SolarMass);
+        Assert.AreEqual(1.5 * rs, Schwarzschild.PhotonSphereRadius(SolarMass), rs * 1e-12);
+        Assert.AreEqual(3.0 * rs, Schwarzschild.Isco(SolarMass), rs * 1e-12);
+    }
+
+    [TestMethod]
+    public void TimeDilationFactor_FarFromMass_ApproachesOne()
+    {
+        double factor = Schwarzschild.TimeDilationFactor(1.0e12, SolarMass);
+        Assert.IsTrue(factor < 1.0);
+        Assert.AreEqual(1.0, factor, 1e-5);
+    }
+
+    [TestMethod]
+    [ExpectedException(typeof(ArgumentOutOfRangeException))]
+    public void TimeDilationFactor_AtHorizon_Throws()
+    {
+        double rs = Schwarzschild.Radius(SolarMass);
+        Schwarzschild.TimeDilationFactor(rs, SolarMass);
+    }
+
+    [TestMethod]
+    public void Redshift_IsPositive_AndSmallFarFromMass()
+    {
+        double z = Schwarzschild.Redshift(1.0e9, SolarMass);
+        Assert.IsTrue(z > 0.0);
+        Assert.IsTrue(z < 1e-5);
+    }
+
+    [TestMethod]
+    public void ClockRateRatio_LowerClockRunsSlow()
+    {
+        // Clock deeper in the well (smaller r) ticks slower than the higher one.
+        double ratio = Schwarzschild.ClockRateRatio(1.0e7, 1.0e9, SolarMass);
+        Assert.IsTrue(ratio < 1.0);
+    }
+
+    [TestMethod]
+    public void EscapeVelocity_AtSchwarzschildRadius_EqualsSpeedOfLight()
+    {
+        double rs = Schwarzschild.Radius(SolarMass);
+        Assert.AreEqual(C, Schwarzschild.EscapeVelocity(rs, SolarMass), 1.0);
+    }
+}

Numerics/NumericTest/SpecialRelativityTests.cs

@@ -0,0 +1,85 @@
+using System;
+using CSharpNumerics.Physics.Constants;
+using CSharpNumerics.Physics.Relativity;
+
+namespace NumericsTests;
+
+[TestClass]
+public class SpecialRelativityTests
+{
+    private const double C = PhysicsConstants.SpeedOfLight;
+
+    [TestMethod]
+    public void LorentzFactor_AtRest_IsOne()
+    {
+        Assert.AreEqual(1.0, SpecialRelativity.LorentzFactor(0.0), 1e-15);
+    }
+
+    [TestMethod]
+    public void LorentzFactor_At080c_IsFiveThirds()
+    {
+        // γ = 1/√(1 − 0.64) = 1/0.6 = 5/3
+        Assert.AreEqual(5.0 / 3.0, SpecialRelativity.LorentzFactor(0.8 * C), 1e-9);
+    }
+
+    [TestMethod]
+    [ExpectedException(typeof(ArgumentOutOfRangeException))]
+    public void LorentzFactor_AtLightSpeed_Throws()
+    {
+        SpecialRelativity.LorentzFactor(C);
+    }
+
+    [TestMethod]
+    public void TimeDilation_MovingClockRunsSlow()
+    {
+        // 1 s of proper time at 0.8c → 5/3 s of coordinate time
+        Assert.AreEqual(5.0 / 3.0, SpecialRelativity.TimeDilation(1.0, 0.8 * C), 1e-9);
+    }
+
+    [TestMethod]
+    public void LengthContraction_At080c_IsThreeFifths()
+    {
+        Assert.AreEqual(0.6, SpecialRelativity.LengthContraction(1.0, 0.8 * C), 1e-9);
+    }
+
+    [TestMethod]
+    public void KineticEnergy_LowSpeed_MatchesClassicalLimit()
+    {
+        // (γ − 1)mc² → ½mv² for v ≪ c
+        double m = 2.0, v = 1000.0;
+        double classical = 0.5 * m * v * v;
+        Assert.AreEqual(classical, SpecialRelativity.KineticEnergy(m, v), 1.0);
+    }
+
+    [TestMethod]
+    public void AddVelocities_HalfPlusHalf_IsBelowLight()
+    {
+        // (0.5c + 0.5c)/(1 + 0.25) = 0.8c
+        Assert.AreEqual(0.8 * C, SpecialRelativity.AddVelocities(0.5 * C, 0.5 * C), 1.0);
+    }
+
+    [TestMethod]
+    public void DopplerFactor_Receding_IsRedshift()
+    {
+        // β = 0.6 → √(0.4/1.6) = 0.5  (observed frequency halved)
+        Assert.AreEqual(0.5, SpecialRelativity.DopplerFactor(0.6 * C), 1e-12);
+    }
+
+    [TestMethod]
+    public void EnergyFromMomentum_MatchesTotalEnergy()
+    {
+        double m = 1.0, v = 0.9 * C;
+        double p = SpecialRelativity.Momentum(m, v);
+        double expected = SpecialRelativity.TotalEnergy(m, v);
+        Assert.AreEqual(expected, SpecialRelativity.EnergyFromMomentum(m, p), expected * 1e-12);
+    }
+
+    [TestMethod]
+    public void Rapidity_AddsLinearlyUnderVelocityComposition()
+    {
+        double u = 0.3 * C, v = 0.4 * C;
+        double combined = SpecialRelativity.Rapidity(SpecialRelativity.AddVelocities(u, v));
+        double sum = SpecialRelativity.Rapidity(u) + SpecialRelativity.Rapidity(v);
+        Assert.AreEqual(sum, combined, 1e-9);
+    }
+}

Numerics/Numerics/Physics/README.md

@@ -692,6 +692,109 @@ double esiVenus = AstronomyExtensions.CalculateEsi(0.95, 0.95, 0.93, 737); // 
 
 ---
 
+## 🌌 Relativity
+
+The `Physics.Relativity` namespace covers **special relativity** (flat-spacetime kinematics and dynamics) and **weak-field general relativity** (Schwarzschild geometry and the classical tests of GR). Speeds are in m/s, masses in kg, lengths in metres, angles in radians.
+
+### Special Relativity
+
+```csharp
+using CSharpNumerics.Physics.Relativity;
+using CSharpNumerics.Physics.Constants;
+
+double c = PhysicsConstants.SpeedOfLight;
+
+// Lorentz factor γ = 1/√(1 − v²/c²)
+double gamma = SpecialRelativity.LorentzFactor(0.8 * c);      // 5/3 ≈ 1.6667
+
+// Time dilation: coordinate time for 1 s of proper time at 0.8c
+double dt = SpecialRelativity.TimeDilation(1.0, 0.8 * c);     // 1.6667 s
+
+// Length contraction: L = L₀/γ
+double L = SpecialRelativity.LengthContraction(1.0, 0.8 * c); // 0.6 m
+
+// Momentum and energy
+double p  = SpecialRelativity.Momentum(mass: 1.0, velocity: 0.9 * c);  // γmv
+double E0 = SpecialRelativity.RestEnergy(1.0);                         // mc²
+double E  = SpecialRelativity.TotalEnergy(1.0, 0.9 * c);               // γmc²
+double Ek = SpecialRelativity.KineticEnergy(1.0, 0.9 * c);            // (γ−1)mc²
+double Ep = SpecialRelativity.EnergyFromMomentum(1.0, p);             // √((pc)²+(mc²)²)
+
+// Relativistic velocity addition — never exceeds c
+double s = SpecialRelativity.AddVelocities(0.5 * c, 0.5 * c);  // 0.8c
+
+// Longitudinal Doppler (positive velocity = receding → redshift)
+double factor = SpecialRelativity.DopplerFactor(0.6 * c);              // 0.5
+double fObs   = SpecialRelativity.ObservedFrequency(1.0e9, 0.6 * c);   // 0.5 GHz
+
+// Rapidity — adds linearly under collinear boosts
+double phi = SpecialRelativity.Rapidity(0.6 * c);
+```
+
+> `KineticEnergy` is evaluated as `mc²·β²/(s(1+s))` with `s = √(1−β²)` — algebraically identical to `(γ−1)mc²` but free of the catastrophic cancellation the direct form suffers at low speeds (where it correctly reduces to ½mv²).
+
+### Schwarzschild Geometry
+
+Geometry around a static, uncharged, spherically symmetric mass.
+
+```csharp
+using CSharpNumerics.Physics.Relativity;
+using CSharpNumerics.Physics.Constants;
+
+double M = PhysicsConstants.SolarMass;
+
+// Characteristic radii
+double rs    = Schwarzschild.Radius(M);              // 2GM/c² ≈ 2953 m (Sun)
+double rPh   = Schwarzschild.PhotonSphereRadius(M);  // 1.5·rs
+double rIsco = Schwarzschild.Isco(M);                // 3·rs (= 6GM/c²)
+
+// Gravitational time dilation & redshift for a static observer at radius r
+double r = 7.0e8;
+double dtau_dt = Schwarzschild.TimeDilationFactor(r, M);  // √(1 − rs/r)
+double z       = Schwarzschild.Redshift(r, M);            // 1/√(1−rs/r) − 1
+
+// Clock-rate ratio between two radii (the lower clock runs slow)
+double ratio = Schwarzschild.ClockRateRatio(radiusLower: 1e7, radiusUpper: 1e9, mass: M);
+
+// Escape velocity √(2GM/r) — equals c exactly at the horizon
+double vEsc = Schwarzschild.EscapeVelocity(rs, M);  // = c
+```
+
+### Classical Tests of General Relativity
+
+Weak-field effects: perihelion precession, light deflection, and the Shapiro delay.
+
+```csharp
+using CSharpNumerics.Physics.Relativity;
+using CSharpNumerics.Physics.Constants;
+
+double M = PhysicsConstants.SolarMass;
+
+// Perihelion precession Δϖ = 6πGM / (c²·a(1−e²))
+double perOrbit = GeneralRelativity.PerihelionPrecessionPerOrbit(
+    semiMajorAxis: 5.79e10, eccentricity: 0.2056, centralMass: M);   // ≈ 5.0e-7 rad/orbit
+
+double arcsecPerCentury = GeneralRelativity.PerihelionPrecessionArcsecPerCentury(
+    5.79e10, 0.2056, M, orbitalPeriodDays: 87.969);                  // ≈ 43″ (Mercury)
+
+// Light deflection α = 4GM / (c²·b)
+double deflection = GeneralRelativity.LightDeflection(
+    impactParameter: 6.957e8, centralMass: M);                       // ≈ 1.75″ at the Sun's limb
+
+// Shapiro time delay for a signal grazing the mass
+double au = PhysicsConstants.AstronomicalUnit;
+double delay = GeneralRelativity.ShapiroDelay(
+    r1: au, r2: au, impactParameter: 6.957e8, centralMass: M);
+```
+
+| Class | Covers |
+|-------|--------|
+| `SpecialRelativity` | Lorentz factor, time dilation, length contraction, momentum/energy, velocity addition, Doppler, rapidity |
+| `Schwarzschild` | gravitational radius, photon sphere, ISCO, time-dilation factor & redshift, clock-rate ratio, escape velocity |
+| `GeneralRelativity` | perihelion precession, light deflection, Shapiro delay |
+
+---
+
 ## 🔄 Oscillations
 
 The `Physics.Oscillations` namespace provides one-dimensional oscillator models with both **analytic** and **numerical** solutions. All oscillators implement `IOscillator` for a consistent API.