Fix/srp atm (#274)

* Refactor spacecraft parameters and add shadow fraction calculations for celestial items

* Update spacecraft documentation and versioning for Astrodynamics framework

* Refactor GeopotentialModelReader and GravitationalField initialization to use FileStream for improved file access

Author: sylvain-guillet

DEVELOPER_GUIDE.md

@@ -1226,15 +1226,18 @@ Represents a spacecraft with components.
 
 | Constructor | Description |
 |------------|-------------|
-| `Spacecraft(int naifId, string name, double dryMass, double maximumThrustPower, Clock clock, OrbitalParameters orbit)` | Create spacecraft |
+| `Spacecraft(int naifId, string name, double dryMass, double maxMass, Clock clock, OrbitalParameters orbit, double sectionalArea = 1.0, double dragCoeff = 2.2, string cosparId = null, double solarRadiationCoeff = 1.0)` | Create spacecraft |
 
 | Property | Description |
 |----------|-------------|
 | `NaifId` | NAIF ID (negative) |
 | `Name` | Spacecraft name |
-| `DryMass` | Mass without fuel (kg) |
-| `MaximumThrustPower` | Maximum thrust power (W) |
+| `DryOperatingMass` | Mass without fuel (kg) |
+| `MaximumOperatingMass` | Maximum operating mass (kg) |
 | `Clock` | Onboard clock |
+| `SectionalArea` | Mean cross-sectional area for drag and SRP (m²) |
+| `DragCoefficient` | Drag coefficient Cd (default 2.2) |
+| `SolarRadiationCoeff` | Solar radiation pressure coefficient Cr (default 1.0, range 1.0–2.0) |
 | `InitialOrbitalParameters` | Initial orbit |
 | `StateVectorsRelativeToICRF` | Propagated states |
 
@@ -1467,8 +1470,8 @@ Numerical orbit propagator using a Velocity-Verlet (symplectic) integrator with
 
 **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
+- **Atmospheric drag** (when `drag` is true): uses the body's atmospheric model with atmosphere-relative velocity (accounts for body co-rotation). Default Cd = 2.2 (free-molecular flow). Mass ratio is dynamic (tracks fuel consumption via `GetTotalMass()`).
+- **Solar radiation pressure** (when `srp` is true): cannonball model with reflectivity coefficient Cr (`Spacecraft.SolarRadiationCoeff`, default 1.0). Uses continuous shadow fraction for partial/annular eclipse geometry instead of binary eclipse detection. Mass ratio is dynamic.
 
 **Geopotential usage:**
 

IO.Astrodynamics.Net/CLAUDE.md

@@ -432,6 +432,45 @@ var propagator = new SpacecraftPropagator(window, spacecraft,
 
 **Thread Safety:** `GeopotentialGravitationalField` is NOT thread-safe (pre-allocated buffers). Create separate `CelestialBody` instances for concurrent propagation.
 
+### Force Models (Propagator Perturbations)
+
+The `IO.Astrodynamics.Propagator.Forces` namespace provides configurable perturbation models used by the `SpacecraftPropagator`.
+
+**Atmospheric Drag (`AtmosphericDrag`)**
+- Uses **atmosphere-relative velocity**: body-centered velocity minus co-rotation (`v_rel = v_body - omega x r_body`)
+- The body's angular velocity is obtained via `CelestialBody.GetOrientation(Frame.ICRF, epoch).AngularVelocity`
+- **Dynamic mass**: area/mass ratio is computed per step using `Spacecraft.GetTotalMass()` (dry + fuel + payload)
+- Drag coefficient default is **2.2** (appropriate for satellites in free-molecular flow)
+
+```csharp
+// Atmospheric drag requires a CelestialBody with an atmospheric model
+var earth = new CelestialBody(PlanetsAndMoons.EARTH, Frames.Frame.ICRF, epoch,
+    geopotentialParams, new EarthStandardAtmosphere());
+var drag = new AtmosphericDrag(spacecraft, earth);
+```
+
+**Solar Radiation Pressure (`SolarRadiationPressure`)**
+- Cannonball model: `F = (L_sun / 4*pi*c) * Cr * (A/m) * r_hat / r^2`
+- Includes **reflectivity coefficient Cr** from `Spacecraft.SolarRadiationCoeff` (range 1.0–2.0, default 1.0)
+- **Continuous shadow model**: uses `ShadowFraction` (two-circle overlap area) instead of binary eclipse check
+  - Handles none, partial, annular, and total eclipse geometries
+  - SRP is scaled by `(1 - maxShadowFraction)` across all occluding bodies
+- **Dynamic mass**: area/mass ratio recomputed each step via `GetTotalMass()`
+- Occluding bodies are materialized to `CelestialBody[]` (avoids IEnumerable re-evaluation)
+
+```csharp
+// SRP with Cr = 1.5
+var spacecraft = new Spacecraft(-1001, "Sat", 100.0, 10000.0, clock, orbit,
+    sectionalArea: 10.0, solarRadiationCoeff: 1.5);
+var srp = new SolarRadiationPressure(spacecraft, [earth]);
+```
+
+**Shadow Fraction (`CelestialItem.ShadowFraction`)**
+- Static method: `ShadowFraction(double angularSeparation, double backSize, double bySize)` → `[0, 1]`
+- Instance method: `ShadowFraction(CelestialItem by, OrbitalParameters from, Aberration aberration)` → `[0, 1]`
+- Returns 0.0 for full illumination, 1.0 for total eclipse
+- Uses the standard two-circle intersection area formula for partial eclipses
+
 ### Attitude Maneuvers
 
 The framework provides a family of attitude maneuvers for spacecraft orientation control. All inherit from the abstract `Attitude` base class.
@@ -574,6 +613,13 @@ Test data files are in `Data/SolarSystem/` and copied to output directory.
    - `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
+11. **Force Models (Drag & SRP)**: When working with atmospheric drag or solar radiation pressure:
+   - Default drag coefficient (Cd) is **2.2** (free-molecular flow for satellites), not 0.3
+   - Atmospheric drag uses **atmosphere-relative velocity** (body-centered velocity minus co-rotation)
+   - SRP includes **Cr coefficient** from `Spacecraft.SolarRadiationCoeff` (default 1.0, range 1.0–2.0)
+   - SRP uses **continuous shadow fraction** (partial/annular eclipses reduce SRP proportionally)
+   - Both forces use **dynamic mass** via `GetTotalMass()` (accounts for fuel consumption during propagation)
+   - Use `CelestialItem.ShadowFraction()` for eclipse geometry calculations
 
 ## 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.1</Version>
+        <Version>0.8.3.2</Version>
         <Title>Astrodynamics command line interface</Title>
         <Authors>Sylvain Guillet</Authors>
         <Description>This CLI allows end user to exploit IO.Astrodynamics framework </Description>

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/Body/CelestialBodyTests.cs

@@ -333,7 +333,7 @@ public void PhaseHelioSynchronousOrbit()
     public void GeopotentialModelReader()
     {
         GeopotentialModelReader geopotentialModelReader =
-            new GeopotentialModelReader(new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+            new GeopotentialModelReader(new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
         Assert.Equal(new GeopotentialCoefficient(4, 1, -0.536157389388867E-06, -0.473567346518086E-06, 0.4568074333E-11, 0.4684043490E-11),
             geopotentialModelReader.ReadCoefficient(4, 1));
     }
@@ -342,7 +342,7 @@ public void GeopotentialModelReader()
     public void GeopotentialModelReaderException()
     {
         GeopotentialModelReader geopotentialModelReader =
-            new GeopotentialModelReader(new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+            new GeopotentialModelReader(new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
         Assert.Throws<ArgumentException>(() => geopotentialModelReader.ReadCoefficient(4, 5));
     }
 
@@ -372,6 +372,46 @@ public void IsOccultedAnnular()
         Assert.Equal(OccultationType.Annular, CelestialItem.IsOcculted(1.0, 4.0, 2.0));
     }
 
+    [Fact]
+    public void ShadowFractionNone()
+    {
+        // No occultation: angular separation larger than sum of radii
+        Assert.Equal(0.0, CelestialItem.ShadowFraction(3.0, 2.0, 4.0));
+        Assert.Equal(0.0, CelestialItem.ShadowFraction(3.0, 4.0, 2.0));
+    }
+
+    [Fact]
+    public void ShadowFractionFull()
+    {
+        // Total eclipse: occluding body completely covers the light source
+        Assert.Equal(1.0, CelestialItem.ShadowFraction(0.5, 2.0, 4.0));
+    }
+
+    [Fact]
+    public void ShadowFractionAnnular()
+    {
+        // Annular eclipse: occluding body fully inside the sun disc
+        // rOcc = 1.0, rSun = 2.0, fraction = (1.0^2)/(2.0^2) = 0.25
+        Assert.Equal(0.25, CelestialItem.ShadowFraction(0.5, 4.0, 2.0));
+    }
+
+    [Fact]
+    public void ShadowFractionPartial()
+    {
+        // Partial eclipse: between no occultation and full/annular
+        var fraction = CelestialItem.ShadowFraction(2.0, 2.0, 4.0);
+        Assert.True(fraction > 0.0);
+        Assert.True(fraction < 1.0);
+    }
+
+    [Fact]
+    public void ShadowFractionEdge()
+    {
+        // Edge case: exactly touching (d = rSun + rOcc)
+        // rSun = 1.0, rOcc = 2.0, d = 3.0 -> no occultation
+        Assert.Equal(0.0, CelestialItem.ShadowFraction(3.0, 2.0, 4.0));
+    }
+
     [Fact]
     public void EarthAirTemperature()
     {

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/Body/GravitationalAccelerationTest.cs

@@ -18,7 +18,7 @@ public GravitationalAccelerationTest()
     public void ComputeGeopotentialGravityAcceleration()
     {
         GeopotentialGravitationalField gravity =
-            new GeopotentialGravitationalField(new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+            new GeopotentialGravitationalField(new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
         StateVector parkingOrbit = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), TestHelpers.EarthWithAtmAndGeoAtJ2000,
             TimeSystem.Time.J2000TDB,
             Frames.Frame.ICRF);
@@ -47,7 +47,7 @@ public void ComputeGeopotentialAt45DegreesLatitude()
         // Position at ~45° latitude exercises tesseral harmonic terms
         GeopotentialGravitationalField gravity =
             new GeopotentialGravitationalField(
-                new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+                new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
 
         // r=6800km, 45° latitude: x=r*cos(45°), y=0, z=r*sin(45°)
         double r = 6800000.0;
@@ -77,7 +77,7 @@ public void ComputeGeopotentialSouthernHemisphere()
         // Southern hemisphere test: odd-degree zonal harmonics (J3, J5) should cause asymmetry
         GeopotentialGravitationalField gravity =
             new GeopotentialGravitationalField(
-                new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+                new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
 
         double r = 6800000.0;
         double component = r / System.Math.Sqrt(2.0);
@@ -117,7 +117,7 @@ public void ComputeGeopotentialWithLongitude()
         // Position with longitude exercises S_nm sine terms
         GeopotentialGravitationalField gravity =
             new GeopotentialGravitationalField(
-                new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+                new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)));
 
         double r = 6800000.0;
         double component = r / System.Math.Sqrt(2.0);
@@ -147,7 +147,7 @@ public void ComputeJ2OnlyAnalytical()
         // Use degree 2 only, which is dominated by J2 (C20 term)
         GeopotentialGravitationalField gravity =
             new GeopotentialGravitationalField(
-                new StreamReader(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")),
+                new StreamReader(new FileStream(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree"), FileMode.Open, FileAccess.Read, FileShare.Read)),
                 maxDegrees: 2);
 
         // Position in body-fixed frame at equator (φ=0)

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/CCSDS/OPM/OpmTests.cs

@@ -787,7 +787,7 @@ public void ToSpacecraft_WithMinimalOpm_CreatesSpacecraftWithDefaults()
             Assert.Equal(1.0, spacecraft.Mass); // Default mass when not specified
             Assert.Equal(500000.0, spacecraft.MaximumOperatingMass);
             Assert.Equal(1.0, spacecraft.SectionalArea); // Default
-            Assert.Equal(0.3, spacecraft.DragCoefficient); // Default
+            Assert.Equal(2.2, spacecraft.DragCoefficient); // Default
             Assert.Equal(1.0, spacecraft.SolarRadiationCoeff); // Default
             Assert.Same(clock, spacecraft.Clock);
         }

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/Propagators/Integrators/Forces/AtmosphericDragTests.cs

@@ -34,6 +34,6 @@ public void ComputeAcceleration()
         StateVector parkingOrbit = new StateVector(new Vector3(7380000.0, 0.0, 0.0), new Vector3(0.0, 9700.0, 0.0), earth, TimeSystem.Time.J2000TDB,
             Frames.Frame.ICRF);
         var res = atmosphericDrag.Apply(parkingOrbit);
-        Assert.Equal(new Vector3(0.0, -1.4911628830275622E-09, 0.0), res);
+        Assert.Equal(new Vector3(-0.0, -1.3302926867698408E-09, 2.184861621946448E-15), res, TestHelpers.VelocityVectorComparer);
     }
 }

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/Propagators/Integrators/Forces/SolarRadiationPressureTests.cs

@@ -1,8 +1,10 @@
-using System.IO;
+using System;
+using System.IO;
 using IO.Astrodynamics.Body;
 using IO.Astrodynamics.Body.Spacecraft;
 using IO.Astrodynamics.OrbitalParameters;
 using IO.Astrodynamics.Propagator.Forces;
+using IO.Astrodynamics.SolarSystemObjects;
 using IO.Astrodynamics.TimeSystem;
 using Xunit;
 using CelestialBody = IO.Astrodynamics.Body.CelestialBody;
@@ -24,11 +26,73 @@ public void ComputeAcceleration()
         Clock clk = new Clock("My clock", 256);
         Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk,
             new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB, Frames.Frame.ICRF));
-        SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure(spc,[earth]);
+        SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure(spc, [earth]);
 
         StateVector parkingOrbit = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB,
             Frames.Frame.ICRF);
         var res = solarRadiationPressure.Apply(parkingOrbit);
-         Assert.Equal(new Vector3(-8.456677063948622E-09, 4.237795744348693E-08, 1.8372880484798083E-08), res, TestHelpers.VectorComparer);
+        Assert.Equal(new Vector3(-8.456677063948622E-09, 4.237795744348693E-08, 1.8372880484798083E-08), res, TestHelpers.VectorComparer);
     }
-}
+
+    [Fact]
+    public void ComputeAccelerationWithCr()
+    {
+        // Cr = 1.8 should scale the SRP by 1.8x compared to Cr = 1.0
+        var earth = new CelestialBody(399, new GeopotentialModelParameters(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+        Clock clk = new Clock("My clock", 256);
+        Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk,
+            new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB, Frames.Frame.ICRF),
+            solarRadiationCoeff: 1.8);
+        SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure(spc, [earth]);
+
+        StateVector parkingOrbit = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB,
+            Frames.Frame.ICRF);
+        var res = solarRadiationPressure.Apply(parkingOrbit);
+
+        // Expected = base result * 1.8
+        var expected = new Vector3(-8.456677063948622E-09 * 1.8, 4.237795744348693E-08 * 1.8, 1.8372880484798083E-08 * 1.8);
+        Assert.Equal(expected, res, TestHelpers.VectorComparer);
+    }
+
+    [Fact]
+    public void ComputeAccelerationWithDynamicMass()
+    {
+        // Spacecraft with fuel tank — GetTotalMass() = dry mass + fuel = 100 + 50 = 150 kg
+        // areaMassRatio = 1.0 / 150 instead of 1.0 / 100
+        // So result should be 100/150 = 2/3 of the base result
+        var earth = new CelestialBody(399, new GeopotentialModelParameters(Path.Combine(Constants.SolarSystemKernelPath.ToString(), "EGM2008_to70_TideFree")));
+        Clock clk = new Clock("My clock", 256);
+        Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk,
+            new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB, Frames.Frame.ICRF));
+        var fuelTank = new FuelTank("tank1", "model1", "sn001", 100.0, 50.0);
+        spc.AddFuelTank(fuelTank);
+
+        SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure(spc, [earth]);
+
+        StateVector parkingOrbit = new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB,
+            Frames.Frame.ICRF);
+        var res = solarRadiationPressure.Apply(parkingOrbit);
+
+        // Expected = base result * (100/150)
+        double ratio = 100.0 / 150.0;
+        var expected = new Vector3(-8.456677063948622E-09 * ratio, 4.237795744348693E-08 * ratio, 1.8372880484798083E-08 * ratio);
+        Assert.Equal(expected, res, TestHelpers.VectorComparer);
+    }
+
+    [Fact]
+    public void ConstructorThrowsOnNullSpacecraft()
+    {
+        var earth = new CelestialBody(399);
+        Assert.Throws<ArgumentNullException>(() => new SolarRadiationPressure(null, [earth]));
+    }
+
+    [Fact]
+    public void ConstructorThrowsOnNullOccultingBodies()
+    {
+        var earth = new CelestialBody(399);
+        Clock clk = new Clock("My clock", 256);
+        Spacecraft spc = new Spacecraft(-1001, "MySpacecraft", 100.0, 10000.0, clk,
+            new StateVector(new Vector3(6800000.0, 0.0, 0.0), new Vector3(0.0, 7656.2204182967143, 0.0), earth, TimeSystem.Time.J2000TDB, Frames.Frame.ICRF));
+        Assert.Throws<ArgumentNullException>(() => new SolarRadiationPressure(spc, null));
+    }
+}