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 latitude arguments used by {@link org.orekit.orbits.CircularOrbit}.
27 * @author Romain Serra
28 * @see org.orekit.orbits.CircularOrbit
29 * @since 12.1
30 */
31 public class CircularLatitudeArgumentUtility {
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 CircularLatitudeArgumentUtility() {
41 // nothing here (utils class)
42 }
43
44 /**
45 * Computes the true latitude argument from the eccentric latitude argument.
46 *
47 * @param ex e cos(ω), first component of circular eccentricity vector
48 * @param ey e sin(ω), second component of circular eccentricity vector
49 * @param alphaE = E + ω eccentric latitude argument (rad)
50 * @return the true latitude argument.
51 */
52 public static double eccentricToTrue(final double ex, final double ey, final double alphaE) {
53 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
54 final SinCos scAlphaE = FastMath.sinCos(alphaE);
55 final double num = ex * scAlphaE.sin() - ey * scAlphaE.cos();
56 final double den = epsilon + 1 - ex * scAlphaE.cos() - ey * scAlphaE.sin();
57 return alphaE + eccentricAndTrueAtan(num, den);
58 }
59
60 /**
61 * Computes the eccentric latitude argument from the true latitude argument.
62 *
63 * @param ex e cos(ω), first component of circular eccentricity vector
64 * @param ey e sin(ω), second component of circular eccentricity vector
65 * @param alphaV = V + ω true latitude argument (rad)
66 * @return the eccentric latitude argument.
67 */
68 public static double trueToEccentric(final double ex, final double ey, final double alphaV) {
69 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
70 final SinCos scAlphaV = FastMath.sinCos(alphaV);
71 final double num = ey * scAlphaV.cos() - ex * scAlphaV.sin();
72 final double den = epsilon + 1 + ex * scAlphaV.cos() + ey * scAlphaV.sin();
73 return alphaV + 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 circular eccentricity vector
80 * @param ey e sin(ω), second component of circular 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 latitude argument from the mean latitude argument.
100 *
101 * @param ex e cos(ω), first component of circular eccentricity vector
102 * @param ey e sin(ω), second component of circular eccentricity vector
103 * @param alphaM = M + ω mean latitude argument (rad)
104 * @return the eccentric latitude argument.
105 */
106 public static double meanToEccentric(final double ex, final double ey, final double alphaM) {
107 // Generalization of Kepler equation to circular parameters
108 // with alphaE = PA + E and
109 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
110 double alphaE = alphaM;
111 double shift;
112 double alphaEMalphaM = 0.0;
113 boolean hasConverged;
114 int iter = 0;
115 do {
116 final SinCos scAlphaE = FastMath.sinCos(alphaE);
117 final double f2 = ex * scAlphaE.sin() - ey * scAlphaE.cos();
118 final double f1 = 1.0 - ex * scAlphaE.cos() - ey * scAlphaE.sin();
119 final double f0 = alphaEMalphaM - f2;
120
121 final double f12 = 2.0 * f1;
122 shift = f0 * f12 / (f1 * f12 - f0 * f2);
123
124 alphaEMalphaM -= shift;
125 alphaE = alphaM + alphaEMalphaM;
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_LATITUDE_ARGUMENT, iter);
132 }
133 return alphaE;
134
135 }
136
137 /**
138 * Computes the mean latitude argument from the eccentric latitude argument.
139 *
140 * @param ex e cos(ω), first component of circular eccentricity vector
141 * @param ey e sin(ω), second component of circular eccentricity vector
142 * @param alphaE = E + ω mean latitude argument (rad)
143 * @return the mean latitude argument.
144 */
145 public static double eccentricToMean(final double ex, final double ey, final double alphaE) {
146 final SinCos scAlphaE = FastMath.sinCos(alphaE);
147 return alphaE + (ey * scAlphaE.cos() - ex * scAlphaE.sin());
148 }
149
150 /**
151 * Computes the mean latitude argument from the eccentric latitude argument.
152 *
153 * @param ex e cos(ω), first component of circular eccentricity vector
154 * @param ey e sin(ω), second component of circular eccentricity vector
155 * @param alphaV = V + ω true latitude argument (rad)
156 * @return the mean latitude argument.
157 */
158 public static double trueToMean(final double ex, final double ey, final double alphaV) {
159 final double alphaE = trueToEccentric(ex, ey, alphaV);
160 return eccentricToMean(ex, ey, alphaE);
161 }
162
163 /**
164 * Computes the true latitude argument from the eccentric latitude argument.
165 *
166 * @param ex e cos(ω), first component of circular eccentricity vector
167 * @param ey e sin(ω), second component of circular eccentricity vector
168 * @param alphaM = M + ω mean latitude argument (rad)
169 * @return the true latitude argument.
170 */
171 public static double meanToTrue(final double ex, final double ey, final double alphaM) {
172 final double alphaE = meanToEccentric(ex, ey, alphaM);
173 return eccentricToTrue(ex, ey, alphaE);
174 }
175
176 /**
177 * Convert argument of latitude.
178 * @param oldType old position angle type
179 * @param alpha old value for argument of latitude
180 * @param ex ex
181 * @param ey ey
182 * @param newType new position angle type
183 * @return convert argument of latitude
184 * @since 12.2
185 */
186 public static double convertAlpha(final PositionAngleType oldType, final double alpha, final double ex,
187 final double ey, final PositionAngleType newType) {
188 if (oldType == newType) {
189 return alpha;
190
191 } else {
192 switch (newType) {
193
194 case ECCENTRIC:
195 if (oldType == PositionAngleType.MEAN) {
196 return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alpha);
197 } else {
198 return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alpha);
199 }
200
201 case MEAN:
202 if (oldType == PositionAngleType.TRUE) {
203 return CircularLatitudeArgumentUtility.trueToMean(ex, ey, alpha);
204 } else {
205 return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alpha);
206 }
207
208 case TRUE:
209 if (oldType == PositionAngleType.MEAN) {
210 return CircularLatitudeArgumentUtility.meanToTrue(ex, ey, alpha);
211 } else {
212 return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alpha);
213 }
214
215 default:
216 throw new OrekitInternalError(null);
217 }
218 }
219 }
220 }