Hi all,
I am trying to solve the classical Lambert’s problem by using the IodLambert class of Orekit.
As I was getting some unexpected results by using this class, I wrote a code in which I am trying to solve two simple scenarios to test IodLambert.
In terms of position-velocity coordinates, let’s identify with S1 :(p1,v1) the initial state at date t1, and with S2 :(p2,v2) the final state at date t2.
Let also K be the Keplerian orbit passing by S1 and S2.
In both the scenarios, S1 corresponds to the state at the periapsis of the orbit K. Moreover, the motion is always consider counter-clockwise.
In the first scenario, the angle swept when going from p1 to p2 is 90 deg.
In the second scenario, instead, the angle swept when going from p1 to p2 corresponds to 270 deg.
As both S1 and S2 belong to the same Keplerian orbit K, the output Keplerian orbit KL of the “iodLambert.estimate” method (with initial and final input given respectively by (p1,t1) and (p2,t2)) should be such that
KL = K
This occurs in the first scenario, but not in the second one.
In fact, in the second scenario it seems – in first approximation - that the Lambert arc follows a clockwise (rather than a counter-clockwise) direction.
However, by having a closer look to K and KL it can be noticed a certain difference between these two orbits, not only given by their opposite motion direction.
In the “iodLambert.estimate” method (and for each test case) the input parameters “posigrade” and “nRev” have been set according to the documentation.
Just to give a try, in the second test case I also attempted to use different combinations of “posigrade” and “nRev”, but giving me in all the cases wrong results (i.e. KL ≠ K).
I am attaching here below the code I wrote to test the IodLambert class: you can set “caseToTest” parameter to 1 to test the first scenario, or to 2 to test the second scenario.
I would be really grateful if you could tell me whether there is something I am setting wrong to use the IodLambert class, or if maybe there is a bug in the implementation of this Orekit class.
Many thanks in advance! 
Nicola.
public class TestLambert {
public static void main(String args[]) throws Exception {
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(args[0])));
double muEarth = Constants.EGM96_EARTH_MU;
// According to the Orekit documentation for IodLambert.estimate:
/*
* "As an example, if t2 is less than half a period after t1, then posigrade should be true and nRev
* should be 0. If t2 is more than half a period after t1 but less than one period after t1,
* posigrade should be false and nRev should be 1."
*/
int caseToTest = 2;
double trueAnomalyDifference = 0;
boolean posigrade = false;
int nRev = 0;
if (caseToTest == 1) {
// This test case works as expected
posigrade = true;
nRev = 0;
trueAnomalyDifference = 0.5*FastMath.PI; // 90 degrees;
} else {
// Gives unexpeted results -> the errors in position are fine, but the errors in velocity are huge.
posigrade = false;
nRev = 1;
trueAnomalyDifference = 1.5*FastMath.PI; // 270 degrees;
}
// Assign position and velocity @t1 (defined as the position and velocity vectors and the periapse)
AbsoluteDate t1 = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
KeplerianOrbit kep1 = new KeplerianOrbit(24000*1e3, 0.72, 0, 0, 0, 0, PositionAngle.TRUE, FramesFactory.getEME2000(), t1, muEarth);
Vector3D p1 = kep1.getPVCoordinates().getPosition();
Vector3D v1 = kep1.getPVCoordinates().getVelocity();
// Assign t2 date (defined as the date after a swept angle of "trueAnomalyDifference" after periapse)
KeplerianOrbit kep2 = new KeplerianOrbit(24000*1e3, 0.72, 0, 0, 0, trueAnomalyDifference, PositionAngle.TRUE, FramesFactory.getEME2000(), t1, muEarth);
double n = kep2.getKeplerianMeanMotion(); //2 * FastMath.PI / kep2.getKeplerianPeriod();
double Me = kep2.getMeanAnomaly();
double delta_t = Me / n; // seconds
AbsoluteDate t2 = t1.shiftedBy(delta_t);
// Assign position and velocity @t2
PVCoordinates pv2 = kep1.getPVCoordinates(t2, FramesFactory.getEME2000());
Vector3D p2 = pv2.getPosition();
Vector3D v2 = pv2.getVelocity();
// // prove trueAnomalyDifference
// CartesianOrbit cart = new CartesianOrbit(pv2, FramesFactory.getEME2000(), t2, muEarth);
// KeplerianOrbit kep3 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(cart);
// System.out.println(FastMath.toDegrees(kep3.getTrueAnomaly()));
// Solve Lambert problem
IodLambert iodLambert = new IodLambert(muEarth);
KeplerianOrbit resultOrbit = iodLambert.estimate(FramesFactory.getEME2000(), posigrade, nRev, p1, t1, p2, t2);
// Get position and velocity coordinates @t1 and @t2 from the output Keplerian orbit of iodLambert.estimate
PVCoordinates resultPv1 = resultOrbit.getPVCoordinates(t1, FramesFactory.getEME2000());
PVCoordinates resultPv2 = resultOrbit.getPVCoordinates(t2, FramesFactory.getEME2000());
Vector3D resultP1 = resultPv1.getPosition();
Vector3D resultV1 = resultPv1.getVelocity();
Vector3D resultP2 = resultPv2.getPosition();
Vector3D resultV2 = resultPv2.getVelocity();
// Print results
System.out.println("Input position vector @t1 to iodLambert.estimate:");
System.out.println("positionVector1 = " + p1 + " m");
System.out.println("Position vector @t1 retrieved from the output Keplerian orbit of iodLambert.estimate:");
System.out.println("lambertPositionVector1 = " + resultP1 + " m\n");
System.out.println("Velocity vector associated to the input position vector @t1 to iodLambert.estimate:");
System.out.println("velocityVector1 = " + v1 + " m/s");
System.out.println("Velocity vector @t1 retrieved from the output Keplerian orbit of iodLambert.estimate:");
System.out.println("lambertVelocityVector1 = " + resultV1 + " m/s");
System.out.println("\n");
System.out.println("||positionVector1 - lambertPositionVector1|| -> should always be numerically close to 0:");
System.out.println("deltaPosition1 = " + p1.subtract(resultP1).getNorm() + " m");
System.out.println("||velocityVector1 - lambertVelocityVector1|| -> for the current test case, should be numerically close to 0:");
System.out.println("deltaV1 = " + v1.subtract(resultV1).getNorm() + " m/s");
System.out.println("\n");
System.out.println("Input position vector @t2 to iodLambert.estimate:");
System.out.println("positionVector2 = " + p2 + " m");
System.out.println("Position vector @t2 retrieved from the output Keplerian orbit of iodLambert.estimate:");
System.out.println("lambertPositionVector2 = " + resultP2 + " m\n");
System.out.println("Velocity vector associated to the input position vector @t2 to iodLambert.estimate:");
System.out.println("velocityVector2 = " + v2 + " m/s");
System.out.println("Velocity vector @t2 retrieved from the output Keplerian orbit of iodLambert.estimate:");
System.out.println("lambertVelocityVector2 = " + resultV2 + " m/s");
System.out.println("\n");
System.out.println("||positionVector2 - lambertPositionVector2|| -> should always be numerically close to 0:");
System.out.println("deltaPosition2 = " + p2.subtract(resultP2).getNorm() + " m");
System.out.println("||velocityVector2 - lambertVelocityVector2|| -> for the current test case, should be numerically close to 0:");
System.out.println("deltaV2 = " + v2.subtract(resultV2).getNorm() + " m/s");
System.out.println("\n");
System.out.println("Angular momentum (per unit of mass) of the Keplerian orbit considered in the current test case:");
System.out.println("h = " + kep1.getPVCoordinates().getMomentum());
System.out.println("Angular momentum (per unit of mass) of the output Keplerian orbit of iodLambert.estimate:");
System.out.println("h = " + resultOrbit.getPVCoordinates().getMomentum());
}
}