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 org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.BeforeEach;
23 import org.junit.jupiter.api.Test;
24
25 import org.orekit.bodies.GeodeticPoint;
26 import org.orekit.bodies.OneAxisEllipsoid;
27 import org.orekit.estimation.measurements.AngularAzEl;
28 import org.orekit.estimation.measurements.AngularRaDec;
29 import org.orekit.estimation.measurements.GroundStation;
30 import org.orekit.estimation.measurements.ObservableSatellite;
31 import org.orekit.frames.TopocentricFrame;
32 import org.orekit.orbits.KeplerianOrbit;
33 import org.orekit.orbits.Orbit;
34 import org.orekit.orbits.PositionAngleType;
35 import org.orekit.propagation.Propagator;
36 import org.orekit.propagation.analytical.KeplerianPropagator;
37 import org.orekit.propagation.analytical.tle.TLE;
38 import org.orekit.propagation.analytical.tle.TLEPropagator;
39 import org.orekit.time.AbsoluteDate;
40 import org.orekit.time.TimeScalesFactory;
41 import org.orekit.utils.Constants;
42 import org.orekit.utils.TimeStampedPVCoordinates;
43
44
45
46
47
48 public class IodLaplaceTest extends AbstractIodTest {
49
50 @BeforeEach
51 public void observerOverride() {
52
53
54 final OneAxisEllipsoid body = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
55 Constants.WGS84_EARTH_FLATTENING, itrf);
56 this.observer = new GroundStation(
57 new TopocentricFrame(body, new GeodeticPoint(0.528253, -1.705768, 0.0), "Austin"));
58 this.observer.getPrimeMeridianOffsetDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
59 this.observer.getPolarOffsetXDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
60 this.observer.getPolarOffsetYDriver().setReferenceDate(AbsoluteDate.J2000_EPOCH);
61
62 }
63
64
65 @Test
66 public void testLaplaceKeplerian1() {
67 final AbsoluteDate date = new AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC());
68 final KeplerianOrbit kep = new KeplerianOrbit(6798938.970424857, 0.0021115522920270016, 0.9008866630545347,
69 1.8278985811406743, -2.7656136723308524,
70 0.8823034512437679, PositionAngleType.MEAN, gcrf,
71 date, Constants.EGM96_EARTH_MU);
72 final KeplerianPropagator prop = new KeplerianPropagator(kep);
73
74
75 final double[] error = estimateOrbit(prop, date, 30.0, 60.0).getErrorNorm();
76 Assertions.assertEquals(0.0, error[0], 275.0);
77 Assertions.assertEquals(0.0, error[1], 0.8);
78 }
79
80
81 @Test
82 public void testLaplaceKeplerian2() {
83 final AbsoluteDate date = new AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC());
84 final KeplerianOrbit kep = new KeplerianOrbit(42165414.60406032, 0.00021743441091199163, 0.0019139259842569903,
85 1.8142608912728584, 1.648821262690012,
86 0.11710513241172144, PositionAngleType.MEAN, gcrf,
87 date, Constants.EGM96_EARTH_MU);
88 final KeplerianPropagator prop = new KeplerianPropagator(kep);
89
90 final double[] error = estimateOrbit(prop, date, 60.0, 120.0).getErrorNorm();
91 Assertions.assertEquals(0.0, error[0], 395.0);
92 Assertions.assertEquals(0.0, error[1], 0.03);
93 }
94
95
96 @Test
97 public void testLaplaceTLE1() {
98 final String tle1 = "1 25544U 98067A 19271.53261574 .00000256 00000-0 12497-4 0 9993";
99 final String tle2 = "2 25544 51.6447 208.7465 0007429 92.6235 253.7389 15.50110361191281";
100
101 final TLE tleParser = new TLE(tle1, tle2);
102 final TLEPropagator tleProp = TLEPropagator.selectExtrapolator(tleParser);
103 final AbsoluteDate obsDate1 = tleParser.getDate();
104
105 final double[] error = estimateOrbit(tleProp, obsDate1, 30.0, 60.0).getErrorNorm();
106
107
108
109 Assertions.assertEquals(0.0, error[0], 5000.0);
110 Assertions.assertEquals(0.0, error[1], 10.0);
111 }
112
113
114 @Test
115 public void testLaplaceTLE2() {
116 final String tle1 = "1 4786U 70103A 19270.85687399 -.00000025 +00000-0 +00000-0 0 9998";
117 final String tle2 = "2 4786 055.8645 163.2517 1329144 016.0116 045.4806 08.42042146501862";
118
119 final TLE tleParser = new TLE(tle1, tle2);
120 final TLEPropagator tleProp = TLEPropagator.selectExtrapolator(tleParser);
121 final AbsoluteDate obsDate1 = tleParser.getDate();
122
123 final double[] error = estimateOrbit(tleProp, obsDate1, 30.0, 60.0).getErrorNorm();
124 Assertions.assertEquals(0.0, error[0], 5000.0);
125 Assertions.assertEquals(0.0, error[1], 10.0);
126 }
127
128
129 @Test
130 public void testLaplaceTLE3() {
131 final String tle1 = "1 28884U 05041A 19270.71989132 .00000058 00000-0 00000+0 0 9991";
132 final String tle2 = "2 28884 0.0023 190.3430 0001786 354.8402 307.2011 1.00272290 51019";
133
134 final TLE tleParser = new TLE(tle1, tle2);
135 final TLEPropagator tleProp = TLEPropagator.selectExtrapolator(tleParser);
136 final AbsoluteDate obsDate1 = tleParser.getDate();
137
138 final double[] error = estimateOrbit(tleProp, obsDate1, 300.0, 600.0).getErrorNorm();
139 Assertions.assertEquals(0.0, error[0], 5000.0);
140 Assertions.assertEquals(0.0, error[1], 10.0);
141 }
142
143 @Test
144 public void testIssue753() {
145
146
147 final AbsoluteDate date = new AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC());
148 final KeplerianOrbit kep = new KeplerianOrbit(6798938.970424857, 0.0021115522920270016, 0.9008866630545347,
149 1.8278985811406743, -2.7656136723308524,
150 0.8823034512437679, PositionAngleType.MEAN, gcrf,
151 date, Constants.EGM96_EARTH_MU);
152 final KeplerianPropagator prop = new KeplerianPropagator(kep);
153
154
155 final AngularRaDec raDec1 = new AngularRaDec(observer, gcrf, date,
156 new double[] {0.5380084720652177, -0.09320078788346774},
157 new double[] {1.0, 1.0}, new double[] {1.0, 1.0},
158 new ObservableSatellite(0));
159 final AngularRaDec raDec2 = new AngularRaDec(observer, gcrf, date.shiftedBy(30.0),
160 new double[] {0.549650227786601, -0.10753788558809535},
161 new double[] {1.0, 1.0}, new double[] {1.0, 1.0},
162 new ObservableSatellite(0));
163 final AngularRaDec raDec3 = new AngularRaDec(observer, gcrf, date.shiftedBy(60.0),
164 new double[] {0.5613500868283529, -0.12182129631017422},
165 new double[] {1.0, 1.0}, new double[] {1.0, 1.0},
166 new ObservableSatellite(0));
167
168
169 final IodLaplace laplace = new IodLaplace(Constants.EGM96_EARTH_MU);
170
171
172 final Orbit orbit = laplace.estimate(gcrf, raDec1, raDec2, raDec3);
173 final TimeStampedPVCoordinates ref = prop.getPVCoordinates(raDec2.getDate(), gcrf);
174
175
176 Assertions.assertEquals(0.0, ref.getPosition().distance(orbit.getPosition()), 275.0);
177 Assertions.assertEquals(0.0, ref.getVelocity().distance(orbit.getPVCoordinates().getVelocity()), 0.8);
178
179 }
180
181 @Test
182 public void testLaplaceKeplerianWithAzEl() {
183
184
185
186 final AbsoluteDate date = new AbsoluteDate(2019, 9, 29, 22, 0, 2.0, TimeScalesFactory.getUTC());
187 final KeplerianOrbit kep = new KeplerianOrbit(6798938.970424857, 0.0021115522920270016, 0.9008866630545347,
188 1.8278985811406743, -2.7656136723308524,
189 0.8823034512437679, PositionAngleType.MEAN, gcrf,
190 date, Constants.EGM96_EARTH_MU);
191 final KeplerianPropagator prop = new KeplerianPropagator(kep);
192
193
194 final AngularAzEl azEl1 = getAzEl(prop, date);
195 final AngularAzEl azEl2 = getAzEl(prop, date.shiftedBy(30.0));
196 final AngularAzEl azEl3 = getAzEl(prop, date.shiftedBy(60.0));
197
198
199 final Orbit estOrbit = new IodLaplace(Constants.EGM96_EARTH_MU).estimate(gcrf, azEl1, azEl2, azEl3);
200
201
202 final TimeStampedPVCoordinates truth = prop.getPVCoordinates(azEl2.getDate(), gcrf);
203 Assertions.assertEquals(0.0, Vector3D.distance(truth.getPosition(), estOrbit.getPosition()), 275.0);
204 Assertions.assertEquals(0.0, Vector3D.distance(truth.getVelocity(), estOrbit.getPVCoordinates().getVelocity()), 0.8);
205 }
206
207
208 public Result estimateOrbit(final Propagator prop, final AbsoluteDate obsDate1,
209 final double t2, final double t3) {
210
211
212 final Vector3D los1 = getLOSAngles(prop,obsDate1);
213
214 final AbsoluteDate obsDate2 = obsDate1.shiftedBy(t2);
215 final Vector3D los2 = getLOSAngles(prop,obsDate2);
216
217 final AbsoluteDate obsDate3 = obsDate1.shiftedBy(t3);
218 final Vector3D los3 = getLOSAngles(prop,obsDate3);
219
220 final TimeStampedPVCoordinates obsPva = observer
221 .getBaseFrame().getPVCoordinates(obsDate2, gcrf);
222
223
224 final TimeStampedPVCoordinates truth = prop.getPVCoordinates(obsDate2, gcrf);
225 final Orbit estOrbit = new IodLaplace(Constants.EGM96_EARTH_MU)
226 .estimate(gcrf, obsPva, obsDate1, los1, obsDate2, los2, obsDate3, los3);
227 return(new Result(truth, estOrbit));
228 }
229
230
231
232 private static class Result {
233 final private double[] errorNorm;
234
235 public Result(final TimeStampedPVCoordinates truth, final Orbit estOrbit)
236 {
237 this.errorNorm = new double[2];
238 this.errorNorm[0] = Vector3D.distance(truth.getPosition(),
239 estOrbit.getPosition());
240 this.errorNorm[1] = Vector3D.distance(truth.getVelocity(),
241 estOrbit.getPVCoordinates().getVelocity());
242 }
243
244 public double[] getErrorNorm()
245 {
246 return(this.errorNorm);
247 }
248 }
249 }