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.bodies.GeodeticPoint;
26 import org.orekit.bodies.OneAxisEllipsoid;
27 import org.orekit.estimation.measurements.AngularRaDec;
28 import org.orekit.estimation.measurements.GroundStation;
29 import org.orekit.estimation.measurements.ObservableSatellite;
30 import org.orekit.frames.Frame;
31 import org.orekit.frames.FramesFactory;
32 import org.orekit.frames.TopocentricFrame;
33 import org.orekit.orbits.CartesianOrbit;
34 import org.orekit.orbits.KeplerianOrbit;
35 import org.orekit.orbits.PositionAngle;
36 import org.orekit.propagation.Propagator;
37 import org.orekit.propagation.SpacecraftState;
38 import org.orekit.propagation.analytical.KeplerianPropagator;
39 import org.orekit.propagation.analytical.tle.TLE;
40 import org.orekit.propagation.analytical.tle.TLEPropagator;
41 import org.orekit.time.AbsoluteDate;
42 import org.orekit.time.TimeScalesFactory;
43 import org.orekit.utils.Constants;
44 import org.orekit.utils.IERSConventions;
45 import org.orekit.utils.AbsolutePVCoordinates;
46 import org.orekit.utils.TimeStampedPVCoordinates;
47 import org.orekit.Utils;
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
157 private Result estimateOrbit(final Propagator prop, final AbsoluteDate obsDate1,
158 final double t2, final double t3)
159 {
160
161 final Vector3D los1 = getLOSAngles(prop, obsDate1, observer);
162
163 final AbsoluteDate obsDate2 = obsDate1.shiftedBy(t2);
164 final Vector3D los2 = getLOSAngles(prop, obsDate2, observer);
165
166 final AbsoluteDate obsDate3 = obsDate1.shiftedBy(t3);
167 final Vector3D los3 = getLOSAngles(prop, obsDate3, observer);
168
169 final TimeStampedPVCoordinates obsPva = observer
170 .getBaseFrame().getPVCoordinates(obsDate2, gcrf);
171
172
173 final TimeStampedPVCoordinates truth = prop.getPVCoordinates(obsDate2, gcrf);
174 final CartesianOrbit estOrbit = new IodLaplace(Constants.EGM96_EARTH_MU)
175 .estimate(gcrf, obsPva, obsDate1, los1, obsDate2, los2, obsDate3, los3);
176 return(new Result(truth, estOrbit));
177 }
178
179
180
181 private Vector3D getLOSAngles(final Propagator prop, final AbsoluteDate date,
182 final GroundStation observer) {
183 final TimeStampedPVCoordinates satPvc = prop.getPVCoordinates(date, gcrf);
184 final AngularRaDec raDec = new AngularRaDec(observer, gcrf, date, new double[] {0.0, 0.0},
185 new double[] {1.0, 1.0},
186 new double[] {1.0, 1.0}, new ObservableSatellite(0));
187
188 final double[] angular = raDec.estimate(0, 0, new SpacecraftState[]{new SpacecraftState(
189 new AbsolutePVCoordinates(gcrf, satPvc))}).getEstimatedValue();
190 final double ra = angular[0];
191 final double dec = angular[1];
192
193 return(new Vector3D(FastMath.cos(dec)*FastMath.cos(ra),
194 FastMath.cos(dec)*FastMath.sin(ra), FastMath.sin(dec)));
195 }
196
197
198
199 private class Result
200 {
201 private double[] errorNorm;
202
203 public Result(final TimeStampedPVCoordinates truth, final CartesianOrbit estOrbit)
204 {
205 this.errorNorm = new double[2];
206 this.errorNorm[0] = Vector3D.distance(truth.getPosition(),
207 estOrbit.getPVCoordinates().getPosition());
208 this.errorNorm[1] = Vector3D.distance(truth.getVelocity(),
209 estOrbit.getPVCoordinates().getVelocity());
210 }
211
212 public double[] getErrorNorm()
213 {
214 return(this.errorNorm);
215 }
216 }
217 }