Commit 1c6184e2 authored by Bryan Cazabonne's avatar Bryan Cazabonne

Merge branch 'issue-752' into 'develop'

Added new method signature in IodLambert using Position measurement.

Closes #752

See merge request orekit/orekit!129
parents 654593df 4b62747f
Pipeline #892 passed with stages
in 33 minutes and 26 seconds
......@@ -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>
......
......@@ -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
......
......@@ -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");
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment