1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.measurements;
18
19 import java.util.List;
20
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.hipparchus.stat.descriptive.StreamingStatistics;
23 import org.hipparchus.stat.descriptive.rank.Median;
24 import org.hipparchus.util.FastMath;
25 import org.junit.jupiter.api.Assertions;
26 import org.junit.jupiter.api.Test;
27 import org.orekit.Utils;
28 import org.orekit.bodies.GeodeticPoint;
29 import org.orekit.bodies.OneAxisEllipsoid;
30 import org.orekit.estimation.Context;
31 import org.orekit.estimation.EstimationTestUtils;
32 import org.orekit.estimation.measurements.generation.AngularRaDecBuilder;
33 import org.orekit.frames.Frame;
34 import org.orekit.frames.FramesFactory;
35 import org.orekit.frames.ITRFVersion;
36 import org.orekit.frames.StaticTransform;
37 import org.orekit.frames.TopocentricFrame;
38 import org.orekit.orbits.CartesianOrbit;
39 import org.orekit.orbits.OrbitType;
40 import org.orekit.orbits.PositionAngleType;
41 import org.orekit.propagation.Propagator;
42 import org.orekit.propagation.SpacecraftState;
43 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
44 import org.orekit.time.AbsoluteDate;
45 import org.orekit.time.DateComponents;
46 import org.orekit.time.TimeScalesFactory;
47 import org.orekit.utils.Constants;
48 import org.orekit.utils.Differentiation;
49 import org.orekit.utils.IERSConventions;
50 import org.orekit.utils.PVCoordinates;
51 import org.orekit.utils.ParameterDriver;
52 import org.orekit.utils.ParameterFunction;
53
54 class AngularRaDecTest {
55
56
57
58
59 @Test
60 void testBug473OnValues() {
61
62 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
63
64 final NumericalPropagatorBuilder propagatorBuilder =
65 context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
66 1.0e-6, 60.0, 0.001);
67
68
69 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
70 propagatorBuilder);
71 final List<ObservedMeasurement<?>> measurements =
72 EstimationTestUtils.createMeasurements(propagator,
73 new AngularRaDecMeasurementCreator(context),
74 0.25, 3.0, 600.0);
75
76 propagator.clearStepHandlers();
77
78
79 final StreamingStatistics raDiffStat = new StreamingStatistics();
80 final StreamingStatistics decDiffStat = new StreamingStatistics();
81
82 for (final ObservedMeasurement<?> measurement : measurements) {
83
84
85 final AbsoluteDate datemeas = measurement.getDate();
86 SpacecraftState state = propagator.propagate(datemeas);
87
88
89 final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
90
91
92 raDiffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
93 decDiffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[1] - measurement.getObservedValue()[1]));
94 }
95
96
97 Assertions.assertEquals(0.0, raDiffStat.getMean(), 6.9e-11);
98 Assertions.assertEquals(0.0, raDiffStat.getStandardDeviation(), 8.5e-11);
99
100 Assertions.assertEquals(0.0, decDiffStat.getMean(), 4.5e-11);
101 Assertions.assertEquals(0.0, decDiffStat.getStandardDeviation(), 3e-11);
102 }
103
104
105
106
107 @Test
108 void testStateDerivatives() {
109
110 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
111
112 final NumericalPropagatorBuilder propagatorBuilder =
113 context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
114 1.0e-6, 60.0, 0.001);
115
116
117 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
118 propagatorBuilder);
119 final List<ObservedMeasurement<?>> measurements =
120 EstimationTestUtils.createMeasurements(propagator,
121 new AngularRaDecMeasurementCreator(context),
122 0.25, 3.0, 600.0);
123
124 propagator.clearStepHandlers();
125
126
127 double[] RaerrorsP = new double[3 * measurements.size()];
128 double[] RaerrorsV = new double[3 * measurements.size()];
129 double[] DecerrorsP = new double[3 * measurements.size()];
130 double[] DecerrorsV = new double[3 * measurements.size()];
131 int RaindexP = 0;
132 int RaindexV = 0;
133 int DecindexP = 0;
134 int DecindexV = 0;
135
136 for (final ObservedMeasurement<?> measurement : measurements) {
137
138
139 final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
140
141
142
143
144
145
146
147
148 final AbsoluteDate datemeas = measurement.getDate();
149 SpacecraftState state = propagator.propagate(datemeas);
150 final Vector3D stationP = stationParameter.getOffsetToInertial(state.getFrame(), datemeas, false).transformPosition(Vector3D.ZERO);
151 final double meanDelay = AbstractMeasurement.signalTimeOfFlightAdjustableEmitter(state.getPVCoordinates(), stationP,
152 datemeas, state.getFrame());
153
154 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
155 state = propagator.propagate(date);
156 final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
157 Assertions.assertEquals(2, estimated.getParticipants().length);
158 final double[][] jacobian = estimated.getStateDerivatives(0);
159
160
161 final double[][] finiteDifferencesJacobian =
162 Differentiation.differentiate(state1 -> measurement.
163 estimateWithoutDerivatives(new SpacecraftState[] {
164 state1
165 }).
166 getEstimatedValue(), measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN,
167 PositionAngleType.TRUE, 250.0, 4).value(state);
168
169 Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
170 Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
171
172 final double smallest = FastMath.ulp(1.0);
173
174 for (int i = 0; i < jacobian.length; ++i) {
175 for (int j = 0; j < jacobian[i].length; ++j) {
176 double relativeError = FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
177 finiteDifferencesJacobian[i][j]);
178
179 if ((FastMath.sqrt(finiteDifferencesJacobian[i][j]) < smallest) && (FastMath.sqrt(jacobian[i][j]) < smallest) ){
180 relativeError = 0.0;
181 }
182
183 if (j < 3) {
184 if (i == 0) {
185 RaerrorsP[RaindexP++] = relativeError;
186 } else {
187 DecerrorsP[DecindexP++] = relativeError;
188 }
189 } else {
190 if (i == 0) {
191 RaerrorsV[RaindexV++] = relativeError;
192 } else {
193 DecerrorsV[DecindexV++] = relativeError;
194 }
195 }
196 }
197 }
198 }
199
200
201 Assertions.assertEquals(0.0, new Median().evaluate(RaerrorsP), 4.8e-11);
202 Assertions.assertEquals(0.0, new Median().evaluate(RaerrorsV), 2.2e-5);
203
204
205 Assertions.assertEquals(0.0, new Median().evaluate(DecerrorsP), 2.2e-11);
206 Assertions.assertEquals(0.0, new Median().evaluate(DecerrorsV), 9.0e-6);
207
208
209 Assertions.assertEquals(AngularRaDec.MEASUREMENT_TYPE, measurements.get(0).getMeasurementType());
210 }
211
212
213
214
215 @Test
216 void testParameterDerivatives() {
217
218 Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
219
220 final NumericalPropagatorBuilder propagatorBuilder =
221 context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
222 1.0e-6, 60.0, 0.001);
223
224
225 for (final GroundStation station : context.stations) {
226 station.getClockOffsetDriver().setSelected(true);
227 station.getEastOffsetDriver().setSelected(true);
228 station.getNorthOffsetDriver().setSelected(true);
229 station.getZenithOffsetDriver().setSelected(true);
230 }
231 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
232 propagatorBuilder);
233 final List<ObservedMeasurement<?>> measurements =
234 EstimationTestUtils.createMeasurements(propagator,
235 new AngularRaDecMeasurementCreator(context),
236 0.25, 3.0, 600.0);
237 propagator.clearStepHandlers();
238
239 for (final ObservedMeasurement<?> measurement : measurements) {
240
241
242 final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
243
244
245
246
247
248
249
250
251 final AbsoluteDate datemeas = measurement.getDate();
252 final SpacecraftState stateini = propagator.propagate(datemeas);
253 final Vector3D stationP = stationParameter.getOffsetToInertial(stateini.getFrame(), datemeas, false).transformPosition(Vector3D.ZERO);
254 final double meanDelay = AbstractMeasurement.signalTimeOfFlightAdjustableEmitter(stateini.getPVCoordinates(), stationP,
255 datemeas, stateini.getFrame());
256
257 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
258 final SpacecraftState state = propagator.propagate(date);
259 final ParameterDriver[] drivers = new ParameterDriver[] {
260 stationParameter.getEastOffsetDriver(),
261 stationParameter.getNorthOffsetDriver(),
262 stationParameter.getZenithOffsetDriver()
263 };
264 for (int i = 0; i < 3; ++i) {
265 final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
266 Assertions.assertEquals(2, measurement.getDimension());
267 Assertions.assertEquals(2, gradient.length);
268
269 for (final int k : new int[] {0, 1}) {
270 final ParameterFunction dMkdP =
271 Differentiation.differentiate(new ParameterFunction() {
272
273 @Override
274 public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
275 return measurement.
276 estimateWithoutDerivatives(new SpacecraftState[] { state }).
277 getEstimatedValue()[k];
278 }
279 }, 3, 50.0 * drivers[i].getScale());
280 final double ref = dMkdP.value(drivers[i], date);
281
282 if (ref > 1.e-12) {
283 Assertions.assertEquals(ref, gradient[k], 3e-9 * FastMath.abs(ref));
284 }
285 }
286 }
287 }
288 }
289
290
291
292
293 @Test
294 void testIssue1026() {
295
296 Utils.setDataRoot("regular-data");
297
298 final double[] pos = { Constants.EGM96_EARTH_EQUATORIAL_RADIUS + 5e5, 1000., 0.};
299 final double[] vel = {0., 10., 0.};
300 final PVCoordinates pvCoordinates = new PVCoordinates(new Vector3D(pos[0], pos[1], pos[2]),
301 new Vector3D(vel[0], vel[1], vel[2]));
302 final AbsoluteDate epoch = new AbsoluteDate(new DateComponents(2000, 1, 1), TimeScalesFactory.getUTC());
303 final Frame gcrf = FramesFactory.getGCRF();
304 final CartesianOrbit orbit = new CartesianOrbit(pvCoordinates, gcrf, epoch, Constants.EGM96_EARTH_MU);
305 final SpacecraftState spacecraftState = new SpacecraftState(orbit);
306
307 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
308 Constants.IERS2010_EARTH_FLATTENING,
309 FramesFactory.getITRF(ITRFVersion.ITRF_2020, IERSConventions.IERS_2010, false));
310
311 final GeodeticPoint point = new GeodeticPoint(0., 0., 100.);
312 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "name");
313 final GroundStation station = new GroundStation(baseFrame);
314
315 final Frame[] frames = {FramesFactory.getEME2000(), FramesFactory.getGCRF(), FramesFactory.getICRF(), FramesFactory.getTOD(false)};
316 final double[][] raDec = new double[frames.length][];
317 for (int i = 0; i < frames.length; i++) {
318
319 final ObservableSatellite os = new ObservableSatellite(0);
320 final AngularRaDecBuilder builder = new AngularRaDecBuilder(null, station, frames[i],
321 new double[]{1., 1.}, new double[]{1., 1.}, os);
322 builder.init(spacecraftState.getDate(), spacecraftState.getDate());
323 final double[] moreRaDec = builder.build(spacecraftState.getDate(), new SpacecraftState[] { spacecraftState })
324 .getObservedValue();
325
326 final StaticTransform transform = frames[i].getStaticTransformTo(orbit.getFrame(), epoch);
327 final Vector3D transformedLoS = transform.transformVector(new Vector3D(moreRaDec[0], moreRaDec[1]));
328 raDec[i] = new double[] {FastMath.toDegrees(transformedLoS.getAlpha()),
329 FastMath.toDegrees(transformedLoS.getDelta())};
330 }
331
332 final double tolAngleDeg = 1e-2 / 3600.;
333 for (int i = 1; i < raDec.length; i++) {
334 Assertions.assertEquals(raDec[i][0], raDec[0][0], tolAngleDeg);
335 Assertions.assertEquals(raDec[i][1], raDec[0][1], tolAngleDeg);
336 }
337
338 }
339
340 }
341