1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.orekit.estimation.iod;
19
20 import java.util.List;
21
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.util.FastMath;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.Test;
26 import org.orekit.bodies.GeodeticPoint;
27 import org.orekit.bodies.OneAxisEllipsoid;
28 import org.orekit.estimation.Context;
29 import org.orekit.estimation.EstimationTestUtils;
30 import org.orekit.estimation.measurements.*;
31 import org.orekit.frames.Frame;
32 import org.orekit.frames.FramesFactory;
33 import org.orekit.frames.TopocentricFrame;
34 import org.orekit.orbits.KeplerianOrbit;
35 import org.orekit.orbits.Orbit;
36 import org.orekit.orbits.OrbitType;
37 import org.orekit.orbits.PositionAngleType;
38 import org.orekit.propagation.Propagator;
39 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
40 import org.orekit.time.AbsoluteDate;
41 import org.orekit.time.Month;
42 import org.orekit.time.TimeScalesFactory;
43 import org.orekit.utils.Constants;
44 import org.orekit.utils.IERSConventions;
45
46
47
48
49
50
51
52
53
54 public class IodGoodingTest extends AbstractIodTest {
55
56
57
58 @Test
59 public void testIssue1166RaDec() {
60 AbsoluteDate t1 = new AbsoluteDate(2023, Month.JUNE, 9, 17, 4,59.10, TimeScalesFactory.getUTC());
61 AbsoluteDate t2 = new AbsoluteDate(2023, Month.JUNE, 9, 17, 10,50.66, TimeScalesFactory.getUTC());
62 AbsoluteDate t3 = new AbsoluteDate(2023, Month.JUNE, 9, 17, 16,21.09, TimeScalesFactory.getUTC());
63
64 Vector3D RA = new Vector3D((15.* (16. + 5./60. + 51.20/3600.)),
65 (15.*(16.+ 11./60. + 43.73/3600.)), (15.*(16.+ 17./60. + 15.1/3600. )));
66
67 Vector3D DEC = new Vector3D((-(6.+ 31./60. + 44.22/3600.)),
68 (-(6. + 31./60. + 52.36/3600.)),
69 (-(6. +32./60. + 0.03/3600.)));
70
71 Frame ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
72 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS ,Constants.WGS84_EARTH_FLATTENING, ITRF);
73 GeodeticPoint stationCoord = new GeodeticPoint(FastMath.toRadians(43.05722), FastMath.toRadians(76.971667), 2735.0);
74 TopocentricFrame stationFrame = new TopocentricFrame(earth, stationCoord, "N42");
75 GroundStation ground_station = new GroundStation(stationFrame);
76
77 double[] angular1 = {FastMath.toRadians(RA.getX()), FastMath.toRadians(DEC.getX())};
78 double[] angular2 = {FastMath.toRadians(RA.getY()), FastMath.toRadians(DEC.getY())};
79 double[] angular3 = {FastMath.toRadians(RA.getZ()), FastMath.toRadians(DEC.getZ())};
80
81 AngularRaDec raDec1 = new AngularRaDec(ground_station, FramesFactory.getEME2000(), t1,
82 angular1, new double[] {1.0, 1.0}, new double[] {1.0, 1.0}, new ObservableSatellite(1));
83 AngularRaDec raDec2 = new AngularRaDec(ground_station, FramesFactory.getEME2000(), t2,
84 angular2, new double[] {1.0, 1.0}, new double[] {1.0, 1.0}, new ObservableSatellite(1));
85 AngularRaDec raDec3 = new AngularRaDec(ground_station, FramesFactory.getEME2000(), t3,
86 angular3, new double[] {1.0, 1.0}, new double[] {1.0, 1.0}, new ObservableSatellite(1));
87
88
89
90
91
92 Orbit estimated_orbit_Gooding = new IodGooding(mu).estimate(eme2000, raDec1,raDec2,raDec3);
93 KeplerianOrbit orbitGooding = new KeplerianOrbit(estimated_orbit_Gooding);
94 Assertions.assertEquals(4.2394187540973224E7, orbitGooding.getA(), 1.0e-10);
95 Assertions.assertEquals(0.004411368860770379, orbitGooding.getE(), 1.0e-10);
96 Assertions.assertEquals(FastMath.toRadians(0.09185983299662298), orbitGooding.getI(), 1.0e-10);
97 Assertions.assertEquals(FastMath.toRadians(169.74389246605776), orbitGooding.getPerigeeArgument(), 1.0e-10);
98 Assertions.assertEquals(FastMath.toRadians(90.92874061328043), orbitGooding.getRightAscensionOfAscendingNode(), 1.0e-10);
99 Assertions.assertEquals(FastMath.toRadians(-18.909215663128727), orbitGooding.getTrueAnomaly(), 1.0e-10);
100 }
101
102 @Test
103 public void testIssue1166AzEl() {
104
105 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
106 final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true, 1.0e-6, 60.0, 0.001);
107 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
108 final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator,
109 new AngularAzElMeasurementCreator(context),
110 0.0, 1.0, 60.0);
111 final AngularAzEl azEl1 = (AngularAzEl) measurements.get(0);
112 final AngularAzEl azEl2 = (AngularAzEl) measurements.get(20);
113 final AngularAzEl azEl3 = (AngularAzEl) measurements.get(40);
114
115
116
117 Orbit estimated_orbit_Gooding = new IodGooding(mu).estimate(eme2000, azEl1,azEl2,azEl3);
118 KeplerianOrbit orbitGooding = new KeplerianOrbit(estimated_orbit_Gooding);
119 Assertions.assertEquals(1.4197961507698055E7, orbitGooding.getA(), 1.0e-10);
120 Assertions.assertEquals(0.16923654961240223, orbitGooding.getE(), 1.0e-10);
121 Assertions.assertEquals(FastMath.toRadians(71.52638181160407), orbitGooding.getI(), 1.0e-10);
122 Assertions.assertEquals(FastMath.toRadians(21.450082668672675), orbitGooding.getPerigeeArgument(), 1.0e-10);
123 Assertions.assertEquals(FastMath.toRadians(78.76324220205018), orbitGooding.getRightAscensionOfAscendingNode(), 1.0e-10);
124 Assertions.assertEquals(FastMath.toRadians(-163.62886990452034), orbitGooding.getTrueAnomaly(), 1.0e-10);
125 }
126
127 @Test
128 public void testGooding() {
129 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
130
131 final double mu = context.initialOrbit.getMu();
132 final Frame frame = context.initialOrbit.getFrame();
133
134 final NumericalPropagatorBuilder propagatorBuilder =
135 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
136 1.0e-6, 60.0, 0.001);
137
138
139 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
140 propagatorBuilder);
141
142 final List<ObservedMeasurement<?>> measurements =
143 EstimationTestUtils.createMeasurements(propagator,
144 new PVMeasurementCreator(),
145 0.0, 1.0, 60.0);
146
147
148 final int idMeasure1 = 0;
149 final AbsoluteDate date1 = measurements.get(idMeasure1).getDate();
150 final Vector3D stapos1 = Vector3D.ZERO;
151
152
153
154 final Vector3D position1 = new Vector3D(measurements.get(idMeasure1).getObservedValue()[0],
155 measurements.get(idMeasure1).getObservedValue()[1],
156 measurements.get(idMeasure1).getObservedValue()[2]);
157 final double r1 = position1.getNorm();
158 final Vector3D lineOfSight1 = position1.normalize();
159
160
161 final int idMeasure2 = 20;
162 final AbsoluteDate date2 = measurements.get(idMeasure2).getDate();
163 final Vector3D stapos2 = Vector3D.ZERO;
164
165
166
167 final Vector3D position2 = new Vector3D(
168 measurements.get(idMeasure2).getObservedValue()[0],
169 measurements.get(idMeasure2).getObservedValue()[1],
170 measurements.get(idMeasure2).getObservedValue()[2]);
171 final Vector3D lineOfSight2 = position2.normalize();
172
173
174 final int idMeasure3 = 40;
175 final AbsoluteDate date3 = measurements.get(idMeasure3).getDate();
176 final Vector3D stapos3 = Vector3D.ZERO;
177
178
179
180 final Vector3D position3 = new Vector3D(
181 measurements.get(idMeasure3).getObservedValue()[0],
182 measurements.get(idMeasure3).getObservedValue()[1],
183 measurements.get(idMeasure3).getObservedValue()[2]);
184 final double r3 = position3.getNorm();
185 final Vector3D lineOfSight3 = position3.normalize();
186
187
188 final IodGooding iod = new IodGooding(mu);
189
190
191
192 final Orbit orbit = iod.estimate(frame,
193 stapos1, stapos2, stapos3,
194 lineOfSight1, date1,
195 lineOfSight2, date2,
196 lineOfSight3, date3,
197 r1 * 1.0, r3 * 1.0);
198 Assertions.assertEquals(orbit.getA(), context.initialOrbit.getA(), 1.0e-6 * context.initialOrbit.getA());
199 Assertions.assertEquals(orbit.getE(), context.initialOrbit.getE(), 1.0e-6 * context.initialOrbit.getE());
200 Assertions.assertEquals(orbit.getI(), context.initialOrbit.getI(), 1.0e-6 * context.initialOrbit.getI());
201
202 Assertions.assertEquals(13127847.99808, iod.getRange1(), 1.0e-3);
203 Assertions.assertEquals(13375711.51931, iod.getRange2(), 1.0e-3);
204 Assertions.assertEquals(13950296.64852, iod.getRange3(), 1.0e-3);
205
206
207 }
208
209 @Test
210 public void testIssue756() {
211 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
212
213 final double mu = context.initialOrbit.getMu();
214 final Frame frame = context.initialOrbit.getFrame();
215
216 final NumericalPropagatorBuilder propagatorBuilder =
217 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
218 1.0e-6, 60.0, 0.001);
219
220
221 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
222 propagatorBuilder);
223
224 final List<ObservedMeasurement<?>> measurements =
225 EstimationTestUtils.createMeasurements(propagator,
226 new AngularRaDecMeasurementCreator(context),
227 0.0, 1.0, 60.0);
228
229
230 final AngularRaDec raDec1 = (AngularRaDec) measurements.get(0);
231 final AngularRaDec raDec2 = (AngularRaDec) measurements.get(20);
232 final AngularRaDec raDec3 = (AngularRaDec) measurements.get(40);
233
234
235 final double rhoInit1 = 1.3127847998082995E7;
236 final double rhoInit3 = 1.3950296648518201E7;
237
238
239 final IodGooding iod = new IodGooding(mu);
240
241 final KeplerianOrbit orbit1 = new KeplerianOrbit(iod.estimate(frame, raDec1, raDec2, raDec3, rhoInit1, rhoInit3));
242 final KeplerianOrbit orbit2 = new KeplerianOrbit(iod.estimate(frame,
243 raDec1.getGroundStationPosition(frame),
244 raDec2.getGroundStationPosition(frame),
245 raDec3.getGroundStationPosition(frame),
246 raDec1.getObservedLineOfSight(frame), raDec1.getDate(),
247 raDec2.getObservedLineOfSight(frame), raDec2.getDate(),
248 raDec3.getObservedLineOfSight(frame), raDec3.getDate(),
249 rhoInit1, rhoInit3));
250
251 Assertions.assertEquals(orbit1.getA(), orbit2.getA(), 1.0e-6 * orbit2.getA());
252 Assertions.assertEquals(orbit1.getE(), orbit2.getE(), 1.0e-6 * orbit2.getE());
253 Assertions.assertEquals(orbit1.getI(), orbit2.getI(), 1.0e-6 * orbit2.getI());
254 Assertions.assertEquals(orbit1.getRightAscensionOfAscendingNode(), orbit2.getRightAscensionOfAscendingNode(), 1.0e-6 * orbit2.getRightAscensionOfAscendingNode());
255 Assertions.assertEquals(orbit1.getPerigeeArgument(), orbit2.getPerigeeArgument(), FastMath.abs(1.0e-6 * orbit2.getPerigeeArgument()));
256 Assertions.assertEquals(orbit1.getMeanAnomaly(), orbit2.getMeanAnomaly(), 1.0e-6 * orbit2.getMeanAnomaly());
257 Assertions.assertEquals(orbit1.getMeanAnomaly(), orbit2.getMeanAnomaly(), 1.0e-6 * orbit2.getMeanAnomaly());
258 }
259
260 @Test
261 public void testIssue1216() {
262 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
263
264 final double mu = context.initialOrbit.getMu();
265 final Frame frame = context.initialOrbit.getFrame();
266
267 final NumericalPropagatorBuilder propagatorBuilder =
268 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
269 1.0e-6, 60.0, 0.001);
270
271
272 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
273 propagatorBuilder);
274
275 final List<ObservedMeasurement<?>> measurements =
276 EstimationTestUtils.createMeasurements(propagator,
277 new AngularAzElMeasurementCreator(context),
278 0.0, 1.0, 60.0);
279
280
281 final AngularAzEl azEl1 = (AngularAzEl) measurements.get(0);
282 final AngularAzEl azEl2 = (AngularAzEl) measurements.get(20);
283 final AngularAzEl azEl3 = (AngularAzEl) measurements.get(40);
284
285
286 final double rhoInit1 = 1.3127847998082995E7;
287 final double rhoInit3 = 1.3950296648518201E7;
288
289
290 final IodGooding iod = new IodGooding(mu);
291
292 final KeplerianOrbit orbit1 = new KeplerianOrbit(iod.estimate(frame, azEl1, azEl2, azEl3, rhoInit1, rhoInit3));
293 final KeplerianOrbit orbit2 = new KeplerianOrbit(iod.estimate(frame,
294 azEl1.getGroundStationPosition(frame),
295 azEl2.getGroundStationPosition(frame),
296 azEl3.getGroundStationPosition(frame),
297 azEl1.getObservedLineOfSight(frame), azEl1.getDate(),
298 azEl2.getObservedLineOfSight(frame), azEl2.getDate(),
299 azEl3.getObservedLineOfSight(frame), azEl3.getDate(),
300 rhoInit1, rhoInit3));
301
302 Assertions.assertEquals(orbit1.getA(), orbit2.getA(), 1.0e-6 * orbit2.getA());
303 Assertions.assertEquals(orbit1.getE(), orbit2.getE(), 1.0e-6 * orbit2.getE());
304 Assertions.assertEquals(orbit1.getI(), orbit2.getI(), 1.0e-6 * orbit2.getI());
305 Assertions.assertEquals(orbit1.getRightAscensionOfAscendingNode(), orbit2.getRightAscensionOfAscendingNode(), 1.0e-6 * orbit2.getRightAscensionOfAscendingNode());
306 Assertions.assertEquals(orbit1.getPerigeeArgument(), orbit2.getPerigeeArgument(), FastMath.abs(1.0e-6 * orbit2.getPerigeeArgument()));
307 Assertions.assertEquals(orbit1.getMeanAnomaly(), orbit2.getMeanAnomaly(), 1.0e-6 * orbit2.getMeanAnomaly());
308 Assertions.assertEquals(orbit1.getMeanAnomaly(), orbit2.getMeanAnomaly(), 1.0e-6 * orbit2.getMeanAnomaly());
309 }
310
311 }