1 /* Copyright 2022-2024 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 latitude arguments used by {@link FieldCircularOrbit}.
28 * @author Romain Serra
29 * @see FieldCircularOrbit
30 * @since 12.1
31 */
32 public class FieldCircularLatitudeArgumentUtility {
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 FieldCircularLatitudeArgumentUtility() {
42 // nothing here (utils class)
43 }
44
45 /**
46 * Computes the true latitude argument from the eccentric latitude argument.
47 *
48 * @param <T> Type of the field elements
49 * @param ex e cos(ω), first component of circular eccentricity vector
50 * @param ey e sin(ω), second component of circular eccentricity vector
51 * @param alphaE = E + ω eccentric latitude argument (rad)
52 * @return the true latitude argument.
53 */
54 public static <T extends CalculusFieldElement<T>> T eccentricToTrue(final T ex, final T ey, final T alphaE) {
55 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
56 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
57 final T cosAlphaE = scAlphaE.cos();
58 final T sinAlphaE = scAlphaE.sin();
59 final T num = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
60 final T den = epsilon.add(1).subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
61 return alphaE.add(eccentricAndTrueAtan(num, den));
62 }
63
64 /**
65 * Computes the eccentric latitude argument from the true latitude argument.
66 *
67 * @param <T> Type of the field elements
68 * @param ex e cos(ω), first component of circular eccentricity vector
69 * @param ey e sin(ω), second component of circular eccentricity vector
70 * @param alphaV = v + ω true latitude argument (rad)
71 * @return the eccentric latitude argument.
72 */
73 public static <T extends CalculusFieldElement<T>> T trueToEccentric(final T ex, final T ey, final T alphaV) {
74 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
75 final FieldSinCos<T> scAlphaV = FastMath.sinCos(alphaV);
76 final T cosAlphaV = scAlphaV.cos();
77 final T sinAlphaV = scAlphaV.sin();
78 final T num = ey.multiply(cosAlphaV).subtract(ex.multiply(sinAlphaV));
79 final T den = epsilon.add(1).add(ex.multiply(cosAlphaV).add(ey.multiply(sinAlphaV)));
80 return alphaV.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 circular eccentricity vector
88 * @param ey e sin(ω), second component of circular 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 latitude argument from the mean latitude argument.
109 *
110 * @param <T> Type of the field elements
111 * @param ex e cos(ω), first component of circular eccentricity vector
112 * @param ey e sin(ω), second component of circular eccentricity vector
113 * @param alphaM = M + ω mean latitude argument (rad)
114 * @return the eccentric latitude argument.
115 */
116 public static <T extends CalculusFieldElement<T>> T meanToEccentric(final T ex, final T ey, final T alphaM) {
117 // Generalization of Kepler equation to circular parameters
118 // with alphaE = PA + E and
119 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
120
121 T alphaE = alphaM;
122 T shift;
123 T alphaEMalphaM = alphaM.getField().getZero();
124 boolean hasConverged;
125 int iter = 0;
126 do {
127 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
128 final T f2 = ex.multiply(scAlphaE.sin()).subtract(ey.multiply(scAlphaE.cos()));
129 final T f1 = ex.negate().multiply(scAlphaE.cos()).subtract(ey.multiply(scAlphaE.sin())).add(1);
130 final T f0 = alphaEMalphaM.subtract(f2);
131
132 final T f12 = f1.multiply(2);
133 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
134
135 alphaEMalphaM = alphaEMalphaM.subtract(shift);
136 alphaE = alphaM.add(alphaEMalphaM);
137
138 hasConverged = FastMath.abs(shift.getReal()) <= TOLERANCE_CONVERGENCE;
139 } while (++iter < MAXIMUM_ITERATION && !hasConverged);
140
141 if (!hasConverged) {
142 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECCENTRIC_LATITUDE_ARGUMENT, iter);
143 }
144 return alphaE;
145
146 }
147
148 /**
149 * Computes the mean latitude argument from the eccentric latitude argument.
150 *
151 * @param <T> Type of the field elements
152 * @param ex e cos(ω), first component of circular eccentricity vector
153 * @param ey e sin(ω), second component of circular eccentricity vector
154 * @param alphaE = E + ω eccentric latitude argument (rad)
155 * @return the mean latitude argument.
156 */
157 public static <T extends CalculusFieldElement<T>> T eccentricToMean(final T ex, final T ey, final T alphaE) {
158 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
159 return alphaE.subtract(ex.multiply(scAlphaE.sin()).subtract(ey.multiply(scAlphaE.cos())));
160 }
161
162 /**
163 * Computes the mean latitude argument from the eccentric latitude argument.
164 *
165 * @param <T> Type of the field elements
166 * @param ex e cos(ω), first component of circular eccentricity vector
167 * @param ey e sin(ω), second component of circular eccentricity vector
168 * @param alphaV = V + ω true latitude argument (rad)
169 * @return the mean latitude argument.
170 */
171 public static <T extends CalculusFieldElement<T>> T trueToMean(final T ex, final T ey, final T alphaV) {
172 final T alphaE = trueToEccentric(ex, ey, alphaV);
173 return eccentricToMean(ex, ey, alphaE);
174 }
175
176 /**
177 * Computes the true latitude argument from the eccentric latitude argument.
178 *
179 * @param <T> Type of the field elements
180 * @param ex e cos(ω), first component of circular eccentricity vector
181 * @param ey e sin(ω), second component of circular eccentricity vector
182 * @param alphaM = M + ω mean latitude argument (rad)
183 * @return the true latitude argument.
184 */
185 public static <T extends CalculusFieldElement<T>> T meanToTrue(final T ex, final T ey, final T alphaM) {
186 final T alphaE = meanToEccentric(ex, ey, alphaM);
187 return eccentricToTrue(ex, ey, alphaE);
188 }
189
190 /**
191 * Convert argument of latitude.
192 * @param oldType old position angle type
193 * @param alpha old value for argument of latitude
194 * @param ex ex
195 * @param ey ey
196 * @param newType new position angle type
197 * @param <T> field type
198 * @return convert argument of latitude
199 * @since 12.2
200 */
201 public static <T extends CalculusFieldElement<T>> T convertAlpha(final PositionAngleType oldType, final T alpha,
202 final T ex, final T ey,
203 final PositionAngleType newType) {
204 if (oldType == newType) {
205 return alpha;
206
207 } else {
208 switch (newType) {
209
210 case ECCENTRIC:
211 if (oldType == PositionAngleType.MEAN) {
212 return FieldCircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alpha);
213 } else {
214 return FieldCircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alpha);
215 }
216
217 case MEAN:
218 if (oldType == PositionAngleType.TRUE) {
219 return FieldCircularLatitudeArgumentUtility.trueToMean(ex, ey, alpha);
220 } else {
221 return FieldCircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alpha);
222 }
223
224 case TRUE:
225 if (oldType == PositionAngleType.MEAN) {
226 return FieldCircularLatitudeArgumentUtility.meanToTrue(ex, ey, alpha);
227 } else {
228 return FieldCircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alpha);
229 }
230
231 default:
232 throw new OrekitInternalError(null);
233 }
234 }
235 }
236 }