1 /* Copyright 2002-2025 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.estimation.iod;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.geometry.euclidean.twod.Vector2D;
21 import org.orekit.control.heuristics.lambert.LambertBoundaryConditions;
22 import org.orekit.control.heuristics.lambert.LambertBoundaryVelocities;
23 import org.orekit.control.heuristics.lambert.LambertSolver;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.errors.OrekitMessages;
26 import org.orekit.estimation.measurements.PV;
27 import org.orekit.estimation.measurements.Position;
28 import org.orekit.frames.Frame;
29 import org.orekit.orbits.CartesianOrbit;
30 import org.orekit.orbits.Orbit;
31 import org.orekit.time.AbsoluteDate;
32 import org.orekit.utils.TimeStampedPVCoordinates;
33
34 /**
35 * Lambert position-based Initial Orbit Determination (IOD) algorithm, assuming Keplerian motion.
36 * <p>
37 * An orbit is determined from two position vectors.
38 * @author Joris Olympio
39 * @see LambertSolver
40 * @since 8.0
41 */
42 public class IodLambert {
43
44 /** gravitational constant. */
45 private final double mu;
46
47 /** Creator.
48 *
49 * @param mu gravitational constant
50 */
51 public IodLambert(final double mu) {
52 this.mu = mu;
53 }
54
55 /** Estimate an initial orbit from two position measurements.
56 * <p>
57 * The logic for setting {@code posigrade} and {@code nRev} is that the
58 * sweep angle Δυ travelled by the object between {@code t1} and {@code t2} is
59 * 2π {@code nRev +1} - α if {@code posigrade} is false and 2π {@code nRev} + α
60 * if {@code posigrade} is true, where α is the separation angle between
61 * {@code p1} and {@code p2}, which is always computed between 0 and π
62 * (because in 3D without a normal reference, vector angles cannot go past π).
63 * </p>
64 * <p>
65 * This implies that {@code posigrade} should be set to true if {@code p2} is
66 * located in the half orbit starting at {@code p1} and it should be set to
67 * false if {@code p2} is located in the half orbit ending at {@code p1},
68 * regardless of the number of periods between {@code t1} and {@code t2},
69 * and {@code nRev} should be set accordingly.
70 * </p>
71 * <p>
72 * As an example, if {@code t2} is less than half a period after {@code t1},
73 * then {@code posigrade} should be {@code true} and {@code nRev} should be 0.
74 * If {@code t2} is more than half a period after {@code t1} but less than
75 * one period after {@code t1}, {@code posigrade} should be {@code false} and
76 * {@code nRev} should be 0.
77 * </p>
78 * @param frame measurements frame
79 * @param posigrade flag indicating the direction of motion
80 * @param nRev number of revolutions
81 * @param p1 first position measurement
82 * @param p2 second position measurement
83 * @return an initial Keplerian orbit estimation at the first observation date t1
84 * @since 11.0
85 */
86 public Orbit estimate(final Frame frame, final boolean posigrade,
87 final int nRev, final Position p1, final Position p2) {
88 return estimate(frame, posigrade, nRev,
89 p1.getPosition(), p1.getDate(), p2.getPosition(), p2.getDate());
90 }
91
92 /** Estimate an initial orbit from two PV measurements.
93 * <p>
94 * The logic for setting {@code posigrade} and {@code nRev} is that the
95 * sweep angle Δυ travelled by the object between {@code t1} and {@code t2} is
96 * 2π {@code nRev +1} - α if {@code posigrade} is false and 2π {@code nRev} + α
97 * if {@code posigrade} is true, where α is the separation angle between
98 * {@code p1} and {@code p2}, which is always computed between 0 and π
99 * (because in 3D without a normal reference, vector angles cannot go past π).
100 * </p>
101 * <p>
102 * This implies that {@code posigrade} should be set to true if {@code p2} is
103 * located in the half orbit starting at {@code p1} and it should be set to
104 * false if {@code p2} is located in the half orbit ending at {@code p1},
105 * regardless of the number of periods between {@code t1} and {@code t2},
106 * and {@code nRev} should be set accordingly.
107 * </p>
108 * <p>
109 * As an example, if {@code t2} is less than half a period after {@code t1},
110 * then {@code posigrade} should be {@code true} and {@code nRev} should be 0.
111 * If {@code t2} is more than half a period after {@code t1} but less than
112 * one period after {@code t1}, {@code posigrade} should be {@code false} and
113 * {@code nRev} should be 0.
114 * </p>
115 * @param frame measurements frame
116 * @param posigrade flag indicating the direction of motion
117 * @param nRev number of revolutions
118 * @param pv1 first PV measurement
119 * @param pv2 second PV measurement
120 * @return an initial Keplerian orbit estimation at the first observation date t1
121 * @since 12.0
122 */
123 public Orbit estimate(final Frame frame, final boolean posigrade,
124 final int nRev, final PV pv1, final PV pv2) {
125 return estimate(frame, posigrade, nRev,
126 pv1.getPosition(), pv1.getDate(), pv2.getPosition(), pv2.getDate());
127 }
128
129 /** Estimate a Keplerian orbit given two position vectors and a duration.
130 * <p>
131 * The logic for setting {@code posigrade} and {@code nRev} is that the
132 * sweep angle Δυ travelled by the object between {@code t1} and {@code t2} is
133 * 2π {@code nRev +1} - α if {@code posigrade} is false and 2π {@code nRev} + α
134 * if {@code posigrade} is true, where α is the separation angle between
135 * {@code p1} and {@code p2}, which is always computed between 0 and π
136 * (because in 3D without a normal reference, vector angles cannot go past π).
137 * </p>
138 * <p>
139 * This implies that {@code posigrade} should be set to true if {@code p2} is
140 * located in the half orbit starting at {@code p1} and it should be set to
141 * false if {@code p2} is located in the half orbit ending at {@code p1},
142 * regardless of the number of periods between {@code t1} and {@code t2},
143 * and {@code nRev} should be set accordingly.
144 * </p>
145 * <p>
146 * As an example, if {@code t2} is less than half a period after {@code t1},
147 * then {@code posigrade} should be {@code true} and {@code nRev} should be 0.
148 * If {@code t2} is more than half a period after {@code t1} but less than
149 * one period after {@code t1}, {@code posigrade} should be {@code false} and
150 * {@code nRev} should be 0.
151 * </p>
152 * @param frame frame
153 * @param posigrade flag indicating the direction of motion
154 * @param nRev number of revolutions
155 * @param p1 position vector 1
156 * @param t1 date of observation 1
157 * @param p2 position vector 2
158 * @param t2 date of observation 2
159 * @return an initial Keplerian orbit estimate at the first observation date t1
160 */
161 public Orbit estimate(final Frame frame, final boolean posigrade,
162 final int nRev,
163 final Vector3D p1, final AbsoluteDate t1,
164 final Vector3D p2, final AbsoluteDate t2) {
165 // Exception if t2 < t1
166 final double tau = t2.durationFrom(t1); // in seconds
167 if (tau < 0.) {
168 throw new OrekitException(OrekitMessages.NON_CHRONOLOGICAL_DATES_FOR_OBSERVATIONS, t1, t2, -tau);
169 }
170 final LambertSolver solver = new LambertSolver(mu);
171 final LambertBoundaryConditions boundaryConditions = new LambertBoundaryConditions(t1, p1, t2, p2,
172 frame);
173 final LambertBoundaryVelocities velocities = solver.solve(posigrade, nRev, boundaryConditions);
174 if (velocities == null) {
175 return null;
176 } else {
177 return new CartesianOrbit(new TimeStampedPVCoordinates(t1, p1, velocities.getInitialVelocity()),
178 frame, mu);
179 }
180 }
181
182 /** Lambert's solver for the historical, planar problem.
183 * Assume mu=1.
184 *
185 * @param r1 radius 1
186 * @param r2 radius 2
187 * @param dth sweep angle
188 * @param tau time of flight
189 * @param mRev number of revs
190 * @param V1 velocity at departure in (T, N) basis
191 * @return exit flag
192 * @deprecated as of 13.1, use {@link LambertSolver}
193 */
194 boolean solveLambertPb(final double r1, final double r2, final double dth, final double tau,
195 final int mRev, final double[] V1) {
196 final Vector2D solution = LambertSolver.solveNormalized2D(r1, r2, dth, tau, mRev);
197 if (solution == Vector2D.NaN) {
198 return false;
199 } else {
200 V1[0] = solution.getX();
201 V1[1] = solution.getY();
202 return true;
203 }
204 }
205
206 }