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