Added new method signature in IodLambert using Position measurement.

4b62747f · Bryan Cazabonne · 654593df · 4b62747f · 4b62747f · 4b62747f
Commit 4b62747f authored Feb 11, 2021 by Bryan Cazabonne
3 changed files
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
  </properties>
  <body>
    <release version="11.0" date="TBD" description="TBD">
+      <action dev="bryan" type="add" issue="752">
+        Added new method signature in IodLambert using Position measurement.
+      </action>
      <action dev="bryan" type="add" issue="751">
        Added new method signature in IodGibbs using Position measurement.
      </action>

--- a/src/main/java/org/orekit/estimation/iod/IodLambert.java
+++ b/src/main/java/org/orekit/estimation/iod/IodLambert.java
@@ -20,6 +20,7 @@ import org.hipparchus.geometry.euclidean.threed.Vector3D;
 import org.hipparchus.util.FastMath;
 import org.orekit.errors.OrekitException;
 import org.orekit.errors.OrekitMessages;
+import org.orekit.estimation.measurements.Position;
 import org.orekit.frames.Frame;
 import org.orekit.orbits.KeplerianOrbit;
 import org.orekit.time.AbsoluteDate;
@@ -49,6 +50,43 @@ public class IodLambert {
        this.mu = mu;
    }

+    /** Estimate an initial orbit from two position measurements.
+     * <p>
+     * The logic for setting {@code posigrade} and {@code nRev} is that the
+     * sweep angle Δυ travelled by the object between {@code t1} and {@code t2} is
+     * 2π {@code nRev +1} - α if {@code posigrade} is false and 2π {@code nRev} + α
+     * if {@code posigrade} is true, where α is the separation angle between
+     * {@code p1} and {@code p2}, which is always computed between 0 and π
+     * (because in 3D without a normal reference, vector angles cannot go past π).
+     * </p>
+     * <p>
+     * This implies that {@code posigrade} should be set to true if {@code p2} is
+     * located in the half orbit starting at {@code p1} and it should be set to
+     * false if {@code p2} is located in the half orbit ending at {@code p1},
+     * regardless of the number of periods between {@code t1} and {@code t2},
+     * and {@code nRev} should be set accordingly.
+     * </p>
+     * <p>
+     * As an example, if {@code t2} is less than half a period after {@code t1},
+     * then {@code posigrade} should be {@code true} and {@code nRev} should be 0.
+     * If {@code t2} is more than half a period after {@code t1} but less than
+     * one period after {@code t1}, {@code posigrade} should be {@code false} and
+     * {@code nRev} should be 0.
+     * </p>
+     * @param frame     measurements frame
+     * @param posigrade flag indicating the direction of motion
+     * @param nRev      number of revolutions
+     * @param p1        first position measurement
+     * @param p2        second position measurement
+     * @return an initial orbit estimation
+     * @since 11.0
+     */
+    public KeplerianOrbit estimate(final Frame frame, final boolean posigrade,
+                                   final int nRev, final Position p1,  final Position p2) {
+        return estimate(frame, posigrade, nRev,
+                        p1.getPosition(), p1.getDate(), p2.getPosition(), p2.getDate());
+    }
+
    /** Estimate a Keplerian orbit given two position vectors and a duration.
     * <p>
     * The logic for setting {@code posigrade} and {@code nRev} is that the

--- a/src/test/java/org/orekit/estimation/iod/IodLambertTest.java
+++ b/src/test/java/org/orekit/estimation/iod/IodLambertTest.java
@@ -29,8 +29,10 @@ import org.orekit.errors.OrekitException;
 import org.orekit.errors.OrekitMessages;
 import org.orekit.estimation.Context;
 import org.orekit.estimation.EstimationTestUtils;
+import org.orekit.estimation.measurements.ObservableSatellite;
 import org.orekit.estimation.measurements.ObservedMeasurement;
 import org.orekit.estimation.measurements.PVMeasurementCreator;
+import org.orekit.estimation.measurements.Position;
 import org.orekit.frames.Frame;
 import org.orekit.frames.FramesFactory;
 import org.orekit.orbits.KeplerianOrbit;
@@ -38,10 +40,12 @@ import org.orekit.orbits.OrbitType;
 import org.orekit.orbits.PositionAngle;
 import org.orekit.propagation.Propagator;
 import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.analytical.KeplerianPropagator;
 import org.orekit.propagation.analytical.tle.TLE;
 import org.orekit.propagation.analytical.tle.TLEPropagator;
 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
 import org.orekit.time.AbsoluteDate;
+import org.orekit.time.TimeScale;
 import org.orekit.time.TimeScalesFactory;
 import org.orekit.utils.Constants;
 import org.orekit.utils.PVCoordinates;
@@ -373,6 +377,48 @@ public class IodLambertTest {
        }
    }

+    @Test
+    public void testIssue752() {
+
+        // Test taken from “Superior Lambert Algorithm” by Gim Der
+        
+        // Initial frame, time scale
+        final Frame inertialFrame = FramesFactory.getEME2000();
+        final TimeScale utc       = TimeScalesFactory.getUTC();
+
+        // Initialisation
+        final IodLambert lambert = new IodLambert(Constants.EGM96_EARTH_MU);
+
+        // Observable satellite to initialize measurements
+        final ObservableSatellite satellite = new ObservableSatellite(0);
+
+        // Observations vector (EME2000)
+        final Vector3D posR1 = new Vector3D(22592145.603, -1599915.239, -19783950.506);
+        final Vector3D posR2 = new Vector3D(1922067.697, 4054157.051, -8925727.465);
+
+        // Epoch corresponding to the observation vector
+        AbsoluteDate dateRef = new AbsoluteDate(2018, 11, 1, 0, 0, 0.0, utc);
+        AbsoluteDate date2 = dateRef.shiftedBy(36000.0);
+
+        // Reference result
+        final Vector3D velR1 = new Vector3D(2000.652697, 387.688615, -2666.947760);
+        final Vector3D velR2 = new Vector3D(-3792.46619, -1777.07641, 6856.81495);
+
+        // Lambert IOD
+        final KeplerianOrbit orbit = lambert.estimate(inertialFrame, true, 0,
+                                                      new Position(dateRef, posR1, 1.0, 1.0, satellite),
+                                                      new Position(date2,   posR2, 1.0, 1.0, satellite));
+
+        // Test for the norm of the first velocity
+        Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR1.getNorm(),  1e-3);
+
+        // Test the norm of the second velocity
+        final KeplerianPropagator kepler = new KeplerianPropagator(orbit);
+        Assert.assertEquals(0.0, kepler.propagate(date2).getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(),  1e-3);
+
+
+    }
+
    @Before
    public void setUp() {
        Utils.setDataRoot("regular-data");