1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.bodies;
18
19 import org.hipparchus.analysis.UnivariateFunction;
20 import org.hipparchus.analysis.UnivariateVectorFunction;
21 import org.hipparchus.analysis.differentiation.DSFactory;
22 import org.hipparchus.analysis.differentiation.DerivativeStructure;
23 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
24 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
25 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
26 import org.hipparchus.geometry.euclidean.threed.Rotation;
27 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.util.Binary64;
30 import org.hipparchus.util.Binary64Field;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathUtils;
33 import org.hipparchus.util.Precision;
34 import org.junit.jupiter.api.Assertions;
35 import org.junit.jupiter.api.BeforeEach;
36 import org.junit.jupiter.api.Test;
37 import org.orekit.Utils;
38 import org.orekit.bodies.IAUPoleFactory.OldIAUPole;
39 import org.orekit.bodies.JPLEphemeridesLoader.EphemerisType;
40 import org.orekit.data.DataContext;
41 import org.orekit.time.AbsoluteDate;
42 import org.orekit.time.FieldAbsoluteDate;
43 import org.orekit.time.TimeScale;
44 import org.orekit.time.TimeScales;
45 import org.orekit.time.TimeScalesFactory;
46 import org.orekit.utils.Constants;
47
48 import java.io.BufferedReader;
49 import java.io.IOException;
50 import java.io.InputStream;
51 import java.io.InputStreamReader;
52 import java.io.UnsupportedEncodingException;
53 import java.nio.charset.StandardCharsets;
54
55 public class PredefinedIAUPolesTest {
56
57 @Test
58 void testGCRFAligned() throws UnsupportedEncodingException, IOException {
59 IAUPole iauPole = PredefinedIAUPoles.getIAUPole(EphemerisType.SOLAR_SYSTEM_BARYCENTER, timeScales);
60 Vector3D pole = iauPole.getPole(AbsoluteDate.J2000_EPOCH);
61 double w = iauPole.getPrimeMeridianAngle(AbsoluteDate.J2000_EPOCH.shiftedBy(3600.0));
62 Assertions.assertEquals(0, Vector3D.distance(pole, Vector3D.PLUS_K), 1.0e-15);
63 Assertions.assertEquals(0.0, w, 1.0e-15);
64 }
65
66 @Test
67 void testSun() throws UnsupportedEncodingException, IOException {
68 IAUPole iauPole = PredefinedIAUPoles.getIAUPole(EphemerisType.SUN, timeScales);
69 Vector3D pole = iauPole.getPole(AbsoluteDate.J2000_EPOCH);
70 final double alphaRef = FastMath.toRadians(286.13);
71 final double deltaRef = FastMath.toRadians(63.87);
72 final double wRef = FastMath.toRadians(84.176);
73 final double rateRef = FastMath.toRadians(14.1844000);
74 double w = iauPole.getPrimeMeridianAngle(new AbsoluteDate(AbsoluteDate.J2000_EPOCH, 3600.0,
75 TimeScalesFactory.getTDB()));
76 Assertions.assertEquals(alphaRef, MathUtils.normalizeAngle(pole.getAlpha(), alphaRef), 1.0e-15);
77 Assertions.assertEquals(deltaRef, pole.getDelta(), 1.0e-15);
78 Assertions.assertEquals(wRef + rateRef / 24.0, w, 1.0e-15);
79 }
80
81 @Test
82 void testNaif() throws UnsupportedEncodingException, IOException {
83 final TimeScale tdb = TimeScalesFactory.getTDB();
84 final InputStream inEntry = getClass().getResourceAsStream("/naif/IAU-pole-NAIF.txt");
85 try (BufferedReader reader = new BufferedReader(new InputStreamReader(inEntry, StandardCharsets.UTF_8))) {
86 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
87 line = line.trim();
88 if (!line.isEmpty() && !line.startsWith("#")) {
89
90
91 String[] fields = line.split("\\s+");
92 final AbsoluteDate date1 = new AbsoluteDate(fields[0], tdb);
93 final AbsoluteDate date2 = new AbsoluteDate(AbsoluteDate.J2000_EPOCH,
94 Double.parseDouble(fields[1]),
95 tdb);
96 final EphemerisType type = EphemerisType.valueOf(fields[2]);
97 final double alphaRef = Double.parseDouble(fields[3]);
98 final double deltaRef = Double.parseDouble(fields[4]);
99 final double wRef = Double.parseDouble(fields[5]);
100 final double[][] m = new double[3][3];
101 int index = 6;
102 for (int i = 0; i < 3; ++i) {
103 for (int j = 0; j < 3; ++j) {
104
105
106 m[j][i] = Double.parseDouble(fields[index++]);
107 }
108 }
109 Rotation rRef = new Rotation(m, 1.0e-10);
110
111
112 IAUPole iauPole = PredefinedIAUPoles.getIAUPole(type, timeScales);
113 Vector3D pole = iauPole.getPole(date2);
114 double w = iauPole.getPrimeMeridianAngle(date2);
115 Assertions.assertEquals(0.0, date2.durationFrom(date1), 8.0e-5);
116 Assertions.assertEquals(alphaRef, MathUtils.normalizeAngle(pole.getAlpha(), alphaRef), 1.8e-15);
117 Assertions.assertEquals(deltaRef, pole.getDelta(), 2.4e-13);
118 Assertions.assertEquals(wRef, MathUtils.normalizeAngle(w, wRef), 2.5e-12);
119
120
121 Vector3D qNode = Vector3D.crossProduct(Vector3D.PLUS_K, pole);
122 if (qNode.getNormSq() < Precision.SAFE_MIN) {
123 qNode = Vector3D.PLUS_I;
124 }
125 final Rotation rotation = new Rotation(Vector3D.PLUS_K, wRef, RotationConvention.FRAME_TRANSFORM).
126 applyTo(new Rotation(pole, qNode, Vector3D.PLUS_K, Vector3D.PLUS_I));
127 Assertions.assertEquals(0.0, Rotation.distance(rRef, rotation), 1.9e-15);
128
129 }
130 }
131 }
132 }
133
134 @Test
135 void testVersus80Implementation() {
136 for (EphemerisType body : EphemerisType.values()) {
137 IAUPole newPole = PredefinedIAUPoles.getIAUPole(body, timeScales);
138 OldIAUPole oldPole = IAUPoleFactory.getIAUPole(body);
139 for (double dt = 0; dt < Constants.JULIAN_YEAR; dt += 3600) {
140 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
141 Assertions.assertEquals(0, Vector3D.angle(newPole.getPole(date), oldPole.getPole(date)), 1.0e-20);
142 Assertions.assertEquals(oldPole.getPrimeMeridianAngle(date), newPole.getPrimeMeridianAngle(date), 5.0e-13);
143 }
144 }
145
146 }
147
148 @Test
149 void testFieldConsistency() {
150 for (IAUPole iaupole : PredefinedIAUPoles.values(timeScales)) {
151 for (double dt = 0; dt < Constants.JULIAN_YEAR; dt += 3600) {
152 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
153 final FieldAbsoluteDate<Binary64> date64 = new FieldAbsoluteDate<>(Binary64Field.getInstance(), date);
154 Assertions.assertEquals(0, Vector3D.angle(iaupole.getPole(date), iaupole.getPole(date64).toVector3D()), 2.0e-15);
155 Assertions.assertEquals(iaupole.getPrimeMeridianAngle(date), iaupole.getPrimeMeridianAngle(date64).getReal(), 1.0e-12);
156 }
157 }
158
159 }
160
161 @Test
162 void testDerivatives() {
163 final DSFactory factory = new DSFactory(1, 1);
164 final AbsoluteDate ref = AbsoluteDate.J2000_EPOCH;
165 final FieldAbsoluteDate<DerivativeStructure> refDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), ref);
166 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 60.0);
167 for (final IAUPole iaupole : PredefinedIAUPoles.values(timeScales)) {
168
169 UnivariateDifferentiableVectorFunction dPole = differentiator.differentiate(new UnivariateVectorFunction() {
170 @Override
171 public double[] value(double t) {
172 return iaupole.getPole(ref.shiftedBy(t)).toArray();
173 }
174 });
175 UnivariateDifferentiableFunction dMeridian = differentiator.differentiate(new UnivariateFunction() {
176 @Override
177 public double value(double t) {
178 return iaupole.getPrimeMeridianAngle(ref.shiftedBy(t));
179 }
180 });
181
182 for (double dt = 0; dt < Constants.JULIAN_YEAR; dt += 3600) {
183
184 final DerivativeStructure dtDS = factory.variable(0, dt);
185
186 final DerivativeStructure[] refPole = dPole.value(dtDS);
187 final DerivativeStructure[] fieldPole = iaupole.getPole(refDS.shiftedBy(dtDS)).toArray();
188 for (int i = 0; i < 3; ++i) {
189 Assertions.assertEquals(refPole[i].getValue(), fieldPole[i].getValue(), 2.0e-15);
190 Assertions.assertEquals(refPole[i].getPartialDerivative(1), fieldPole[i].getPartialDerivative(1), 4.0e-17);
191 }
192
193 final DerivativeStructure refMeridian = dMeridian.value(dtDS);
194 final DerivativeStructure fieldMeridian = iaupole.getPrimeMeridianAngle(refDS.shiftedBy(dtDS));
195 Assertions.assertEquals(refMeridian.getValue(), fieldMeridian.getValue(), 4.0e-12);
196 Assertions.assertEquals(refMeridian.getPartialDerivative(1), fieldMeridian.getPartialDerivative(1), 9.0e-14);
197
198 }
199 }
200
201 }
202
203 private TimeScales timeScales;
204
205 @BeforeEach
206 public void setUp() {
207 Utils.setDataRoot("regular-data");
208 timeScales = DataContext.getDefault().getTimeScales();
209 }
210
211 }
212