Fix/harmonics (#272)

* Refactor GeopotentialGravitationalField for improved performance and accuracy in gravitational calculations

* Add test for 24h propagation with degree-10 geopotential and validate position error

* Enhance GeopotentialGravitationalField with detailed comments and improve Legendre function handling at poles; add new WebFetch domains in settings

Author: sylvain-guillet

DEVELOPER_GUIDE.md

@@ -65,6 +65,7 @@ public const double Rad2Deg = 180.0 / Math.PI; // Radians to degrees
 - **Time**: IERS Conventions, leap seconds from NAIF LSK kernels
 - **Ephemerides**: JPL Development Ephemeris (DE series)
 - **Reference Frames**: IAU/IAG standards, IERS conventions
+- **Gravity Models**: EGM2008 (Earth Gravitational Model 2008), geodesy-normalized spherical harmonics up to degree/order 70
 - **Atmospheric Models**: U.S. Standard Atmosphere 1976, NRLMSISE-00
 - **TLE Propagation**: SGP4/SDP4 (AFSPC compatibility)
 
@@ -302,6 +303,7 @@ Represents a natural celestial body (planet, moon, star).
 | `Flattening` | Flattening coefficient |
 | `Frame` | Body-fixed reference frame |
 | `SphereOfInfluence` | Sphere of influence radius (m) |
+| `GravitationalField` | Gravity model (point-mass or geopotential) |
 | `InitialOrbitalParameters` | Initial orbital state |
 
 | Method | Description |
@@ -1453,7 +1455,7 @@ Console.WriteLine($"Delta-V at arrival: {solution.ArrivalVelocity.Magnitude():F1
 
 #### SpacecraftPropagator
 
-Numerical orbit propagator with perturbations.
+Numerical orbit propagator using a Velocity-Verlet (symplectic) integrator with configurable perturbation models.
 
 | Constructor | Description |
 |------------|-------------|
@@ -1463,6 +1465,37 @@ Numerical orbit propagator with perturbations.
 |--------|-------------|
 | `Propagate()` | Execute propagation |
 
+**Force Models:**
+- **Gravitational acceleration** from each body in the `bodies` list. If a body has a `GeopotentialModelParameters`, the full spherical harmonic model (EGM2008) is used; otherwise point-mass gravity is applied.
+- **Atmospheric drag** (when `drag` is true): uses the body's atmospheric model
+- **Solar radiation pressure** (when `srp` is true): cannonball model with eclipse detection
+
+**Geopotential usage:**
+
+```csharp
+// Create Earth with degree-10 EGM2008 geopotential
+var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    new GeopotentialModelParameters("Data/SolarSystem/EGM2008_to70_TideFree", 10));
+
+// Initial state in ICRF
+var orbit = new StateVector(
+    new Vector3(6800000.0, 0, 0),
+    new Vector3(0, 7500.0, 0),
+    earth, epoch, Frames.Frame.ICRF);
+
+var spacecraft = new Spacecraft(-1001, "Sat", 100.0, 10000.0, clock, orbit);
+
+// Propagate with geopotential, Moon, and Sun perturbations
+var propagator = new SpacecraftPropagator(
+    new Window(epoch, epoch.AddDays(1)),
+    spacecraft,
+    [earth, PlanetsAndMoons.MOON_BODY, Stars.SUN_BODY],
+    false, false, TimeSpan.FromSeconds(1.0));
+propagator.Propagate();
+```
+
+**Accuracy:** With degree-10 EGM2008, Moon, and Sun perturbations, 24h LEO propagation achieves sub-kilometer position accuracy compared to STK HPOP reference solutions.
+
 #### TLEPropagator
 
 SGP4/SDP4 propagator for TLE elements.
@@ -1477,6 +1510,74 @@ SGP4/SDP4 propagator for TLE elements.
 
 ---
 
+### Geopotential Gravity Models
+
+The framework supports high-fidelity Earth gravity modeling using spherical harmonic expansion with the EGM2008 model (up to degree/order 70).
+
+#### GeopotentialModelParameters
+
+Configuration for enabling geopotential gravity on a `CelestialBody`.
+
+| Constructor | Description |
+|------------|-------------|
+| `GeopotentialModelParameters(string path, ushort degree = 60)` | Load model from file path |
+| `GeopotentialModelParameters(Stream stream, ushort degree = 60)` | Load model from stream |
+
+| Property | Description |
+|----------|-------------|
+| `GeopotentialModelPath` | Stream reader for the model file |
+| `GeopotentialDegree` | Maximum degree/order to use |
+
+#### GeopotentialGravitationalField
+
+Computes the full 3D gravitational acceleration including spherical harmonic perturbations using the Montenbruck & Gill spherical coordinate gradient formulation.
+
+| Property | Description |
+|----------|-------------|
+| `MaxDegree` | Maximum harmonic degree |
+
+| Method | Description |
+|--------|-------------|
+| `ComputeGravitationalAcceleration(StateVector sv)` | Compute full 3D acceleration including harmonics |
+
+**Algorithm:**
+1. Transform position to body-fixed frame
+2. Compute geocentric latitude, longitude, and radius
+3. Evaluate geodesy-normalized associated Legendre functions and derivatives via stable recursion
+4. Sum radial, latitudinal, and longitudinal gradient components over all harmonic degrees/orders
+5. Transform spherical acceleration to body-fixed Cartesian, then rotate to inertial frame
+
+**Normalization convention:** Geodesy (fully-normalized) without Condon-Shortley phase:
+`P_bar_nm = sqrt((2 - delta_0m)(2n+1)(n-m)!/(n+m)!) * P_nm`
+
+**Supported model file:** EGM2008 (tide-free), included as `Data/SolarSystem/EGM2008_to70_TideFree`. The file provides coefficients C_nm and S_nm from degree 2 to degree 70.
+
+**Thread safety:** `GeopotentialGravitationalField` is **not** thread-safe. Pre-allocated buffers (Legendre tables, trig arrays) are reused across calls. Create a separate `CelestialBody` instance per thread when propagating concurrently.
+
+**Degree selection guidelines:**
+
+| Degree | Use Case | Typical LEO Accuracy (24h) |
+|--------|----------|---------------------------|
+| 2 | J2-only, fast preliminary analysis | ~10 km |
+| 10 | Good accuracy for most LEO missions | < 1 km vs STK HPOP |
+| 20-30 | High-fidelity geodesy applications | Sub-100 m |
+| 70 | Maximum available (EGM2008_to70) | Highest fidelity |
+
+```csharp
+// Degree-10 geopotential (good balance of accuracy and speed)
+var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    new GeopotentialModelParameters("Data/SolarSystem/EGM2008_to70_TideFree", 10));
+
+// Degree-2 (J2-only) for quick analysis
+var earthJ2 = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    new GeopotentialModelParameters("Data/SolarSystem/EGM2008_to70_TideFree", 2));
+
+// Point-mass only (no geopotential)
+var earthPointMass = new CelestialBody(PlanetsAndMoons.EARTH);
+```
+
+---
+
 ### IO.Astrodynamics.Atmosphere
 
 The atmosphere subsystem provides a unified interface for atmospheric modeling across different planets and model complexities.
@@ -2455,6 +2556,43 @@ Console.WriteLine($"  Position error: {positionError:F1} m");
 Console.WriteLine($"  Velocity error: {velocityError:F4} m/s");
 ```
 
+### High-Fidelity LEO Propagation with Geopotential
+
+Propagate a LEO satellite with EGM2008 spherical harmonics, Moon, and Sun perturbations:
+
+```csharp
+var epoch = new Time(2025, 8, 25, 11, 55, 44, frame: TimeFrame.UTCFrame);
+
+// Earth with degree-10 EGM2008 geopotential
+var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    new GeopotentialModelParameters("Data/SolarSystem/EGM2008_to70_TideFree", 10));
+
+// LEO orbit (~400 km altitude)
+var orbit = new StateVector(
+    new Vector3(5442162.59, -4068949.85, -13456.85),   // Position (m)
+    new Vector3(2858.20, 3809.79, 6002.13),             // Velocity (m/s)
+    earth, epoch, Frames.Frame.ICRF);
+
+var clock = new Clock("Clock", 256);
+var spacecraft = new Spacecraft(-1001, "LEOSat", 100.0, 10000.0, clock, orbit);
+
+var propagator = new SpacecraftPropagator(
+    new Window(epoch, epoch.AddDays(1)),
+    spacecraft,
+    [earth, PlanetsAndMoons.MOON_BODY, Stars.SUN_BODY],
+    false, false,
+    TimeSpan.FromSeconds(1.0));
+
+propagator.Propagate();
+
+// Get final state relative to Earth
+var finalState = spacecraft.StateVectorsRelativeToICRF.Values.Last()
+    .RelativeTo(earth, Aberration.None) as StateVector;
+
+Console.WriteLine($"Position after 24h: ({finalState.Position.X/1000:F3}, " +
+    $"{finalState.Position.Y/1000:F3}, {finalState.Position.Z/1000:F3}) km");
+```
+
 ### Earth Observation with Attitude Maneuvers
 
 Point instruments at ground targets during observation passes:

IO.Astrodynamics.Net/.claude/settings.local.json

@@ -31,7 +31,10 @@
       "WebFetch(domain:en.wikipedia.org)",
       "WebFetch(domain:johannesbuchner.github.io)",
       "WebFetch(domain:www.scratchapixel.com)",
-      "Bash(tee:*)"
+      "Bash(tee:*)",
+      "Bash(xargs cat:*)",
+      "WebFetch(domain:shtools.github.io)",
+      "WebFetch(domain:degenerateconic.com)"
     ],
     "deny": [],
     "ask": []

IO.Astrodynamics.Net/CLAUDE.md

@@ -71,8 +71,11 @@ dotnet tool install --global --add-source ./IO.Astrodynamics.CLI/bin/Debug IO.As
 - `IO.Astrodynamics.Maneuver`: Lambert solvers, launch windows, maneuver planning, attitude maneuvers
 - `IO.Astrodynamics.Frames`: Reference frames and transformations
 - `IO.Astrodynamics.TimeSystem`: Time frames (UTC, TDB, TAI, etc.)
-- `IO.Astrodynamics.Propagator`: Orbital propagation and integration
+- `IO.Astrodynamics.Propagator`: Orbital propagation and integration (Velocity-Verlet symplectic integrator)
+- `IO.Astrodynamics.Propagator.Forces`: Force models (gravitational, atmospheric drag, solar radiation pressure)
 - `IO.Astrodynamics.Atmosphere`: Atmospheric density, temperature, and pressure models for Earth and Mars
+- `IO.Astrodynamics.Math`: Vectors, matrices, quaternions, Legendre functions
+- `IO.Astrodynamics.Physics`: Geopotential model reader, coefficients
 
 **External Dependencies**
 - MathNet.Numerics: Linear algebra operations
@@ -399,6 +402,36 @@ var atm = earth.GetAtmosphere(context);  // Uses NRLMSISE-00 automatically
 - Use `CelestialBody.GetAtmosphere(context)` with full context for automatic NRLMSISE-00 selection
 - Use `MarsStandardAtmosphere` for preliminary Mars mission analysis
 
+### Geopotential Gravity Model
+
+The framework supports spherical harmonic gravity modeling using EGM2008 coefficients (up to degree/order 70).
+
+**Key Classes**
+- `GeopotentialModelParameters`: Configuration (model file path + max degree)
+- `GeopotentialGravitationalField`: Computes full 3D acceleration via Montenbruck & Gill formulation
+- `GeopotentialModelReader`: Parses EGM2008 coefficient files
+- `GeopotentialCoefficient`: Holds C_nm, S_nm coefficients for a single (n,m) pair
+- `LegendreFunctions`: Geodesy-normalized associated Legendre functions and derivatives
+
+**Conventions**
+- Geodesy normalization: `sqrt((2-delta_0m)(2n+1)(n-m)!/(n+m)!)` — no Condon-Shortley phase
+- Coefficients: fully-normalized C_nm and S_nm from EGM2008 (tide-free)
+- Model file starts at degree 2 (degrees 0 and 1 are absent; handled as zeros)
+
+**Usage**
+```csharp
+// Create Earth with degree-10 geopotential
+var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    new GeopotentialModelParameters("Data/SolarSystem/EGM2008_to70_TideFree", 10));
+
+// The propagator automatically uses the geopotential model when present
+var propagator = new SpacecraftPropagator(window, spacecraft,
+    [earth, PlanetsAndMoons.MOON_BODY, Stars.SUN_BODY],
+    false, false, TimeSpan.FromSeconds(1.0));
+```
+
+**Thread Safety:** `GeopotentialGravitationalField` is NOT thread-safe (pre-allocated buffers). Create separate `CelestialBody` instances for concurrent propagation.
+
 ### Attitude Maneuvers
 
 The framework provides a family of attitude maneuvers for spacecraft orientation control. All inherit from the abstract `Attitude` base class.
@@ -534,6 +567,13 @@ Test data files are in `Data/SolarSystem/` and copied to output directory.
    - Ensure body vectors and reference vectors are not collinear (minimum 5 degrees separation)
    - Use `Spacecraft.Front`, `Spacecraft.Up`, etc. for standard body frame directions
    - Use `Instrument.GetBoresightInSpacecraftFrame()` and `GetRefVectorInSpacecraftFrame()` for instrument-based pointing
+10. **Geopotential Gravity**: When working with spherical harmonic gravity models:
+   - Pass `GeopotentialModelParameters` to `CelestialBody` constructor to enable geopotential
+   - Use degree 10 for typical LEO accuracy; degree 2 for J2-only analysis
+   - EGM2008 model file: `Data/SolarSystem/EGM2008_to70_TideFree` (degrees 2-70)
+   - `GeopotentialGravitationalField` is NOT thread-safe — one instance per thread
+   - Legendre functions use geodesy normalization without Condon-Shortley phase
+   - The model file starts at degree 2; degrees 0 and 1 are handled as zeros internally
 
 ## Code Quality Standards
 

IO.Astrodynamics.Net/IO.Astrodynamics.CLI/IO.Astrodynamics.CLI.csproj

@@ -9,7 +9,7 @@
         <FileVersion>0.0.1</FileVersion>
         <PackAsTool>true</PackAsTool>
         <ToolCommandName>astro</ToolCommandName>
-        <Version>0.8.3.0</Version>
+        <Version>0.8.3.1</Version>
         <Title>Astrodynamics command line interface</Title>
         <Authors>Sylvain Guillet</Authors>
         <Description>This CLI allows end user to exploit IO.Astrodynamics framework </Description>