1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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.Utils;
26  import org.orekit.errors.OrekitException;
27  import org.orekit.errors.OrekitMessages;
28  import org.orekit.estimation.Context;
29  import org.orekit.estimation.EstimationTestUtils;
30  import org.orekit.estimation.measurements.ObservableSatellite;
31  import org.orekit.estimation.measurements.ObservedMeasurement;
32  import org.orekit.estimation.measurements.PV;
33  import org.orekit.estimation.measurements.PVMeasurementCreator;
34  import org.orekit.estimation.measurements.Position;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.orbits.KeplerianOrbit;
38  import org.orekit.orbits.OrbitType;
39  import org.orekit.orbits.PositionAngle;
40  import org.orekit.propagation.Propagator;
41  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
42  import org.orekit.time.AbsoluteDate;
43  import org.orekit.time.TimeScalesFactory;
44  import org.orekit.utils.Constants;
45  
46  /**
47   *
48   * Source: http://ccar.colorado.edu/asen5050/projects/projects_2012/kemble/gibbs_derivation.htm
49   *
50   * @author Joris Olympio
51   * @since 7.1
52   *
53   */
54  public class IodGibbsTest {
55  
56      @Test
57      public void testGibbs1() {
58          final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
59          final double mu = context.initialOrbit.getMu();
60          final Frame frame = context.initialOrbit.getFrame();
61  
62          final NumericalPropagatorBuilder propagatorBuilder =
63                          context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
64                                                1.0e-6, 60.0, 0.001);
65  
66          // create perfect range measurements
67          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
68                                                                             propagatorBuilder);
69          final ObservableSatellite satellite = new ObservableSatellite(0);
70  
71          final List<ObservedMeasurement<?>> measurements =
72                          EstimationTestUtils.createMeasurements(propagator,
73                                                                 new PVMeasurementCreator(),
74                                                                 0.0, 1.0, 60.0);
75  
76          final Vector3D position1 = new Vector3D(measurements.get(0).getObservedValue()[0],
77                                                  measurements.get(0).getObservedValue()[1],
78                                                  measurements.get(0).getObservedValue()[2]);
79          final PV pv1 = new PV(measurements.get(0).getDate(), position1, Vector3D.ZERO, 0., 0., 1., satellite);
80  
81          final Vector3D position2 = new Vector3D(measurements.get(1).getObservedValue()[0],
82                                                  measurements.get(1).getObservedValue()[1],
83                                                  measurements.get(1).getObservedValue()[2]);
84          final PV pv2 = new PV(measurements.get(1).getDate(), position2, Vector3D.ZERO, 0., 0., 1., satellite);
85  
86          final Vector3D position3 = new Vector3D(measurements.get(2).getObservedValue()[0],
87                                                  measurements.get(2).getObservedValue()[1],
88                                                  measurements.get(2).getObservedValue()[2]);
89          final PV pv3 = new PV(measurements.get(2).getDate(), position3, Vector3D.ZERO, 0., 0., 1., satellite);
90  
91          // instantiate the IOD method
92          final IodGibbs gibbs = new IodGibbs(mu);
93          final KeplerianOrbit orbit = gibbs.estimate(frame, pv1, pv2, pv3);
94  
95          Assert.assertEquals(context.initialOrbit.getA(), orbit.getA(), 1.0e-9 * context.initialOrbit.getA());
96          Assert.assertEquals(context.initialOrbit.getE(), orbit.getE(),  1.0e-9 * context.initialOrbit.getE());
97          Assert.assertEquals(context.initialOrbit.getI(),  orbit.getI(), 1.0e-9 * context.initialOrbit.getI());
98      }
99  
100     @Test
101     public void testGibbs2() {
102 
103         // test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-3, p457
104 
105         //extraction of the context.
106         final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
107         final double mu = context.initialOrbit.getMu();
108 
109         //initialisation
110         final IodGibbs gibbs = new IodGibbs(mu);
111 
112         // Observation  vector (EME2000)
113         final Vector3D posR1= new Vector3D(0.0, 0.0, 6378137.0);
114         final Vector3D posR2= new Vector3D(0.0, -4464696.0, -5102509.0);
115         final Vector3D posR3= new Vector3D(0.0, 5740323.0, 3189068);
116 
117         //epoch corresponding to the observation vector
118         AbsoluteDate dateRef = new AbsoluteDate(2000, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC());
119         AbsoluteDate date2 = dateRef.shiftedBy(76.48);
120         AbsoluteDate date3 = dateRef.shiftedBy(153.04);
121 
122         // Reference result (cf. Vallado)
123         final Vector3D velR2 = new Vector3D(0.0, 5531.148, -5191.806);
124 
125         //Gibbs IOD
126         final KeplerianOrbit orbit = gibbs.estimate(FramesFactory.getEME2000(),
127                                                     posR1, dateRef, posR2, date2, posR3, date3);
128 
129         //test
130         Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(), 1e-3);
131     }
132 
133     @Test
134     public void testGibbs3() {
135 
136         // test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-4, p463
137         // Remark: the test value in Vallado is performed with an Herrick-Gibbs methods but results are very close with Gibbs method.
138 
139         //extraction of context
140         final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
141         final double mu = context.initialOrbit.getMu();
142 
143         //Initialisation
144         final IodGibbs gibbs = new IodGibbs(mu);
145 
146         // Observations vector (EME2000)
147         final Vector3D posR1 = new Vector3D(3419855.64, 6019826.02, 2784600.22);
148         final Vector3D posR2 = new Vector3D(2935911.95, 6326183.24, 2660595.84);
149         final Vector3D posR3 = new Vector3D(2434952.02, 6597386.74, 2521523.11);
150 
151         //epoch corresponding to the observation vector
152         AbsoluteDate dateRef = new AbsoluteDate(2000, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC());
153         AbsoluteDate date2 = dateRef.shiftedBy(76.48);
154         AbsoluteDate date3 = dateRef.shiftedBy(153.04);
155 
156         // Reference result
157         final Vector3D velR2 = new Vector3D(-6441.632, 3777.625, -1720.582);
158 
159         //Gibbs IOD
160         final KeplerianOrbit orbit = gibbs.estimate(FramesFactory.getEME2000(),
161                                                     posR1, dateRef, posR2, date2, posR3, date3);
162 
163         //test for the norm of the velocity
164         Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(),  1e-3);
165 
166     }
167 
168     @Test
169     public void testIssue751() {
170 
171         // test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-4, p463
172         // Remark: the test value in Vallado is performed with an Herrick-Gibbs methods but results are very close with Gibbs method.
173 
174         Utils.setDataRoot("regular-data");
175 
176         // Initialisation
177         final IodGibbs gibbs = new IodGibbs(Constants.WGS84_EARTH_MU);
178 
179         // Observable satellite to initialize measurements
180         final ObservableSatellite satellite = new ObservableSatellite(0);
181 
182         // Observations vector (EME2000)
183         final Vector3D posR1 = new Vector3D(3419855.64, 6019826.02, 2784600.22);
184         final Vector3D posR2 = new Vector3D(2935911.95, 6326183.24, 2660595.84);
185         final Vector3D posR3 = new Vector3D(2434952.02, 6597386.74, 2521523.11);
186 
187         // Epoch corresponding to the observation vector
188         AbsoluteDate dateRef = new AbsoluteDate(2000, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC());
189         AbsoluteDate date2 = dateRef.shiftedBy(76.48);
190         AbsoluteDate date3 = dateRef.shiftedBy(153.04);
191 
192         // Reference result
193         final Vector3D velR2 = new Vector3D(-6441.632, 3777.625, -1720.582);
194 
195         // Gibbs IOD
196         final KeplerianOrbit orbit = gibbs.estimate(FramesFactory.getEME2000(),
197                                                     new Position(dateRef, posR1, 1.0, 1.0, satellite),
198                                                     new Position(date2,   posR2, 1.0, 1.0, satellite), 
199                                                     new Position(date3,   posR3, 1.0, 1.0, satellite));
200 
201         // Test for the norm of the velocity
202         Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(),  1e-3);
203 
204     }
205 
206     @Test
207     public void testNonDifferentDatesForObservations() {
208         Utils.setDataRoot("regular-data:potential:tides");
209 
210         // Initialization
211         final double mu = Constants.WGS84_EARTH_MU;
212         final Frame frame = FramesFactory.getEME2000();
213         final AbsoluteDate date = new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, TimeScalesFactory.getUTC());
214         final ObservableSatellite satellite = new ObservableSatellite(0);
215         final PV pv1 = new PV(date, Vector3D.PLUS_I, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
216         final PV pv2 = new PV(date, Vector3D.PLUS_J, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
217         final PV pv3 = new PV(date, Vector3D.PLUS_K, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
218 
219         // Create the IOD method
220         final IodGibbs gibbs = new IodGibbs(mu);
221 
222         try {
223             gibbs.estimate(frame, pv1, pv2, pv3);
224         } catch (OrekitException oe) {
225             Assert.assertEquals(OrekitMessages.NON_DIFFERENT_DATES_FOR_OBSERVATIONS, oe.getSpecifier());
226         }
227     }
228 
229     @Test
230     public void testNonCoplanarPoints() {
231         Utils.setDataRoot("regular-data:potential:tides");
232 
233         // Initialization
234         final double mu = Constants.WGS84_EARTH_MU;
235         final Frame frame = FramesFactory.getEME2000();
236         final AbsoluteDate date = new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, TimeScalesFactory.getUTC());
237         final ObservableSatellite satellite = new ObservableSatellite(0);
238         final PV pv1 = new PV(date, Vector3D.PLUS_I, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
239         final PV pv2 = new PV(date.shiftedBy(1.0), Vector3D.PLUS_J, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
240         final PV pv3 = new PV(date.shiftedBy(2.0), Vector3D.PLUS_K, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
241 
242         // Create the IOD method
243         final IodGibbs gibbs = new IodGibbs(mu);
244 
245         try {
246             gibbs.estimate(frame, pv1, pv2, pv3);
247         } catch (OrekitException oe) {
248             Assert.assertEquals(OrekitMessages.NON_COPLANAR_POINTS, oe.getSpecifier());
249         }
250     }
251 
252 }