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.OrekitMessages;
24
25 /**
26 * Utility methods for converting between different latitude arguments used by {@link FieldCircularOrbit}.
27 * @author Romain Serra
28 * @see FieldCircularOrbit
29 * @since 12.1
30 */
31 public class FieldCircularLatitudeArgumentUtility {
32
33 /** Tolerance for stopping criterion in iterative conversion from mean to eccentric angle. */
34 private static final double TOLERANCE_CONVERGENCE = 1.0e-12;
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 FieldCircularLatitudeArgumentUtility() {
41 // nothing here (utils class)
42 }
43
44 /**
45 * Computes the true latitude argument from the eccentric latitude argument.
46 *
47 * @param <T> Type of the field elements
48 * @param ex e cos(ω), first component of circular eccentricity vector
49 * @param ey e sin(ω), second component of circular eccentricity vector
50 * @param alphaE = E + ω eccentric latitude argument (rad)
51 * @return the true latitude argument.
52 */
53 public static <T extends CalculusFieldElement<T>> T eccentricToTrue(final T ex, final T ey, final T alphaE) {
54 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
55 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
56 final T cosAlphaE = scAlphaE.cos();
57 final T sinAlphaE = scAlphaE.sin();
58 final T num = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
59 final T den = epsilon.add(1).subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
60 return alphaE.add(eccentricAndTrueAtan(num, den));
61 }
62
63 /**
64 * Computes the eccentric latitude argument from the true latitude argument.
65 *
66 * @param <T> Type of the field elements
67 * @param ex e cos(ω), first component of circular eccentricity vector
68 * @param ey e sin(ω), second component of circular eccentricity vector
69 * @param alphaV = v + ω true latitude argument (rad)
70 * @return the eccentric latitude argument.
71 */
72 public static <T extends CalculusFieldElement<T>> T trueToEccentric(final T ex, final T ey, final T alphaV) {
73 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
74 final FieldSinCos<T> scAlphaV = FastMath.sinCos(alphaV);
75 final T cosAlphaV = scAlphaV.cos();
76 final T sinAlphaV = scAlphaV.sin();
77 final T num = ey.multiply(cosAlphaV).subtract(ex.multiply(sinAlphaV));
78 final T den = epsilon.add(1).add(ex.multiply(cosAlphaV).add(ey.multiply(sinAlphaV)));
79 return alphaV.add(eccentricAndTrueAtan(num, den));
80 }
81
82 /**
83 * Computes an intermediate quantity for conversions between true and eccentric.
84 *
85 * @param <T> Type of the field elements
86 * @param ex e cos(ω), first component of circular eccentricity vector
87 * @param ey e sin(ω), second component of circular eccentricity vector
88 * @return intermediate variable referred to as epsilon.
89 */
90 private static <T extends CalculusFieldElement<T>> T eccentricAndTrueEpsilon(final T ex, final T ey) {
91 return (ex.square().negate().subtract(ey.square()).add(1.)).sqrt();
92 }
93
94 /**
95 * Computes another intermediate quantity for conversions between true and eccentric.
96 *
97 * @param <T> Type of the field elements
98 * @param num numerator for angular conversion
99 * @param den denominator for angular conversion
100 * @return arc-tangent of ratio of inputs times two.
101 */
102 private static <T extends CalculusFieldElement<T>> T eccentricAndTrueAtan(final T num, final T den) {
103 return (num.divide(den)).atan().multiply(2);
104 }
105
106 /**
107 * Computes the eccentric latitude argument from the mean latitude argument.
108 *
109 * @param <T> Type of the field elements
110 * @param ex e cos(ω), first component of circular eccentricity vector
111 * @param ey e sin(ω), second component of circular eccentricity vector
112 * @param alphaM = M + ω mean latitude argument (rad)
113 * @return the eccentric latitude argument.
114 */
115 public static <T extends CalculusFieldElement<T>> T meanToEccentric(final T ex, final T ey, final T alphaM) {
116 // Generalization of Kepler equation to circular parameters
117 // with alphaE = PA + E and
118 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
119
120 T alphaE = alphaM;
121 T shift;
122 T alphaEMalphaM = alphaM.getField().getZero();
123 boolean hasConverged;
124 int iter = 0;
125 do {
126 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
127 final T f2 = ex.multiply(scAlphaE.sin()).subtract(ey.multiply(scAlphaE.cos()));
128 final T f1 = ex.negate().multiply(scAlphaE.cos()).subtract(ey.multiply(scAlphaE.sin())).add(1);
129 final T f0 = alphaEMalphaM.subtract(f2);
130
131 final T f12 = f1.multiply(2);
132 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
133
134 alphaEMalphaM = alphaEMalphaM.subtract(shift);
135 alphaE = alphaM.add(alphaEMalphaM);
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_LATITUDE_ARGUMENT, iter);
142 }
143 return alphaE;
144
145 }
146
147 /**
148 * Computes the mean latitude argument from the eccentric latitude argument.
149 *
150 * @param <T> Type of the field elements
151 * @param ex e cos(ω), first component of circular eccentricity vector
152 * @param ey e sin(ω), second component of circular eccentricity vector
153 * @param alphaE = E + ω eccentric latitude argument (rad)
154 * @return the mean latitude argument.
155 */
156 public static <T extends CalculusFieldElement<T>> T eccentricToMean(final T ex, final T ey, final T alphaE) {
157 final FieldSinCos<T> scAlphaE = FastMath.sinCos(alphaE);
158 return alphaE.subtract(ex.multiply(scAlphaE.sin()).subtract(ey.multiply(scAlphaE.cos())));
159 }
160
161 /**
162 * Computes the mean latitude argument from the eccentric latitude argument.
163 *
164 * @param <T> Type of the field elements
165 * @param ex e cos(ω), first component of circular eccentricity vector
166 * @param ey e sin(ω), second component of circular eccentricity vector
167 * @param alphaV = V + ω true latitude argument (rad)
168 * @return the mean latitude argument.
169 */
170 public static <T extends CalculusFieldElement<T>> T trueToMean(final T ex, final T ey, final T alphaV) {
171 final T alphaE = trueToEccentric(ex, ey, alphaV);
172 return eccentricToMean(ex, ey, alphaE);
173 }
174
175 /**
176 * Computes the true latitude argument from the eccentric latitude argument.
177 *
178 * @param <T> Type of the field elements
179 * @param ex e cos(ω), first component of circular eccentricity vector
180 * @param ey e sin(ω), second component of circular eccentricity vector
181 * @param alphaM = M + ω mean latitude argument (rad)
182 * @return the true latitude argument.
183 */
184 public static <T extends CalculusFieldElement<T>> T meanToTrue(final T ex, final T ey, final T alphaM) {
185 final T alphaE = meanToEccentric(ex, ey, alphaM);
186 return eccentricToTrue(ex, ey, alphaE);
187 }
188
189 }