Merge pull request #231

* Compute time to given radius

Author: sylvain-guillet

IO.Astrodynamics.Net/IO.Astrodynamics.Tests/OrbitalParameters/KeplerianElementsTests.cs

@@ -458,7 +458,7 @@ public void ParabolicToState()
             var ke2 = sv.ToKeplerianElements();
             Assert.Equal(ke, ke2, TestHelpers.KeplerComparer);
         }
-        
+
         [Fact]
         public void EllipticToState()
         {
@@ -474,7 +474,7 @@ public void EllipticToState()
             Assert.Equal(ke.Epoch, ke2.Epoch);
             Assert.Equal(ke.Observer, ke2.Observer);
         }
-        
+
         [Fact]
         public void HyperbolicToEquinoctial()
         {
@@ -493,14 +493,14 @@ public void ParabolicToEquinoctial()
             var ke2 = eq.ToKeplerianElements();
             Assert.Equal(ke, ke2, TestHelpers.KeplerComparer);
         }
-        
+
         [Fact]
         public void EllipticToEquinoctial()
         {
             var ke = new KeplerianElements(6800000.0, 0.2, 0.2, 0.3, 0.4, 0.5, TestHelpers.EarthAtJ2000, TimeSystem.Time.J2000TDB, Frames.Frame.ICRF);
             var eq = ke.ToEquinoctial();
             var ke2 = eq.ToKeplerianElements();
-            Assert.Equal(ke,ke2, TestHelpers.KeplerComparer);
+            Assert.Equal(ke, ke2, TestHelpers.KeplerComparer);
         }
 
         [Fact]
@@ -522,5 +522,31 @@ public void ParabolicPerigeeRadius()
                 new TimeSystem.Time(DateTime.Now, TimeFrame.TDBFrame), Frames.Frame.ICRF, perigeeRadius: 6800000.0);
             Assert.Equal(6800000.0, ke.PerigeeRadius(), 6);
         }
+
+        [Fact]
+        public void TimeToApogeeRadius()
+        {
+            var epoch = TimeSystem.Time.J2000TDB;
+            KeplerianElements ke = new KeplerianElements(8000000.0, 0.5, 30.0 * IO.Astrodynamics.Constants.Deg2Rad,
+                40.0 * IO.Astrodynamics.Constants.Deg2Rad, 50.0 * IO.Astrodynamics.Constants.Deg2Rad, 0.0, TestHelpers.EarthAtJ2000,
+                epoch, Frames.Frame.ICRF);
+            var orbitalPeriod = ke.Period();
+            var expectedApogeeTime = epoch + 0.5 * orbitalPeriod;
+            var apogeeRadius = ke.ApogeeVector().Magnitude();
+            Assert.Equal(expectedApogeeTime, ke.TimeToRadius(apogeeRadius), TestHelpers.TimeComparer);
+        }
+
+        [Fact]
+        public void TimeToPerigeeRadius()
+        {
+            var epoch = TimeSystem.Time.J2000TDB;
+            KeplerianElements ke = new KeplerianElements(8000000.0, 0.5, 30.0 * IO.Astrodynamics.Constants.Deg2Rad,
+                40.0 * IO.Astrodynamics.Constants.Deg2Rad, 50.0 * IO.Astrodynamics.Constants.Deg2Rad, Astrodynamics.Constants.PI, TestHelpers.EarthAtJ2000,
+                epoch, Frames.Frame.ICRF);
+            var orbitalPeriod = ke.Period();
+            var expectedApogeeTime = epoch + 0.5 * orbitalPeriod;
+            var perigeeRadius = ke.PerigeeVector().Magnitude();
+            Assert.Equal(expectedApogeeTime, ke.TimeToRadius(perigeeRadius), TestHelpers.TimeComparer);
+        }
     }
 }

IO.Astrodynamics.Net/IO.Astrodynamics/OrbitalParameters/KeplerianElements.cs

@@ -165,6 +165,8 @@ public override double MeanAnomaly()
             return M;
         }
 
+        
+        
         #region Operator
 
         public bool Equals(KeplerianElements other)

IO.Astrodynamics.Net/IO.Astrodynamics/OrbitalParameters/OrbitalParameters.cs

@@ -935,7 +935,41 @@ public double SemiLatusRectum()
         return hNorm * hNorm / Observer.GM;
     }
 
-    
+    /// <summary>
+    /// Computes the time to reach a given radius.
+    /// </summary>
+    /// <param name="radius"></param>
+    /// <returns></returns>
+    public Time TimeToRadius(double radius)
+    {
+        // Extract orbital elements
+        double e = Eccentricity(); // Eccentricity
+
+        // Semi-latus rectum
+        double p = SemiLatusRectum();
+
+        // Compute the true anomaly at the radius
+        double cosNuRad = System.Math.Clamp((p / radius - 1) / e, -1.0, 1.0);
+        double nuRad = System.Math.Acos(cosNuRad);
+
+        return EpochAtMeanAnomaly(TrueAnomalyToMeanAnomaly(nuRad, e, EccentricAnomaly(nuRad)));
+    }
+
+    /// <summary>
+    /// Computes the time to reach the sphere of influence of the center of motion.
+    /// </summary>
+    /// <returns></returns>
+    /// <exception cref="InvalidOperationException">The observer must be a celestial body to compute the time to the sphere of influence.</exception>
+    public Time TimeToSphereOfInfluence()
+    {
+        var centerOfMotion = Observer as CelestialBody;
+        if (centerOfMotion is null)
+        {
+            throw new InvalidOperationException("The observer must be a celestial body to compute the time to the sphere of influence.");
+        }
+
+        return TimeToRadius(centerOfMotion.SphereOfInfluence);
+    }
 
     #region Operators