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.util.FastMath;
20 import org.hipparchus.util.SinCos;
21 import org.orekit.errors.OrekitException;
22 import org.orekit.errors.OrekitMessages;
23
24 /**
25 * Utility methods for converting between different latitude arguments used by {@link org.orekit.orbits.CircularOrbit}.
26 * @author Romain Serra
27 * @see org.orekit.orbits.CircularOrbit
28 * @since 12.1
29 */
30 public class CircularLatitudeArgumentUtility {
31
32 /** Tolerance for stopping criterion in iterative conversion from mean to eccentric angle. */
33 private static final double TOLERANCE_CONVERGENCE = 1.0e-12;
34
35 /** Maximum number of iterations in iterative conversion from mean to eccentric angle. */
36 private static final int MAXIMUM_ITERATION = 50;
37
38 /** Private constructor for utility class. */
39 private CircularLatitudeArgumentUtility() {
40 // nothing here (utils class)
41 }
42
43 /**
44 * Computes the true latitude argument from the eccentric latitude argument.
45 *
46 * @param ex e cos(ω), first component of circular eccentricity vector
47 * @param ey e sin(ω), second component of circular eccentricity vector
48 * @param alphaE = E + ω eccentric latitude argument (rad)
49 * @return the true latitude argument.
50 */
51 public static double eccentricToTrue(final double ex, final double ey, final double alphaE) {
52 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
53 final SinCos scAlphaE = FastMath.sinCos(alphaE);
54 final double num = ex * scAlphaE.sin() - ey * scAlphaE.cos();
55 final double den = epsilon + 1 - ex * scAlphaE.cos() - ey * scAlphaE.sin();
56 return alphaE + eccentricAndTrueAtan(num, den);
57 }
58
59 /**
60 * Computes the eccentric latitude argument from the true latitude argument.
61 *
62 * @param ex e cos(ω), first component of circular eccentricity vector
63 * @param ey e sin(ω), second component of circular eccentricity vector
64 * @param alphaV = V + ω true latitude argument (rad)
65 * @return the eccentric latitude argument.
66 */
67 public static double trueToEccentric(final double ex, final double ey, final double alphaV) {
68 final double epsilon = eccentricAndTrueEpsilon(ex, ey);
69 final SinCos scAlphaV = FastMath.sinCos(alphaV);
70 final double num = ey * scAlphaV.cos() - ex * scAlphaV.sin();
71 final double den = epsilon + 1 + ex * scAlphaV.cos() + ey * scAlphaV.sin();
72 return alphaV + eccentricAndTrueAtan(num, den);
73 }
74
75 /**
76 * Computes an intermediate quantity for conversions between true and eccentric.
77 *
78 * @param ex e cos(ω), first component of circular eccentricity vector
79 * @param ey e sin(ω), second component of circular eccentricity vector
80 * @return intermediate variable referred to as epsilon.
81 */
82 private static double eccentricAndTrueEpsilon(final double ex, final double ey) {
83 return FastMath.sqrt(1 - ex * ex - ey * ey);
84 }
85
86 /**
87 * Computes another intermediate quantity for conversions between true and eccentric.
88 *
89 * @param num numerator for angular conversion
90 * @param den denominator for angular conversion
91 * @return arc-tangent of ratio of inputs times two.
92 */
93 private static double eccentricAndTrueAtan(final double num, final double den) {
94 return 2. * FastMath.atan(num / den);
95 }
96
97 /**
98 * Computes the eccentric latitude argument from the mean latitude argument.
99 *
100 * @param ex e cos(ω), first component of circular eccentricity vector
101 * @param ey e sin(ω), second component of circular eccentricity vector
102 * @param alphaM = M + ω mean latitude argument (rad)
103 * @return the eccentric latitude argument.
104 */
105 public static double meanToEccentric(final double ex, final double ey, final double alphaM) {
106 // Generalization of Kepler equation to circular parameters
107 // with alphaE = PA + E and
108 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
109 double alphaE = alphaM;
110 double shift;
111 double alphaEMalphaM = 0.0;
112 boolean hasConverged;
113 int iter = 0;
114 do {
115 final SinCos scAlphaE = FastMath.sinCos(alphaE);
116 final double f2 = ex * scAlphaE.sin() - ey * scAlphaE.cos();
117 final double f1 = 1.0 - ex * scAlphaE.cos() - ey * scAlphaE.sin();
118 final double f0 = alphaEMalphaM - f2;
119
120 final double f12 = 2.0 * f1;
121 shift = f0 * f12 / (f1 * f12 - f0 * f2);
122
123 alphaEMalphaM -= shift;
124 alphaE = alphaM + alphaEMalphaM;
125
126 hasConverged = FastMath.abs(shift) <= TOLERANCE_CONVERGENCE;
127 } while (++iter < MAXIMUM_ITERATION && !hasConverged);
128
129 if (!hasConverged) {
130 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECCENTRIC_LATITUDE_ARGUMENT, iter);
131 }
132 return alphaE;
133
134 }
135
136 /**
137 * Computes the mean latitude argument from the eccentric latitude argument.
138 *
139 * @param ex e cos(ω), first component of circular eccentricity vector
140 * @param ey e sin(ω), second component of circular eccentricity vector
141 * @param alphaE = E + ω mean latitude argument (rad)
142 * @return the mean latitude argument.
143 */
144 public static double eccentricToMean(final double ex, final double ey, final double alphaE) {
145 final SinCos scAlphaE = FastMath.sinCos(alphaE);
146 return alphaE + (ey * scAlphaE.cos() - ex * scAlphaE.sin());
147 }
148
149 /**
150 * Computes the mean latitude argument from the eccentric latitude argument.
151 *
152 * @param ex e cos(ω), first component of circular eccentricity vector
153 * @param ey e sin(ω), second component of circular eccentricity vector
154 * @param alphaV = V + ω true latitude argument (rad)
155 * @return the mean latitude argument.
156 */
157 public static double trueToMean(final double ex, final double ey, final double alphaV) {
158 final double alphaE = trueToEccentric(ex, ey, alphaV);
159 return eccentricToMean(ex, ey, alphaE);
160 }
161
162 /**
163 * Computes the true latitude argument from the eccentric latitude argument.
164 *
165 * @param ex e cos(ω), first component of circular eccentricity vector
166 * @param ey e sin(ω), second component of circular eccentricity vector
167 * @param alphaM = M + ω mean latitude argument (rad)
168 * @return the true latitude argument.
169 */
170 public static double meanToTrue(final double ex, final double ey, final double alphaM) {
171 final double alphaE = meanToEccentric(ex, ey, alphaM);
172 return eccentricToTrue(ex, ey, alphaE);
173 }
174
175 }