Strenghten atmospheric drag constructor. (#180)

All bodies in scenario can be involved in atmospheric drag and solar radiation pressure

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.2.1</Version>
+        <Version>0.2.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.Performance/Scenario.cs

@@ -42,14 +42,14 @@ public Scenario()
         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, DateTimeExtension.J2000, Frames.Frame.ICRF));
-        _srp = new SolarRadiationPressure(spc);
+        _srp = new SolarRadiationPressure(spc,[_earth]);
         _atm = new AtmosphericDrag(spc, _earth);
         List<ForceBase> forces = new List<ForceBase>();
         forces.Add(new GravitationalAcceleration(_sun));
         forces.Add(new GravitationalAcceleration(_moon));
         forces.Add(new GravitationalAcceleration(_earth));
         forces.Add(new AtmosphericDrag(spc, _earth));
-        forces.Add(new SolarRadiationPressure(spc));
+        forces.Add(new SolarRadiationPressure(spc,[_earth]));
         _integrator = new VVIntegrator(forces, TimeSpan.FromSeconds(1.0), new StateVector(new Vector3(6800000.0 - Random.Shared.NextDouble(), 0.0, 0.0),
             new Vector3(0.0, 8000.0 - Random.Shared.NextDouble(), 0.0), _earth,
             DateTimeExtension.J2000, Frame.ICRF));

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

@@ -24,7 +24,7 @@ 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, DateTimeExtension.J2000, Frames.Frame.ICRF));
-        SolarRadiationPressure solarRadiationPressure = new SolarRadiationPressure(spc);
+        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, DateTimeExtension.J2000,
             Frames.Frame.ICRF);

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/Propagators/Integrators/IntegratorTests.cs

@@ -34,7 +34,7 @@ public void IntegrateWithPerturbations()
         forces.Add(new GravitationalAcceleration(moon));
         forces.Add(new GravitationalAcceleration(earth));
         forces.Add(new AtmosphericDrag(spc, earth));
-        forces.Add(new SolarRadiationPressure(spc));
+        forces.Add(new SolarRadiationPressure(spc,[earth]));
         VVIntegrator vvIntegrator = new VVIntegrator(forces, TimeSpan.FromSeconds(1.0), spc.InitialOrbitalParameters.ToStateVector());
         StateVector[] data = new StateVector[2];
         Array.Fill(data, spc.InitialOrbitalParameters.ToStateVector(), 0, 2);

IO.Astrodynamics.Net/IO.Astrodynamics/Body/CelestialBody.cs

@@ -24,8 +24,9 @@ public class CelestialBody : CelestialItem, IOrientable
     public double J3 { get; }
     public double J4 { get; }
 
-
     protected AtmosphericModel AtmosphericModel { get; }
+
+    public bool HasAtmosphericModel => AtmosphericModel != null;
 
     /// <summary>
     /// Instantiate celestial body from naif object with default parameters (Ecliptic J2000 at J2000 epoch)

IO.Astrodynamics.Net/IO.Astrodynamics/IO.Astrodynamics.nuspec

@@ -4,7 +4,7 @@
         <id>IO.Astrodynamics</id>
         <authors>Sylvain Guillet</authors>
         <copyright>Sylvain Guillet</copyright>
-        <version>3.3.1</version>
+        <version>3.4.0</version>
         <title>Astrodynamics framework</title>
         <icon>images\dragonfly-dark-trans.png</icon>
         <readme>docs\README.md</readme>

IO.Astrodynamics.Net/IO.Astrodynamics/Propagator/Forces/AtmosphericDrag.cs

@@ -16,12 +16,18 @@ public AtmosphericDrag(Spacecraft spacecraft, CelestialBody celestialBody)
     {
         _spacecraft = spacecraft ?? throw new ArgumentNullException(nameof(spacecraft));
         _celestialBody = celestialBody ?? throw new ArgumentNullException(nameof(celestialBody));
+        if (!_celestialBody.HasAtmosphericModel)
+        {
+            throw new ArgumentException($"The celestial body {_celestialBody.Name} does not have an atmospheric model.");
+        }
+
         _areaMassRatio = _spacecraft.SectionalArea / _spacecraft.Mass;
     }
 
     public override Vector3 Apply(StateVector stateVector)
     {
-        var planetodetic = stateVector.RelativeTo(_celestialBody,Aberration.None).ToPlanetocentric(Aberration.None).ToPlanetodetic(_celestialBody!.Flattening, _celestialBody.EquatorialRadius);
+        var planetodetic = stateVector.RelativeTo(_celestialBody, Aberration.None).ToPlanetocentric(Aberration.None)
+            .ToPlanetodetic(_celestialBody!.Flattening, _celestialBody.EquatorialRadius);
         var density = _celestialBody.GetAirDensity(planetodetic.Altitude);
         return stateVector.Velocity * -0.5 * density * _areaMassRatio * _spacecraft.DragCoefficient * stateVector.Velocity.Magnitude();
     }

IO.Astrodynamics.Net/IO.Astrodynamics/Propagator/Forces/SolarRadiationPressure.cs

@@ -1,4 +1,6 @@
-using IO.Astrodynamics.Body;
+using System;
+using System.Collections.Generic;
+using IO.Astrodynamics.Body;
 using IO.Astrodynamics.Body.Spacecraft;
 using IO.Astrodynamics.OrbitalParameters;
 using Vector3 = IO.Astrodynamics.Math.Vector3;
@@ -10,18 +12,23 @@ public class SolarRadiationPressure : ForceBase
     private readonly CelestialBody _sun = new CelestialBody(10);
     private readonly double _areaMassRatio;
     private readonly double _term1;
+    private readonly IEnumerable<CelestialBody> _occultingBodies;
 
-    public SolarRadiationPressure(Spacecraft spacecraft)
+    public SolarRadiationPressure(Spacecraft spacecraft, IEnumerable<CelestialBody> occultingBodies)
     {
+        _occultingBodies = occultingBodies ?? throw new ArgumentNullException(nameof(occultingBodies));
         _areaMassRatio = spacecraft.SectionalArea / spacecraft.Mass;
         _term1 = Constants.SolarMeanRadiativeLuminosity / (4.0 * System.Math.PI * Constants.C);
     }
 
     public override Vector3 Apply(StateVector stateVector)
     {
-        if (_sun.IsOcculted(stateVector.Observer as CelestialItem, stateVector) == OccultationType.Full)
+        foreach (var occultingBody in _occultingBodies)
         {
-            return Vector3.Zero;
+            if (_sun.IsOcculted(occultingBody, stateVector) == OccultationType.Full)
+            {
+                return Vector3.Zero;
+            }
         }
 
         var position = stateVector.RelativeTo(_sun, Aberration.LT).ToStateVector().Position;

IO.Astrodynamics.Net/IO.Astrodynamics/Propagator/SpacecraftPropagator.cs

@@ -77,15 +77,15 @@ private List<ForceBase> InitializeForces(bool includeAtmosphericDrag, bool inclu
 
         if (includeAtmosphericDrag)
         {
-            if (_originalObserver is CelestialBody body)
+            foreach (var celestialBody in CelestialItems.OfType<CelestialBody>().Where(x=>x.HasAtmosphericModel).ToArray())
             {
-                forces.Add(new AtmosphericDrag(Spacecraft, body));
+                forces.Add(new AtmosphericDrag(Spacecraft, celestialBody));
             }
         }
 
         if (includeSolarRadiationPressure)
         {
-            forces.Add(new SolarRadiationPressure(Spacecraft));
+            forces.Add(new SolarRadiationPressure(Spacecraft,CelestialItems.OfType<CelestialBody>()));
         }
 
         return forces;