1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.measurements.modifiers;
18
19 import org.hipparchus.Field;
20 import org.hipparchus.analysis.UnivariateFunction;
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.Gradient;
25 import org.hipparchus.analysis.differentiation.GradientField;
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.hipparchus.util.FastMath;
28 import org.junit.jupiter.api.Assertions;
29 import org.junit.jupiter.api.BeforeAll;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.Utils;
32 import org.orekit.bodies.GeodeticPoint;
33 import org.orekit.data.DataContext;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.errors.OrekitMessages;
36 import org.orekit.estimation.measurements.AngularRaDec;
37 import org.orekit.estimation.measurements.EstimatedMeasurement;
38 import org.orekit.estimation.measurements.GroundStation;
39 import org.orekit.estimation.measurements.ObservableSatellite;
40 import org.orekit.frames.FieldTransform;
41 import org.orekit.frames.Frame;
42 import org.orekit.frames.FramesFactory;
43 import org.orekit.frames.TopocentricFrame;
44 import org.orekit.models.earth.ReferenceEllipsoid;
45 import org.orekit.orbits.CartesianOrbit;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.time.AbsoluteDate;
48 import org.orekit.time.FieldAbsoluteDate;
49 import org.orekit.time.TimeScalesFactory;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.IERSConventions;
52 import org.orekit.utils.PVCoordinates;
53 import org.orekit.utils.ParameterDriver;
54
55 import java.util.ArrayList;
56 import java.util.HashMap;
57 import java.util.List;
58 import java.util.Map;
59
60 public class AberrationModifierTest {
61
62 static GroundStation groundStation;
63
64 @BeforeAll
65 static void setup() {
66 Utils.setDataRoot("cr3bp:regular-data");
67
68
69 Frame fixedFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
70
71
72 GeodeticPoint stationLocation = new GeodeticPoint(FastMath.toRadians(-31.511673),
73 FastMath.toRadians(139.343695),
74 51.8);
75 TopocentricFrame stationFrame = new TopocentricFrame(ReferenceEllipsoid.getWgs84(fixedFrame),
76 stationLocation, "station");
77 groundStation = new GroundStation(stationFrame);
78
79
80 List<ParameterDriver> parameterDrivers = new ArrayList<>();
81 parameterDrivers.add(groundStation.getPrimeMeridianOffsetDriver());
82 parameterDrivers.add(groundStation.getPolarOffsetXDriver());
83 parameterDrivers.add(groundStation.getPolarOffsetYDriver());
84 for (ParameterDriver driver : parameterDrivers) {
85 driver.setReferenceDate(AbsoluteDate.ARBITRARY_EPOCH);
86 driver.setSelected(true);
87 }
88
89 }
90
91 @Test
92 void testAberration() {
93
94
95 Frame measurementFrame = FramesFactory.getGCRF();
96
97
98 double[] raDec = new double[]{-0.094838051874766, 0.157294151619018};
99
100
101 AbsoluteDate epoch = new AbsoluteDate(2022, 11, 8, 5, 59, 42.0, TimeScalesFactory.getUTC());
102
103
104 double[] proper = AberrationModifier.naturalToProper(raDec, groundStation, epoch, measurementFrame);
105
106
107 double[] expected = new double[]{-0.0947807299849455, 0.157333364635137};
108 Assertions.assertArrayEquals(expected, proper, 1e-10);
109
110
111 double[] natural = AberrationModifier.properToNatural(proper, groundStation, epoch, measurementFrame);
112 Assertions.assertArrayEquals(raDec, natural, 1e-12);
113
114 }
115
116 @FunctionalInterface
117 public interface AberrationTransform {
118 Gradient[] apply(Gradient[] raDecDS, FieldTransform<Gradient> stationToInertial, Frame obsFrame);
119 }
120
121 private void checkDerivative(AberrationTransform aberrationTransform,
122 int raDecIndex,
123 double gradient,
124 Gradient[] raDecDS,
125 GroundStation groundStation,
126 ParameterDriver driver,
127 Frame obsFrame,
128 AbsoluteDate epoch,
129 int nbParams,
130 Map<String, Integer> indices) {
131
132 FiniteDifferencesDifferentiator differentiator =
133 new FiniteDifferencesDifferentiator(3, 50.0 * driver.getScale());
134 UnivariateFunction uf = new UnivariateFunction() {
135
136
137
138 @Override
139 public double value(final double value) {
140 double current = driver.getValue();
141
142 driver.setValue(value);
143 FieldTransform<Gradient> stationToInertial =
144 groundStation.getOffsetToInertial(obsFrame, epoch, nbParams, indices);
145 Gradient[] naturalDS = aberrationTransform.apply(raDecDS, stationToInertial, obsFrame);
146
147 driver.setValue(current);
148 return naturalDS[raDecIndex].getValue();
149 }
150 };
151
152 DSFactory factory = new DSFactory(1, 1);
153 final DerivativeStructure dsParam = factory.variable(0, driver.getValue());
154 final DerivativeStructure dsValue = differentiator.differentiate(uf).value(dsParam);
155
156 double deriv = dsValue.getPartialDerivative(1);
157
158 Assertions.assertEquals(gradient, deriv, 1e-3 * FastMath.abs(gradient));
159 }
160
161 private void checkNaturalToProperDerivative(int raDecIndex,
162 double gradient,
163 Gradient[] raDecDS,
164 GroundStation groundStation,
165 ParameterDriver driver,
166 Frame obsFrame,
167 AbsoluteDate epoch,
168 int nbParams,
169 Map<String, Integer> indices) {
170
171 checkDerivative(AberrationModifier::fieldNaturalToProper, raDecIndex,
172 gradient, raDecDS, groundStation, driver,
173 obsFrame, epoch, nbParams, indices);
174 }
175
176 private void checkProperToNaturalDerivative(int raDecIndex,
177 double gradient,
178 Gradient[] raDecDS,
179 GroundStation groundStation,
180 ParameterDriver driver,
181 Frame obsFrame,
182 AbsoluteDate epoch,
183 int nbParams,
184 Map<String, Integer> indices) {
185
186 checkDerivative(AberrationModifier::fieldProperToNatural, raDecIndex,
187 gradient, raDecDS, groundStation, driver,
188 obsFrame, epoch, nbParams, indices);
189 }
190
191 @Test
192 void testFieldAberration() {
193
194
195 Frame measurementFrame = FramesFactory.getGCRF();
196
197
198 double[] raDec = new double[]{-0.094838051874766, 0.157294151619018};
199
200
201 AbsoluteDate epoch = new AbsoluteDate(2022, 11, 8, 5, 59, 42.0, TimeScalesFactory.getUTC());
202
203
204 List<ParameterDriver> parameterDrivers = new ArrayList<>();
205 parameterDrivers.add(groundStation.getPrimeMeridianOffsetDriver());
206 parameterDrivers.add(groundStation.getPolarOffsetXDriver());
207 parameterDrivers.add(groundStation.getPolarOffsetYDriver());
208
209
210 int nbParams = 6;
211 final Map<String, Integer> indices = new HashMap<>();
212
213 for (ParameterDriver driver : parameterDrivers) {
214 driver.setReferenceDate(epoch);
215 driver.setSelected(true);
216 indices.put(driver.getNameSpan(epoch), nbParams++);
217 }
218
219
220 final FieldTransform<Gradient> stationToInertial =
221 groundStation.getOffsetToInertial(measurementFrame, epoch, nbParams, indices);
222
223
224 final Field<Gradient> field = GradientField.getField(nbParams);
225 Gradient[] raDecDS = new Gradient[]{
226 field.getZero().add(raDec[0]),
227 field.getZero().add(raDec[1])
228 };
229 Gradient[] properDS = AberrationModifier.fieldNaturalToProper(raDecDS, stationToInertial, measurementFrame);
230 double[] proper = new double[]{properDS[0].getValue(), properDS[1].getValue()};
231
232
233 double[] expected = new double[]{-0.0947807299849455, 0.157333364635137};
234 Assertions.assertArrayEquals(expected, proper, 1e-10);
235
236
237 for (ParameterDriver driver : parameterDrivers) {
238 double raGradient = properDS[0].getGradient()[indices.get(driver.getNameSpan(epoch))];
239 checkNaturalToProperDerivative(0, raGradient, raDecDS, groundStation, driver,
240 measurementFrame, epoch, nbParams, indices);
241
242 double decGradient = properDS[1].getGradient()[indices.get(driver.getNameSpan(epoch))];
243 checkNaturalToProperDerivative(1, decGradient, raDecDS, groundStation, driver,
244 measurementFrame, epoch, nbParams, indices);
245 }
246
247
248 Gradient[] expectedDS = new Gradient[]{
249 field.getZero().add(expected[0]),
250 field.getZero().add(expected[1])
251 };
252 Gradient[] naturalDS = AberrationModifier.fieldProperToNatural(expectedDS, stationToInertial, measurementFrame);
253 double[] natural = new double[]{naturalDS[0].getValue(), naturalDS[1].getValue()};
254
255
256 Assertions.assertArrayEquals(raDec, natural, 1e-10);
257
258
259 for (ParameterDriver driver : parameterDrivers) {
260 double raGradient = naturalDS[0].getGradient()[indices.get(driver.getNameSpan(epoch))];
261 checkProperToNaturalDerivative(0, raGradient, expectedDS, groundStation, driver,
262 measurementFrame, epoch, nbParams, indices);
263
264 double decGradient = naturalDS[1].getGradient()[indices.get(driver.getNameSpan(epoch))];
265 checkProperToNaturalDerivative(1, decGradient, expectedDS, groundStation, driver,
266 measurementFrame, epoch, nbParams, indices);
267 }
268
269 }
270
271
272 @Test
273 void testAberrationModifier() {
274
275
276 AbsoluteDate epoch = new AbsoluteDate(2022, 11, 9, 3, 15, 18.454, TimeScalesFactory.getUTC());
277
278
279 Vector3D position = new Vector3D(-4350341.308092136, 2.5233509978715107E7, -8187957.452234574);
280 Vector3D velocity = new Vector3D(-2097.9025703889993, -1293.1199759574952, -2928.2553383447744);
281 PVCoordinates coordinates = new PVCoordinates(position, velocity);
282 CartesianOrbit orbit = new CartesianOrbit(coordinates, FramesFactory.getGCRF(), epoch, Constants.IERS2010_EARTH_MU);
283 SpacecraftState state = new SpacecraftState(orbit);
284
285
286 AngularRaDec raDec = defaultRaDec(FramesFactory.getGCRF(), epoch);
287
288
289 double[] estimatedRaDec = raDec
290 .estimate(0, 0, new SpacecraftState[]{state})
291 .getEstimatedValue();
292
293
294 AngularRaDec modifiedRaDec = defaultRaDec(FramesFactory.getGCRF(), epoch);
295 modifiedRaDec.addModifier(new AberrationModifier());
296
297
298 EstimatedMeasurement<AngularRaDec> modEstimated = modifiedRaDec
299 .estimate(0, 0, new SpacecraftState[]{state});
300 double[] estModRaDec = modEstimated.getEstimatedValue();
301
302
303 double[] unmodRaDec = AberrationModifier.naturalToProper(estModRaDec, groundStation, epoch, FramesFactory.getGCRF());
304 Assertions.assertArrayEquals(estimatedRaDec, unmodRaDec, 1e-12);
305
306 Assertions.assertEquals(1,
307 modEstimated.getAppliedEffects().entrySet().stream().
308 filter(e -> e.getKey().getEffectName().equals("aberration")).count());
309
310 }
311
312 @Test
313 public void testIssue1230() {
314
315
316 final AbsoluteDate epoch = new AbsoluteDate(2022, 11, 9, 3, 15, 18.454, TimeScalesFactory.getUTC());
317 final Frame frame = FramesFactory.getGCRF();
318 final Vector3D position = new Vector3D(-4350341.308092136, 2.5233509978715107E7, -8187957.452234574);
319 final Vector3D velocity = new Vector3D(-2097.9025703889993, -1293.1199759574952, -2928.2553383447744);
320 final PVCoordinates coordinates = new PVCoordinates(position, velocity);
321 final CartesianOrbit orbit = new CartesianOrbit(coordinates, frame, epoch, Constants.IERS2010_EARTH_MU);
322 final SpacecraftState state = new SpacecraftState(orbit);
323
324
325 final AngularRaDec raDec1 = defaultRaDec(frame, epoch);
326 final AngularRaDec raDec2 = defaultRaDec(frame, epoch);
327
328
329 final EstimatedMeasurement<AngularRaDec> estimated1 = raDec1.estimate(0, 0, new SpacecraftState[] { state });
330 final EstimatedMeasurement<AngularRaDec> estimated2 = raDec2.estimate(0, 0, new SpacecraftState[] { state });
331
332
333 final DataContext dataContext = DataContext.getDefault();
334
335
336 new AberrationModifier().modify(estimated1);
337 new AberrationModifier(dataContext).modify(estimated2);
338 Assertions.assertEquals(estimated1.getEstimatedValue()[0], estimated2.getEstimatedValue()[0]);
339 Assertions.assertEquals(estimated1.getEstimatedValue()[1], estimated2.getEstimatedValue()[1]);
340
341
342 new AberrationModifier().modifyWithoutDerivatives(estimated1);
343 new AberrationModifier(dataContext).modifyWithoutDerivatives(estimated2);
344 Assertions.assertEquals(estimated1.getEstimatedValue()[0], estimated2.getEstimatedValue()[0]);
345 Assertions.assertEquals(estimated1.getEstimatedValue()[1], estimated2.getEstimatedValue()[1]);
346
347
348 final double[] proper1 = AberrationModifier.naturalToProper(raDec1.getObservedValue(), groundStation, epoch, frame);
349 final double[] natural1 = AberrationModifier.properToNatural(proper1, groundStation, epoch, frame);
350 final double[] proper2 = AberrationModifier.naturalToProper(raDec2.getObservedValue(), groundStation, epoch, frame, dataContext);
351 final double[] natural2 = AberrationModifier.properToNatural(proper2, groundStation, epoch, frame, dataContext);
352 Assertions.assertEquals(natural1[0], natural2[0]);
353 Assertions.assertEquals(natural1[1], natural2[1]);
354
355
356 final Field<Gradient> field = GradientField.getField(6);
357 final Gradient[] raDecG = new Gradient[] { field.getZero().add(raDec1.getObservedValue()[0]),
358 field.getZero().add(raDec1.getObservedValue()[1]) };
359 final FieldTransform<Gradient> stationToInertial =
360 groundStation.getBaseFrame().getTransformTo(frame, new FieldAbsoluteDate<>(field, epoch));
361
362 final Gradient[] proper1G = AberrationModifier.fieldNaturalToProper(raDecG, stationToInertial, frame);
363 final Gradient[] natural1G = AberrationModifier.fieldProperToNatural(proper1G, stationToInertial, frame);
364 final Gradient[] proper2G = AberrationModifier.fieldNaturalToProper(raDecG, stationToInertial, frame, dataContext);
365 final Gradient[] natural2G = AberrationModifier.fieldProperToNatural(proper2G, stationToInertial, frame, dataContext);
366
367 Assertions.assertEquals(natural1G[0].getValue(), natural2G[0].getValue());
368 Assertions.assertEquals(natural1G[1].getValue(), natural2G[1].getValue());
369
370 }
371
372 @Test
373 public void testExceptionIfNonInertialFrame() {
374
375 final AbsoluteDate epoch = new AbsoluteDate(2022, 11, 9, 3, 15, 18.454, TimeScalesFactory.getUTC());
376 final double[] raDec1 = new double[] { 1.0, 1.0 };
377 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
378
379
380
381 final OrekitException exceptionNToP = Assertions.assertThrows(OrekitException.class,
382 () -> AberrationModifier.naturalToProper(raDec1,
383 groundStation,
384 epoch, itrf));
385 final OrekitException exceptionPToN = Assertions.assertThrows(OrekitException.class,
386 () -> AberrationModifier.properToNatural(raDec1,
387 groundStation,
388 epoch, itrf));
389
390 Assertions.assertEquals(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, exceptionNToP.getSpecifier());
391 Assertions.assertEquals(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, exceptionPToN.getSpecifier());
392 }
393
394 private static AngularRaDec defaultRaDec(Frame frame, AbsoluteDate date) {
395 return new AngularRaDec(groundStation, frame, date, new double[]{0.0, 0.0},
396 new double[]{1.0, 1.0}, new double[]{1.0, 1.0}, new ObservableSatellite(0));
397 }
398 }