diff --git a/.gitignore b/.gitignore index b02de1b..f93f341 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,4 @@ _ReSharper*/ .vs/ #Nuget packages folder packages/ +/.claude/settings.json diff --git a/Numerics/NumericTest/GeneralRelativityTests.cs b/Numerics/NumericTest/GeneralRelativityTests.cs new file mode 100644 index 0000000..433d62e --- /dev/null +++ b/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); + } +} diff --git a/Numerics/NumericTest/SchwarzschildTests.cs b/Numerics/NumericTest/SchwarzschildTests.cs new file mode 100644 index 0000000..58c5dcc --- /dev/null +++ b/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); + } +} diff --git a/Numerics/NumericTest/SpecialRelativityTests.cs b/Numerics/NumericTest/SpecialRelativityTests.cs new file mode 100644 index 0000000..4a686fb --- /dev/null +++ b/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); + } +} diff --git a/Numerics/Numerics/Physics/README.md b/Numerics/Numerics/Physics/README.md index 5df54fe..bd2649a 100644 --- a/Numerics/Numerics/Physics/README.md +++ b/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. diff --git a/Numerics/Numerics/Physics/Relativity/GeneralRelativity.cs b/Numerics/Numerics/Physics/Relativity/GeneralRelativity.cs new file mode 100644 index 0000000..5b7c256 --- /dev/null +++ b/Numerics/Numerics/Physics/Relativity/GeneralRelativity.cs @@ -0,0 +1,99 @@ +using System; +using CSharpNumerics.Physics.Constants; + +namespace CSharpNumerics.Physics.Relativity; + +/// +/// Weak-field general-relativistic effects — the classical tests of GR: +/// perihelion precession, gravitational light deflection, and the Shapiro +/// time delay. Masses are in kg, lengths in metres, angles in radians. +/// +public static class GeneralRelativity +{ + private const double G = PhysicsConstants.GravitationalConstant; + private const double C = PhysicsConstants.SpeedOfLight; + private const double C2 = C * C; + private const double C3 = C * C * C; + + private const double RadiansToArcseconds = 180.0 / Math.PI * 3600.0; + + /// + /// Relativistic perihelion precession per orbit (radians): + /// Δϖ = 6πGM / (c²·a·(1 − e²)). + /// For Mercury this is ≈ 5.0×10⁻⁷ rad per orbit (≈ 43″ per century). + /// + /// Orbital semi-major axis a (m). + /// Orbital eccentricity e, 0 ≤ e < 1. + /// Mass of the central body (kg). + public static double PerihelionPrecessionPerOrbit(double semiMajorAxis, double eccentricity, double centralMass) + { + if (semiMajorAxis <= 0) + throw new ArgumentOutOfRangeException(nameof(semiMajorAxis), "Semi-major axis must be positive."); + if (eccentricity < 0 || eccentricity >= 1) + throw new ArgumentOutOfRangeException(nameof(eccentricity), "Eccentricity must be in [0, 1)."); + + double latusFactor = semiMajorAxis * (1.0 - eccentricity * eccentricity); + return 6.0 * Math.PI * G * centralMass / (C2 * latusFactor); + } + + /// + /// Relativistic perihelion precession in arcseconds per century, given the + /// orbital period. Convenience wrapper around + /// . + /// + /// Semi-major axis a (m). + /// Eccentricity e, 0 ≤ e < 1. + /// Central mass (kg). + /// Orbital period (days). + public static double PerihelionPrecessionArcsecPerCentury( + double semiMajorAxis, double eccentricity, double centralMass, double orbitalPeriodDays) + { + if (orbitalPeriodDays <= 0) + throw new ArgumentOutOfRangeException(nameof(orbitalPeriodDays), "Period must be positive."); + + double perOrbit = PerihelionPrecessionPerOrbit(semiMajorAxis, eccentricity, centralMass); + double orbitsPerCentury = 36525.0 / orbitalPeriodDays; + return perOrbit * orbitsPerCentury * RadiansToArcseconds; + } + + /// + /// Gravitational deflection of light passing a mass at impact parameter + /// : α = 4GM / (c²·b) radians. + /// For a ray grazing the Sun's limb this is ≈ 1.75″. + /// + /// Closest approach distance b (m). + /// Deflecting mass (kg). + public static double LightDeflection(double impactParameter, double centralMass) + { + if (impactParameter <= 0) + throw new ArgumentOutOfRangeException(nameof(impactParameter), "Impact parameter must be positive."); + + return 4.0 * G * centralMass / (C2 * impactParameter); + } + + /// + /// Shapiro time delay (seconds) — the extra one-way light-travel time for a + /// signal passing a mass at impact parameter + /// , between a source at radial distance + /// and a receiver at : + /// Δt = (2GM/c³)·ln[ (r₁+√(r₁²−b²))(r₂+√(r₂²−b²)) / b² ]. + /// + /// Radial distance of the source from the mass (m). + /// Radial distance of the receiver from the mass (m). + /// Impact parameter b of the signal path (m). + /// Mass causing the delay (kg). + public static double ShapiroDelay(double r1, double r2, double impactParameter, double centralMass) + { + if (impactParameter <= 0) + throw new ArgumentOutOfRangeException(nameof(impactParameter), "Impact parameter must be positive."); + if (r1 < impactParameter || r2 < impactParameter) + throw new ArgumentOutOfRangeException(nameof(r1), + "Both radial distances must be at least the impact parameter."); + + double x1 = Math.Sqrt(r1 * r1 - impactParameter * impactParameter); + double x2 = Math.Sqrt(r2 * r2 - impactParameter * impactParameter); + double argument = (r1 + x1) * (r2 + x2) / (impactParameter * impactParameter); + + return 2.0 * G * centralMass / C3 * Math.Log(argument); + } +} diff --git a/Numerics/Numerics/Physics/Relativity/Schwarzschild.cs b/Numerics/Numerics/Physics/Relativity/Schwarzschild.cs new file mode 100644 index 0000000..1ab079e --- /dev/null +++ b/Numerics/Numerics/Physics/Relativity/Schwarzschild.cs @@ -0,0 +1,92 @@ +using System; +using CSharpNumerics.Physics.Constants; + +namespace CSharpNumerics.Physics.Relativity; + +/// +/// Geometry of the Schwarzschild solution — the spacetime around a static, +/// uncharged, spherically symmetric mass. Provides the characteristic radii +/// and the gravitational time-dilation / redshift experienced by a static +/// observer outside the mass. Masses are in kg, radii in metres. +/// +public static class Schwarzschild +{ + private const double G = PhysicsConstants.GravitationalConstant; + private const double C = PhysicsConstants.SpeedOfLight; + private const double C2 = C * C; + + /// + /// Schwarzschild (gravitational) radius r_s = 2GM/c². The event horizon of + /// a non-rotating black hole of this mass. + /// + public static double Radius(double mass) + { + if (mass < 0) throw new ArgumentOutOfRangeException(nameof(mass), "Mass must be non-negative."); + return 2.0 * G * mass / C2; + } + + /// + /// Photon-sphere radius r = 3GM/c² = 1.5·r_s, where light can orbit on + /// (unstable) circular null geodesics. + /// + public static double PhotonSphereRadius(double mass) + { + return 1.5 * Radius(mass); + } + + /// + /// Innermost stable circular orbit (ISCO) for a massive test particle, + /// r = 6GM/c² = 3·r_s. Inside this radius no stable circular orbit exists. + /// + public static double Isco(double mass) + { + return 3.0 * Radius(mass); + } + + /// + /// Gravitational time-dilation factor dτ/dt = √(1 − r_s/r) for a static + /// observer at : their proper time runs slow + /// relative to a clock at infinity. + /// + /// Areal radius (m), must be greater than r_s. + /// Central mass (kg). + /// If r ≤ r_s. + public static double TimeDilationFactor(double radius, double mass) + { + double rs = Radius(mass); + if (radius <= rs) + throw new ArgumentOutOfRangeException(nameof(radius), + "Radius must be outside the Schwarzschild radius."); + + return Math.Sqrt(1.0 - rs / radius); + } + + /// + /// Gravitational redshift z = 1/√(1 − r_s/r) − 1 of light emitted by a + /// static source at and received at infinity. + /// + public static double Redshift(double radius, double mass) + { + return 1.0 / TimeDilationFactor(radius, mass) - 1.0; + } + + /// + /// Ratio of clock rates between two static radii, dτ(r1)/dτ(r2) + /// = √((1 − r_s/r1)/(1 − r_s/r2)). Values < 1 mean the clock at + /// runs slow relative to the upper one. + /// + public static double ClockRateRatio(double radiusLower, double radiusUpper, double mass) + { + return TimeDilationFactor(radiusLower, mass) / TimeDilationFactor(radiusUpper, mass); + } + + /// + /// Newtonian escape velocity √(2GM/r). In the Schwarzschild geometry this + /// equals c exactly at the event horizon (r = r_s). + /// + public static double EscapeVelocity(double radius, double mass) + { + if (radius <= 0) throw new ArgumentOutOfRangeException(nameof(radius), "Radius must be positive."); + return Math.Sqrt(2.0 * G * mass / radius); + } +} diff --git a/Numerics/Numerics/Physics/Relativity/SpecialRelativity.cs b/Numerics/Numerics/Physics/Relativity/SpecialRelativity.cs new file mode 100644 index 0000000..f843ef0 --- /dev/null +++ b/Numerics/Numerics/Physics/Relativity/SpecialRelativity.cs @@ -0,0 +1,150 @@ +using System; +using CSharpNumerics.Physics.Constants; + +namespace CSharpNumerics.Physics.Relativity; + +/// +/// Special-relativistic kinematics and dynamics in flat spacetime. +/// All speeds are in m/s and must be strictly below the speed of light. +/// These quantities are the foundation for the general-relativistic models +/// (Schwarzschild geometry, GPS-style clock corrections, etc.). +/// +public static class SpecialRelativity +{ + private const double C = PhysicsConstants.SpeedOfLight; + private const double C2 = C * C; + + /// + /// Lorentz factor γ = 1 / √(1 − v²/c²). + /// + /// Speed in m/s (sign is irrelevant). + /// γ ≥ 1. + /// If |v| ≥ c. + public static double LorentzFactor(double velocity) + { + double beta2 = velocity * velocity / C2; + if (beta2 >= 1.0) + throw new ArgumentOutOfRangeException(nameof(velocity), "Speed must be below the speed of light."); + + return 1.0 / Math.Sqrt(1.0 - beta2); + } + + /// + /// Time dilation: the coordinate time Δt elapsed in the rest frame while a + /// clock moving at measures proper time + /// . Δt = γ·Δτ (the moving clock runs slow). + /// + public static double TimeDilation(double properTime, double velocity) + { + return LorentzFactor(velocity) * properTime; + } + + /// + /// Length contraction: the length L = L₀/γ measured for an object whose + /// proper (rest-frame) length is when it + /// moves at along the measured direction. + /// + public static double LengthContraction(double properLength, double velocity) + { + return properLength / LorentzFactor(velocity); + } + + /// + /// Relativistic momentum p = γ·m·v. + /// + public static double Momentum(double mass, double velocity) + { + return LorentzFactor(velocity) * mass * velocity; + } + + /// + /// Rest energy E₀ = m·c². + /// + public static double RestEnergy(double mass) + { + return mass * C2; + } + + /// + /// Total relativistic energy E = γ·m·c². + /// + public static double TotalEnergy(double mass, double velocity) + { + return LorentzFactor(velocity) * mass * C2; + } + + /// + /// Relativistic kinetic energy E_k = (γ − 1)·m·c². + /// Evaluated as m·c²·β²/(s·(1 + s)) with s = √(1 − β²), which is + /// algebraically identical to (γ − 1)·m·c² but avoids the catastrophic + /// cancellation of computing γ − 1 directly at low speeds (where it + /// correctly reduces to the classical ½·m·v²). + /// + public static double KineticEnergy(double mass, double velocity) + { + double beta2 = velocity * velocity / C2; + if (beta2 >= 1.0) + throw new ArgumentOutOfRangeException(nameof(velocity), "Speed must be below the speed of light."); + + double s = Math.Sqrt(1.0 - beta2); + return mass * C2 * beta2 / (s * (1.0 + s)); + } + + /// + /// Total energy from the energy–momentum relation E = √((p·c)² + (m·c²)²). + /// + public static double EnergyFromMomentum(double mass, double momentum) + { + double pc = momentum * C; + double mc2 = mass * C2; + return Math.Sqrt(pc * pc + mc2 * mc2); + } + + /// + /// Relativistic velocity addition: combines two collinear velocities + /// s = (u + v) / (1 + u·v/c²). The result is always below c when both + /// inputs are. + /// + public static double AddVelocities(double u, double v) + { + return (u + v) / (1.0 + u * v / C2); + } + + /// + /// Relativistic longitudinal Doppler factor f_observed / f_source for purely + /// radial motion: √((1 − β)/(1 + β)), with β = v/c. + /// + /// Radial speed in m/s. Positive = receding (redshift), + /// negative = approaching (blueshift). + public static double DopplerFactor(double velocity) + { + double beta = velocity / C; + if (Math.Abs(beta) >= 1.0) + throw new ArgumentOutOfRangeException(nameof(velocity), "Speed must be below the speed of light."); + + return Math.Sqrt((1.0 - beta) / (1.0 + beta)); + } + + /// + /// Observed frequency under the relativistic longitudinal Doppler effect. + /// + /// Emitted frequency (Hz). + /// Radial speed in m/s. Positive = receding. + public static double ObservedFrequency(double sourceFrequency, double velocity) + { + return sourceFrequency * DopplerFactor(velocity); + } + + /// + /// Rapidity φ = artanh(v/c), the additive measure of velocity + /// (rapidities add linearly under collinear boosts). + /// + public static double Rapidity(double velocity) + { + double beta = velocity / C; + if (Math.Abs(beta) >= 1.0) + throw new ArgumentOutOfRangeException(nameof(velocity), "Speed must be below the speed of light."); + + return 0.5 * Math.Log((1.0 + beta) / (1.0 - beta)); + } +}