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