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 longitude arguments used by {@link FieldEquinoctialOrbit}.
27 * @author Romain Serra
28 * @see FieldEquinoctialOrbit
29 * @since 12.1
30 */
31 public class FieldEquinoctialLongitudeArgumentUtility {
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 FieldEquinoctialLongitudeArgumentUtility() {
41 // nothing here (utils class)
42 }
43
44 /**
45 * Computes the true longitude argument from the eccentric longitude argument.
46 *
47 * @param <T> Type of the field elements
48 * @param ex e cos(ω), first component of eccentricity vector
49 * @param ey e sin(ω), second component of eccentricity vector
50 * @param lE = E + ω + Ω eccentric longitude argument (rad)
51 * @return the true longitude argument.
52 */
53 public static <T extends CalculusFieldElement<T>> T eccentricToTrue(final T ex, final T ey, final T lE) {
54 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
55 final FieldSinCos<T> scLE = FastMath.sinCos(lE);
56 final T cosLE = scLE.cos();
57 final T sinLE = scLE.sin();
58 final T num = ex.multiply(sinLE).subtract(ey.multiply(cosLE));
59 final T den = epsilon.add(1).subtract(ex.multiply(cosLE)).subtract(ey.multiply(sinLE));
60 return lE.add(eccentricAndTrueAtan(num, den));
61 }
62
63 /**
64 * Computes the eccentric longitude argument from the true longitude argument.
65 *
66 * @param <T> Type of the field elements
67 * @param ex e cos(ω), first component of eccentricity vector
68 * @param ey e sin(ω), second component of eccentricity vector
69 * @param lV = V + ω + Ω true longitude argument (rad)
70 * @return the eccentric longitude argument.
71 */
72 public static <T extends CalculusFieldElement<T>> T trueToEccentric(final T ex, final T ey, final T lV) {
73 final T epsilon = eccentricAndTrueEpsilon(ex, ey);
74 final FieldSinCos<T> scLv = FastMath.sinCos(lV);
75 final T cosLv = scLv.cos();
76 final T sinLv = scLv.sin();
77 final T num = ey.multiply(cosLv).subtract(ex.multiply(sinLv));
78 final T den = epsilon.add(1).add(ex.multiply(cosLv)).add(ey.multiply(sinLv));
79 return lV.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 eccentricity vector
87 * @param ey e sin(ω), second component of 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 longitude argument from the mean longitude argument.
108 *
109 * @param <T> Type of the field elements
110 * @param ex e cos(ω), first component of eccentricity vector
111 * @param ey e sin(ω), second component of eccentricity vector
112 * @param lM = M + ω + Ω mean longitude argument (rad)
113 * @return the eccentric longitude argument.
114 */
115 public static <T extends CalculusFieldElement<T>> T meanToEccentric(final T ex, final T ey, final T lM) {
116 // Generalization of Kepler equation to equinoctial parameters
117 // with lE = PA + RAAN + E and
118 // lM = PA + RAAN + M = lE - ex.sin(lE) + ey.cos(lE)
119 T lE = lM;
120 T shift;
121 T lEmlM = lM.getField().getZero();
122 boolean hasConverged;
123 int iter = 0;
124 do {
125 final FieldSinCos<T> scLE = FastMath.sinCos(lE);
126 final T f2 = ex.multiply(scLE.sin()).subtract(ey.multiply(scLE.cos()));
127 final T f1 = ex.multiply(scLE.cos()).add(ey.multiply(scLE.sin())).negate().add(1);
128 final T f0 = lEmlM.subtract(f2);
129
130 final T f12 = f1.multiply(2.0);
131 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
132
133 lEmlM = lEmlM.subtract(shift);
134 lE = lM.add(lEmlM);
135
136 hasConverged = FastMath.abs(shift.getReal()) <= TOLERANCE_CONVERGENCE;
137 } while (++iter < MAXIMUM_ITERATION && !hasConverged);
138
139 if (!hasConverged) {
140 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECCENTRIC_LONGITUDE_ARGUMENT, iter);
141 }
142 return lE;
143
144 }
145
146 /**
147 * Computes the mean longitude argument from the eccentric longitude argument.
148 *
149 * @param <T> Type of the field elements
150 * @param ex e cos(ω), first component of eccentricity vector
151 * @param ey e sin(ω), second component of eccentricity vector
152 * @param lE = E + ω + Ω mean longitude argument (rad)
153 * @return the mean longitude argument.
154 */
155 public static <T extends CalculusFieldElement<T>> T eccentricToMean(final T ex, final T ey, final T lE) {
156 final FieldSinCos<T> scLE = FastMath.sinCos(lE);
157 return lE.subtract(ex.multiply(scLE.sin())).add(ey.multiply(scLE.cos()));
158 }
159
160 /**
161 * Computes the mean longitude argument from the eccentric longitude argument.
162 *
163 * @param <T> Type of the field elements
164 * @param ex e cos(ω), first component of eccentricity vector
165 * @param ey e sin(ω), second component of eccentricity vector
166 * @param lV = V + ω + Ω true longitude argument (rad)
167 * @return the mean longitude argument.
168 */
169 public static <T extends CalculusFieldElement<T>> T trueToMean(final T ex, final T ey, final T lV) {
170 final T alphaE = trueToEccentric(ex, ey, lV);
171 return eccentricToMean(ex, ey, alphaE);
172 }
173
174 /**
175 * Computes the true longitude argument from the eccentric longitude argument.
176 *
177 * @param <T> Type of the field elements
178 * @param ex e cos(ω), first component of eccentricity vector
179 * @param ey e sin(ω), second component of eccentricity vector
180 * @param lM = M + ω + Ω mean longitude argument (rad)
181 * @return the true longitude argument.
182 */
183 public static <T extends CalculusFieldElement<T>> T meanToTrue(final T ex, final T ey, final T lM) {
184 final T alphaE = meanToEccentric(ex, ey, lM);
185 return eccentricToTrue(ex, ey, alphaE);
186 }
187
188 }