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.junit.Assert;
24 import org.junit.Test;
25 import org.orekit.estimation.Context;
26 import org.orekit.estimation.EstimationTestUtils;
27 import org.orekit.estimation.measurements.AngularRaDec;
28 import org.orekit.estimation.measurements.AngularRaDecMeasurementCreator;
29 import org.orekit.estimation.measurements.ObservedMeasurement;
30 import org.orekit.estimation.measurements.PVMeasurementCreator;
31 import org.orekit.frames.Frame;
32 import org.orekit.orbits.KeplerianOrbit;
33 import org.orekit.orbits.OrbitType;
34 import org.orekit.orbits.PositionAngle;
35 import org.orekit.propagation.Propagator;
36 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
37 import org.orekit.time.AbsoluteDate;
38
39
40
41
42
43
44
45
46
47 public class IodGoodingTest {
48
49 @Test
50 public void testGooding() {
51 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
52
53 final double mu = context.initialOrbit.getMu();
54 final Frame frame = context.initialOrbit.getFrame();
55
56 final NumericalPropagatorBuilder propagatorBuilder =
57 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
58 1.0e-6, 60.0, 0.001);
59
60
61 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
62 propagatorBuilder);
63
64 final List<ObservedMeasurement<?>> measurements =
65 EstimationTestUtils.createMeasurements(propagator,
66 new PVMeasurementCreator(),
67 0.0, 1.0, 60.0);
68
69
70 final int idMeasure1 = 0;
71 final AbsoluteDate date1 = measurements.get(idMeasure1).getDate();
72 final Vector3D stapos1 = Vector3D.ZERO;
73
74
75
76 final Vector3D position1 = new Vector3D(measurements.get(idMeasure1).getObservedValue()[0],
77 measurements.get(idMeasure1).getObservedValue()[1],
78 measurements.get(idMeasure1).getObservedValue()[2]);
79 final double r1 = position1.getNorm();
80 final Vector3D lineOfSight1 = position1.normalize();
81
82
83 final int idMeasure2 = 20;
84 final AbsoluteDate date2 = measurements.get(idMeasure2).getDate();
85 final Vector3D stapos2 = Vector3D.ZERO;
86
87
88
89 final Vector3D position2 = new Vector3D(
90 measurements.get(idMeasure2).getObservedValue()[0],
91 measurements.get(idMeasure2).getObservedValue()[1],
92 measurements.get(idMeasure2).getObservedValue()[2]);
93 final Vector3D lineOfSight2 = position2.normalize();
94
95
96 final int idMeasure3 = 40;
97 final AbsoluteDate date3 = measurements.get(idMeasure3).getDate();
98 final Vector3D stapos3 = Vector3D.ZERO;
99
100
101
102 final Vector3D position3 = new Vector3D(
103 measurements.get(idMeasure3).getObservedValue()[0],
104 measurements.get(idMeasure3).getObservedValue()[1],
105 measurements.get(idMeasure3).getObservedValue()[2]);
106 final double r3 = position3.getNorm();
107 final Vector3D lineOfSight3 = position3.normalize();
108
109
110 final IodGooding iod = new IodGooding(mu);
111
112
113
114 final KeplerianOrbit orbit = iod.estimate(frame,
115 stapos1, stapos2, stapos3,
116 lineOfSight1, date1,
117 lineOfSight2, date2,
118 lineOfSight3, date3,
119 r1 * 1.0, r3 * 1.0);
120 Assert.assertEquals(orbit.getA(), context.initialOrbit.getA(), 1.0e-6 * context.initialOrbit.getA());
121 Assert.assertEquals(orbit.getE(), context.initialOrbit.getE(), 1.0e-6 * context.initialOrbit.getE());
122 Assert.assertEquals(orbit.getI(), context.initialOrbit.getI(), 1.0e-6 * context.initialOrbit.getI());
123
124 Assert.assertEquals(13127847.99808, iod.getRange1(), 1.0e-3);
125 Assert.assertEquals(13375711.51931, iod.getRange2(), 1.0e-3);
126 Assert.assertEquals(13950296.64852, iod.getRange3(), 1.0e-3);
127
128 }
129
130 @Test
131 public void testIssue756() {
132 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
133
134 final double mu = context.initialOrbit.getMu();
135 final Frame frame = context.initialOrbit.getFrame();
136
137 final NumericalPropagatorBuilder propagatorBuilder =
138 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
139 1.0e-6, 60.0, 0.001);
140
141
142 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
143 propagatorBuilder);
144
145 final List<ObservedMeasurement<?>> measurements =
146 EstimationTestUtils.createMeasurements(propagator,
147 new AngularRaDecMeasurementCreator(context),
148 0.0, 1.0, 60.0);
149
150
151 final AngularRaDec raDec1 = (AngularRaDec) measurements.get(0);
152 final AngularRaDec raDec2 = (AngularRaDec) measurements.get(20);
153 final AngularRaDec raDec3 = (AngularRaDec) measurements.get(40);
154
155
156 final double rhoInit1 = 1.3127847998082995E7;
157 final double rhoInit3 = 1.3950296648518201E7;
158
159
160 final IodGooding iod = new IodGooding(mu);
161
162
163
164 final KeplerianOrbit orbit1 = iod.estimate(frame, raDec1, raDec2, raDec3, rhoInit1, rhoInit3);
165 final KeplerianOrbit orbit2 = iod.estimate(frame,
166 stationPosition(frame, raDec1),
167 stationPosition(frame, raDec2),
168 stationPosition(frame, raDec3),
169 IodGooding.lineOfSight(raDec1), raDec1.getDate(),
170 IodGooding.lineOfSight(raDec2), raDec2.getDate(),
171 IodGooding.lineOfSight(raDec3), raDec3.getDate(),
172 rhoInit1, rhoInit3);
173
174 Assert.assertEquals(orbit1.getA(), orbit2.getA(), 1.0e-6 * orbit2.getA());
175 Assert.assertEquals(orbit1.getE(), orbit2.getE(), 1.0e-6 * orbit2.getE());
176 Assert.assertEquals(orbit1.getI(), orbit2.getI(), 1.0e-6 * orbit2.getI());
177 Assert.assertEquals(orbit1.getRightAscensionOfAscendingNode(), orbit2.getRightAscensionOfAscendingNode(), 1.0e-6 * orbit2.getRightAscensionOfAscendingNode());
178 Assert.assertEquals(orbit1.getPerigeeArgument(), orbit2.getPerigeeArgument(), 1.0e-6 * orbit2.getPerigeeArgument());
179 Assert.assertEquals(orbit1.getMeanAnomaly(), orbit2.getMeanAnomaly(), 1.0e-6 * orbit2.getMeanAnomaly());
180
181 }
182
183 private static Vector3D stationPosition(final Frame frame, final AngularRaDec raDec) {
184 return raDec.getStation().getBaseFrame().getPVCoordinates(raDec.getDate(), frame).getPosition();
185 }
186
187 }