1 /* Copyright 2002-2025 CS GROUP
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.bodies;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.orekit.time.AbsoluteDate;
22 import org.orekit.time.TimeScalesFactory;
23 import org.orekit.utils.Constants;
24
25 import java.io.Serializable;
26
27 /** Old Factory class for IAU poles, replaced as of Orekit 9.0 by a new implementation.
28 * <p>The pole models provided here come from the <a
29 * href="http://astropedia.astrogeology.usgs.gov/alfresco/d/d/workspace/SpacesStore/28fd9e81-1964-44d6-a58b-fbbf61e64e15/WGCCRE2009reprint.pdf">
30 * 2009 report</a> and the <a href="http://astropedia.astrogeology.usgs.gov/alfresco/d/d/workspace/SpacesStore/04d348b0-eb2b-46a2-abe9-6effacb37763/WGCCRE-Erratum-2011reprint.pdf">
31 * 2011 erratum</a> of the IAU/IAG Working Group on Cartographic Coordinates
32 * and Rotational Elements of the Planets and Satellites (WGCCRE). Note that these value
33 * differ from earliest reports (before 2005).
34 *</p>
35 * @author Luc Maisonobe
36 * @since 5.1
37 */
38 class IAUPoleFactory {
39
40 interface OldIAUPole extends Serializable {
41
42 /** Get the body North pole direction in ICRF frame.
43 * @param date current date
44 * @return body North pole direction in ICRF frame
45 */
46 Vector3D getPole(AbsoluteDate date);
47
48 /** Get the prime meridian angle.
49 * <p>
50 * The prime meridian angle is the angle between the Q node and the
51 * prime meridian. represents the body rotation.
52 * </p>
53 * @param date current date
54 * @return prime meridian vector
55 */
56 double getPrimeMeridianAngle(AbsoluteDate date);
57
58 }
59
60 /** Private constructor.
61 * <p>This class is a utility class, it should neither have a public
62 * nor a default constructor. This private constructor prevents
63 * the compiler from generating one automatically.</p>
64 */
65 private IAUPoleFactory() {
66 }
67
68 /** Get an IAU pole.
69 * @param body body for which the pole is requested
70 * @return IAU pole for the body, or dummy GCRF aligned pole
71 * for barycenters
72 */
73 public static OldIAUPole getIAUPole(final JPLEphemeridesLoader.EphemerisType body) {
74 switch (body) {
75 case SUN:
76 return new OldIAUPole() {
77
78 /** Serializable UID. */
79 private static final long serialVersionUID = 5715331729495237139L;
80
81 /** {@inheritDoc }*/
82 public Vector3D getPole(final AbsoluteDate date) {
83 return new Vector3D(FastMath.toRadians(286.13),
84 FastMath.toRadians(63.87));
85 }
86
87 /** {@inheritDoc }*/
88 public double getPrimeMeridianAngle(final AbsoluteDate date) {
89 return FastMath.toRadians(84.176 + 14.1844000 * d(date));
90 }
91
92 };
93 case MERCURY:
94 return new OldIAUPole() {
95
96 /** Serializable UID. */
97 private static final long serialVersionUID = -5769710119654037007L;
98
99 /** {@inheritDoc }*/
100 public Vector3D getPole(final AbsoluteDate date) {
101 final double t = t(date);
102 return new Vector3D(FastMath.toRadians(281.0097 - 0.0328 * t),
103 FastMath.toRadians( 61.4143 - 0.0049 * t));
104 }
105
106 /** {@inheritDoc }*/
107 public double getPrimeMeridianAngle(final AbsoluteDate date) {
108 final double[] m = computeMi(date);
109 return FastMath.toRadians(329.5469 + 6.1385025 * d(date) +
110 0.00993822 * FastMath.sin(m[0]) -
111 0.00104581 * FastMath.sin(m[1]) -
112 0.00010280 * FastMath.sin(m[2]) -
113 0.00002364 * FastMath.sin(m[3]) -
114 0.00000532 * FastMath.sin(m[4]));
115 }
116
117 /** Compute the Mercury angles M<sub>i</sub>.
118 * @param date date
119 * @return array of Mercury angles, with M<sub>i</sub> stored at index i-1
120 */
121 private double[] computeMi(final AbsoluteDate date) {
122 final double d = d(date);
123 return new double[] {
124 FastMath.toRadians(174.791086 + 4.092335 * d), // M1
125 FastMath.toRadians(349.582171 + 8.184670 * d), // M2
126 FastMath.toRadians(164.373257 + 12.277005 * d), // M3
127 FastMath.toRadians(339.164343 + 16.369340 * d), // M4
128 FastMath.toRadians(153.955429 + 20.461675 * d), // M5
129 };
130 }
131 };
132 case VENUS:
133 return new OldIAUPole() {
134
135 /** Serializable UID. */
136 private static final long serialVersionUID = 7030506277976648896L;
137
138 /** {@inheritDoc }*/
139 public Vector3D getPole(final AbsoluteDate date) {
140 return new Vector3D(FastMath.toRadians(272.76),
141 FastMath.toRadians(67.16));
142 }
143
144 /** {@inheritDoc }*/
145 public double getPrimeMeridianAngle(final AbsoluteDate date) {
146 return FastMath.toRadians(160.20 - 1.4813688 * d(date));
147 }
148
149 };
150 case EARTH:
151 return new OldIAUPole() {
152
153 /** Serializable UID. */
154 private static final long serialVersionUID = 6912325697192667056L;
155
156 /** {@inheritDoc }*/
157 public Vector3D getPole(final AbsoluteDate date) {
158 final double t = t(date);
159 return new Vector3D(FastMath.toRadians( 0.00 - 0.641 * t),
160 FastMath.toRadians(90.00 - 0.557 * t));
161 }
162
163 /** {@inheritDoc }*/
164 public double getPrimeMeridianAngle(final AbsoluteDate date) {
165 return FastMath.toRadians(190.147 + 360.9856235 * d(date));
166 }
167
168 };
169 case MOON:
170 return new OldIAUPole() {
171
172
173 /** Serializable UID. */
174 private static final long serialVersionUID = -1310155975084976571L;
175
176 /** {@inheritDoc }*/
177 public Vector3D getPole(final AbsoluteDate date) {
178 final double[] e = computeEi(date);
179 final double t = t(date);
180 return new Vector3D(FastMath.toRadians(269.9949 + 0.0031 * t - 3.8787 * FastMath.sin(e[0]) -
181 0.1204 * FastMath.sin(e[1]) + 0.0700 * FastMath.sin(e[2]) -
182 0.0172 * FastMath.sin(e[3]) + 0.0072 * FastMath.sin(e[5]) -
183 0.0052 * FastMath.sin(e[9]) + 0.0043 * FastMath.sin(e[12])),
184 FastMath.toRadians( 66.5392 + 0.0130 * t + 1.5419 * FastMath.cos(e[0]) +
185 0.0239 * FastMath.cos(e[1]) - 0.0278 * FastMath.cos(e[2]) +
186 0.0068 * FastMath.cos(e[3]) - 0.0029 * FastMath.cos(e[5]) +
187 0.0009 * FastMath.cos(e[6]) + 0.0008 * FastMath.cos(e[9]) -
188 0.0009 * FastMath.cos(e[12])));
189 }
190
191 /** {@inheritDoc }*/
192 public double getPrimeMeridianAngle(final AbsoluteDate date) {
193 final double[] e = computeEi(date);
194 final double d = d(date);
195 return FastMath.toRadians(38.3213 + (13.17635815 - 1.4e-12 * d) * d + 3.5610 * FastMath.sin(e[0]) +
196 0.1208 * FastMath.sin(e[1]) - 0.0642 * FastMath.sin(e[2]) +
197 0.0158 * FastMath.sin(e[3]) + 0.0252 * FastMath.sin(e[4]) -
198 0.0066 * FastMath.sin(e[5]) - 0.0047 * FastMath.sin(e[6]) -
199 0.0046 * FastMath.sin(e[7]) + 0.0028 * FastMath.sin(e[8]) +
200 0.0052 * FastMath.sin(e[9]) + 0.0040 * FastMath.sin(e[10]) +
201 0.0019 * FastMath.sin(e[11]) - 0.0044 * FastMath.sin(e[12]));
202 }
203
204 /** Compute the Moon angles E<sub>i</sub>.
205 * @param date date
206 * @return array of Moon angles, with E<sub>i</sub> stored at index i-1
207 */
208 private double[] computeEi(final AbsoluteDate date) {
209 final double d = d(date);
210 return new double[] {
211 FastMath.toRadians(125.045 - 0.0529921 * d), // E1
212 FastMath.toRadians(250.089 - 0.1059842 * d), // E2
213 FastMath.toRadians(260.008 + 13.0120009 * d), // E3
214 FastMath.toRadians(176.625 + 13.3407154 * d), // E4
215 FastMath.toRadians(357.529 + 0.9856003 * d), // E5
216 FastMath.toRadians(311.589 + 26.4057084 * d), // E6
217 FastMath.toRadians(134.963 + 13.0649930 * d), // E7
218 FastMath.toRadians(276.617 + 0.3287146 * d), // E8
219 FastMath.toRadians( 34.226 + 1.7484877 * d), // E9
220 FastMath.toRadians( 15.134 - 0.1589763 * d), // E10
221 FastMath.toRadians(119.743 + 0.0036096 * d), // E11
222 FastMath.toRadians(239.961 + 0.1643573 * d), // E12
223 FastMath.toRadians( 25.053 + 12.9590088 * d) // E13
224 };
225 }
226
227 };
228 case MARS:
229 return new OldIAUPole() {
230
231 /** Serializable UID. */
232 private static final long serialVersionUID = 1471983418540015411L;
233
234 /** {@inheritDoc }*/
235 public Vector3D getPole(final AbsoluteDate date) {
236 final double t = t(date);
237 return new Vector3D(FastMath.toRadians(317.68143 - 0.1061 * t),
238 FastMath.toRadians( 52.88650 - 0.0609 * t));
239 }
240
241 /** {@inheritDoc }*/
242 public double getPrimeMeridianAngle(final AbsoluteDate date) {
243 return FastMath.toRadians(176.630 + 350.89198226 * d(date));
244 }
245
246 };
247 case JUPITER:
248 return new OldIAUPole() {
249
250 /** Serializable UID. */
251 private static final long serialVersionUID = 6959753758673537524L;
252
253 /** {@inheritDoc }*/
254 public Vector3D getPole(final AbsoluteDate date) {
255
256 final double t = t(date);
257 final double ja = FastMath.toRadians( 99.360714 + 4850.4046 * t);
258 final double jb = FastMath.toRadians(175.895369 + 1191.9605 * t);
259 final double jc = FastMath.toRadians(300.323162 + 262.5475 * t);
260 final double jd = FastMath.toRadians(114.012305 + 6070.2476 * t);
261 final double je = FastMath.toRadians( 49.511251 + 64.3000 * t);
262
263 return new Vector3D(FastMath.toRadians(268.056595 - 0.006499 * t +
264 0.000117 * FastMath.sin(ja) +
265 0.000938 * FastMath.sin(jb) +
266 0.001432 * FastMath.sin(jc) +
267 0.000030 * FastMath.sin(jd) +
268 0.002150 * FastMath.sin(je)),
269 FastMath.toRadians( 64.495303 + 0.002413 * t +
270 0.000050 * FastMath.cos(ja) +
271 0.000404 * FastMath.cos(jb) +
272 0.000617 * FastMath.cos(jc) -
273 0.000013 * FastMath.cos(jd) +
274 0.000926 * FastMath.cos(je)));
275 }
276
277 /** {@inheritDoc }*/
278 public double getPrimeMeridianAngle(final AbsoluteDate date) {
279 return FastMath.toRadians(284.95 + 870.5360000 * d(date));
280 }
281
282 };
283 case SATURN:
284 return new OldIAUPole() {
285
286 /** Serializable UID. */
287 private static final long serialVersionUID = -1082211873912149774L;
288
289 /** {@inheritDoc }*/
290 public Vector3D getPole(final AbsoluteDate date) {
291 final double t = t(date);
292 return new Vector3D(FastMath.toRadians(40.589 - 0.036 * t),
293 FastMath.toRadians(83.537 - 0.004 * t));
294 }
295
296 /** {@inheritDoc }*/
297 public double getPrimeMeridianAngle(final AbsoluteDate date) {
298 return FastMath.toRadians(38.90 + 810.7939024 * d(date));
299 }
300
301 };
302 case URANUS:
303 return new OldIAUPole() {
304
305 /** Serializable UID. */
306 private static final long serialVersionUID = 362792230470085154L;
307
308 /** {@inheritDoc }*/
309 public Vector3D getPole(final AbsoluteDate date) {
310 return new Vector3D(FastMath.toRadians(257.311),
311 FastMath.toRadians(-15.175));
312 }
313
314 /** {@inheritDoc }*/
315 public double getPrimeMeridianAngle(final AbsoluteDate date) {
316 return FastMath.toRadians(203.81 - 501.1600928 * d(date));
317 }
318
319 };
320 case NEPTUNE:
321 return new OldIAUPole() {
322
323 /** Serializable UID. */
324 private static final long serialVersionUID = 560614555734665287L;
325
326 /** {@inheritDoc }*/
327 public Vector3D getPole(final AbsoluteDate date) {
328 final double n = FastMath.toRadians(357.85 + 52.316 * t(date));
329 return new Vector3D(FastMath.toRadians(299.36 + 0.70 * FastMath.sin(n)),
330 FastMath.toRadians( 43.46 - 0.51 * FastMath.cos(n)));
331 }
332
333 /** {@inheritDoc }*/
334 public double getPrimeMeridianAngle(final AbsoluteDate date) {
335 final double n = FastMath.toRadians(357.85 + 52.316 * t(date));
336 return FastMath.toRadians(253.18 + 536.3128492 * d(date) - 0.48 * FastMath.sin(n));
337 }
338
339 };
340 case PLUTO:
341 return new OldIAUPole() {
342
343 /** Serializable UID. */
344 private static final long serialVersionUID = -1277113129327018062L;
345
346 /** {@inheritDoc }*/
347 public Vector3D getPole(final AbsoluteDate date) {
348 return new Vector3D(FastMath.toRadians(132.993),
349 FastMath.toRadians(-6.163));
350 }
351
352 /** {@inheritDoc }*/
353 public double getPrimeMeridianAngle(final AbsoluteDate date) {
354 return FastMath.toRadians(302.695 + 56.3625225 * d(date));
355 }
356
357 };
358 default:
359 return new GCRFAligned();
360 }
361 }
362
363 /** Compute the interval in julian centuries from standard epoch.
364 * @param date date
365 * @return interval between date and standard epoch in julian centuries
366 */
367 private static double t(final AbsoluteDate date) {
368 return date.offsetFrom(AbsoluteDate.J2000_EPOCH, TimeScalesFactory.getTDB()) / Constants.JULIAN_CENTURY;
369 }
370
371 /** Compute the interval in julian days from standard epoch.
372 * @param date date
373 * @return interval between date and standard epoch in julian days
374 */
375 private static double d(final AbsoluteDate date) {
376 return date.offsetFrom(AbsoluteDate.J2000_EPOCH, TimeScalesFactory.getTDB()) / Constants.JULIAN_DAY;
377 }
378
379 /** Default IAUPole implementation for barycenters.
380 * <p>
381 * This implementation defines directions such that the inertially oriented and body
382 * oriented frames are identical and aligned with GCRF. It is used for example
383 * to define the ICRF.
384 * </p>
385 */
386 private static class GCRFAligned implements OldIAUPole {
387
388 /** Serializable UID. */
389 private static final long serialVersionUID = 20130327L;
390
391 /** {@inheritDoc} */
392 public Vector3D getPole(final AbsoluteDate date) {
393 return Vector3D.PLUS_K;
394 }
395
396 /** {@inheritDoc} */
397 public double getPrimeMeridianAngle(final AbsoluteDate date) {
398 return 0;
399 }
400
401 }
402
403 }