1 /* Copyright 2022-2025 Romain Serra
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.orbits;
18
19 import org.hipparchus.util.FastMath;
20 import org.hipparchus.util.SinCos;
21 import org.orekit.errors.OrekitException;
22 import org.orekit.errors.OrekitInternalError;
23 import org.orekit.errors.OrekitMessages;
24
25 /**
26 * Utility methods for converting between different longitude arguments used by {@link EquinoctialOrbit}.
27 * @author Romain Serra
28 * @see EquinoctialOrbit
29 * @since 12.1
30 */
31 public class EquinoctialLongitudeArgumentUtility {
32
33 /** Tolerance for stopping criterion in iterative conversion from mean to eccentric angle. */
34 private static final double TOLERANCE_CONVERGENCE = 1.0e-11;
35
36 /** Maximum number of iterations in iterative conversion from mean to eccentric angle. */
37 private static final int MAXIMUM_ITERATION = 50;
38
39 /** Private constructor for utility class. */
40 private EquinoctialLongitudeArgumentUtility() {
41 // nothing here (utils class)
42 }
43
44 /**
45 * Computes the true longitude argument from the eccentric longitude argument.
46 *
47 * @param ex e cos(ω), first component of eccentricity vector
48 * @param ey e sin(ω), second component of eccentricity vector
49 * @param lE = E + ω + Ω eccentric longitude argument (rad)
50 * @return the true longitude argument.
51 */
52 public static double eccentricToTrue(final double ex, final double ey, final double lE) {
53 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
54 final SinCos scLE = FastMath.sinCos(lE);
55 final double num = ex * scLE.sin() - ey * scLE.cos();
56 final double den = epsilon + 1 - ex * scLE.cos() - ey * scLE.sin();
57 return lE + eccentricAndTrueAtan(num, den);
58 }
59
60 /**
61 * Computes the eccentric longitude argument from the true longitude argument.
62 *
63 * @param ex e cos(ω), first component of eccentricity vector
64 * @param ey e sin(ω), second component of eccentricity vector
65 * @param lV = V + ω + Ω true longitude argument (rad)
66 * @return the eccentric longitude argument.
67 */
68 public static double trueToEccentric(final double ex, final double ey, final double lV) {
69 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
70 final SinCos scLv = FastMath.sinCos(lV);
71 final double num = ey * scLv.cos() - ex * scLv.sin();
72 final double den = epsilon + 1 + ex * scLv.cos() + ey * scLv.sin();
73 return lV + eccentricAndTrueAtan(num, den);
74 }
75
76 /**
77 * Computes an intermediate quantity for conversions between true and eccentric.
78 *
79 * @param ex e cos(ω), first component of eccentricity vector
80 * @param ey e sin(ω), second component of eccentricity vector
81 * @return intermediate variable referred to as epsilon.
82 */
83 private static double eccentricAndTrueEpsilon(final double ex, final double ey) {
84 return FastMath.sqrt(1 - ex * ex - ey * ey);
85 }
86
87 /**
88 * Computes another intermediate quantity for conversions between true and eccentric.
89 *
90 * @param num numerator for angular conversion
91 * @param den denominator for angular conversion
92 * @return arc-tangent of ratio of inputs times two.
93 */
94 private static double eccentricAndTrueAtan(final double num, final double den) {
95 return 2. * FastMath.atan(num / den);
96 }
97
98 /**
99 * Computes the eccentric longitude argument from the mean longitude argument.
100 *
101 * @param ex e cos(ω), first component of eccentricity vector
102 * @param ey e sin(ω), second component of eccentricity vector
103 * @param lM = M + ω + Ω mean longitude argument (rad)
104 * @return the eccentric longitude argument.
105 */
106 public static double meanToEccentric(final double ex, final double ey, final double lM) {
107 // Generalization of Kepler equation to equinoctial parameters
108 // with lE = PA + RAAN + E and
109 // lM = PA + RAAN + M = lE - ex.sin(lE) + ey.cos(lE)
110 double lE = lM;
111 double shift;
112 double lEmlM = 0.0;
113 boolean hasConverged;
114 int iter = 0;
115 do {
116 final SinCos scLE = FastMath.sinCos(lE);
117 final double f2 = ex * scLE.sin() - ey * scLE.cos();
118 final double f1 = 1.0 - ex * scLE.cos() - ey * scLE.sin();
119 final double f0 = lEmlM - f2;
120
121 final double f12 = 2.0 * f1;
122 shift = f0 * f12 / (f1 * f12 - f0 * f2);
123
124 lEmlM -= shift;
125 lE = lM + lEmlM;
126
127 hasConverged = FastMath.abs(shift) <= TOLERANCE_CONVERGENCE;
128 } while (++iter < MAXIMUM_ITERATION && !hasConverged);
129
130 if (!hasConverged) {
131 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECCENTRIC_LONGITUDE_ARGUMENT, iter);
132 }
133 return lE;
134
135 }
136
137 /**
138 * Computes the mean longitude argument from the eccentric longitude argument.
139 *
140 * @param ex e cos(ω), first component of eccentricity vector
141 * @param ey e sin(ω), second component of eccentricity vector
142 * @param lE = E + ω + Ω mean longitude argument (rad)
143 * @return the mean longitude argument.
144 */
145 public static double eccentricToMean(final double ex, final double ey, final double lE) {
146 final SinCos scLE = FastMath.sinCos(lE);
147 return lE - ex * scLE.sin() + ey * scLE.cos();
148 }
149
150 /**
151 * Computes the mean longitude argument from the eccentric longitude argument.
152 *
153 * @param ex e cos(ω), first component of eccentricity vector
154 * @param ey e sin(ω), second component of eccentricity vector
155 * @param lV = V + ω + Ω true longitude argument (rad)
156 * @return the mean longitude argument.
157 */
158 public static double trueToMean(final double ex, final double ey, final double lV) {
159 final double lE = trueToEccentric(ex, ey, lV);
160 return eccentricToMean(ex, ey, lE);
161 }
162
163 /**
164 * Computes the true longitude argument from the eccentric longitude argument.
165 *
166 * @param ex e cos(ω), first component of eccentricity vector
167 * @param ey e sin(ω), second component of eccentricity vector
168 * @param lM = M + ω + Ω mean longitude argument (rad)
169 * @return the true longitude argument.
170 */
171 public static double meanToTrue(final double ex, final double ey, final double lM) {
172 final double lE = meanToEccentric(ex, ey, lM);
173 return eccentricToTrue(ex, ey, lE);
174 }
175
176 /**
177 * Convert argument of longitude.
178 * @param oldType old position angle type
179 * @param l old value for argument of longitude
180 * @param ex ex
181 * @param ey ey
182 * @param newType new position angle type
183 * @return converted argument of longitude
184 * @since 12.2
185 */
186 public static double convertL(final PositionAngleType oldType, final double l, final double ex,
187 final double ey, final PositionAngleType newType) {
188 if (oldType == newType) {
189 return l;
190
191 } else {
192 switch (newType) {
193
194 case ECCENTRIC:
195 if (oldType == PositionAngleType.MEAN) {
196 return EquinoctialLongitudeArgumentUtility.meanToEccentric(ex, ey, l);
197 } else {
198 return EquinoctialLongitudeArgumentUtility.trueToEccentric(ex, ey, l);
199 }
200
201 case MEAN:
202 if (oldType == PositionAngleType.TRUE) {
203 return EquinoctialLongitudeArgumentUtility.trueToMean(ex, ey, l);
204 } else {
205 return EquinoctialLongitudeArgumentUtility.eccentricToMean(ex, ey, l);
206 }
207
208 case TRUE:
209 if (oldType == PositionAngleType.MEAN) {
210 return EquinoctialLongitudeArgumentUtility.meanToTrue(ex, ey, l);
211 } else {
212 return EquinoctialLongitudeArgumentUtility.eccentricToTrue(ex, ey, l);
213 }
214
215 default:
216 throw new OrekitInternalError(null);
217 }
218 }
219 }
220 }