1   /* Copyright 2022-2025 Bryan Cazabonne
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    * Bryan Cazabonne 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  package org.orekit.estimation.iod;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.Binary64Field;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeAll;
26  import org.junit.jupiter.api.Test;
27  import org.orekit.Utils;
28  import org.orekit.errors.OrekitException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.estimation.measurements.ObservableSatellite;
31  import org.orekit.estimation.measurements.PV;
32  import org.orekit.estimation.measurements.Position;
33  import org.orekit.frames.FramesFactory;
34  import org.orekit.orbits.FieldOrbit;
35  import org.orekit.orbits.Orbit;
36  import org.orekit.time.AbsoluteDate;
37  import org.orekit.time.FieldAbsoluteDate;
38  import org.orekit.time.TimeScalesFactory;
39  import org.orekit.utils.Constants;
40  
41  class IodHerrickGibbsTest {
42  
43      @Test
44      public void testCompareWithVallado() {
45  
46          // Test extracted from Vallado D., Fundamentals of astrodynamics & applications, 4rd Edition, Example 7-4, p467.
47  
48          // Initialisation
49          final IodHerrickGibbs gibbs = new IodHerrickGibbs(398600.4418);
50  
51          // Observable satellite to initialize measurements
52          final ObservableSatellite satellite = new ObservableSatellite(0);
53  
54          // Observations vector (EME2000)
55          final Vector3D r1 = new Vector3D(3419.85564, 6019.82602, 2784.60022);
56          final Vector3D r2 = new Vector3D(2935.91195, 6326.18324, 2660.59584);
57          final Vector3D r3 = new Vector3D(2434.95202, 6597.38674, 2521.52311);
58  
59          // Epoch corresponding to the observation vector
60          AbsoluteDate t1 = AbsoluteDate.J2000_EPOCH;
61          AbsoluteDate t2 = t1.shiftedBy(76.48);
62          AbsoluteDate t3 = t1.shiftedBy(153.04);
63  
64          // Reference result (from Vallado's example)
65          //
66          // IMPORTANT: Reference results in the 4th Edition of the book are not correct! The reference below comes from
67          // Vallado's validation test available in its Github repository (see test_hgibbs()):
68          // https://github.com/CelesTrak/fundamentals-of-astrodynamics/blob/main/software/python/tests/astro/iod/test_utils.py
69          final Vector3D referenceV2 = new Vector3D(-6.441557227511062, 3.777559606719521, -1.7205675602414345);
70  
71          // Herrick-Gibbs IOD
72          final Orbit orbit = gibbs.estimate(FramesFactory.getEME2000(),
73                                             new Position(t1, r1, 1.0, 1.0, satellite),
74                                             new Position(t2, r2, 1.0, 1.0, satellite),
75                                             new Position(t3, r3, 1.0, 1.0, satellite));
76  
77          // Verify
78          Assertions.assertEquals(0.0, orbit.durationFrom(t2));
79          Assertions.assertEquals(r2.getNorm(), orbit.getPosition().getNorm(), 1.0e-10);
80          Assertions.assertEquals(referenceV2.getNorm(), orbit.getVelocity().getNorm(),  1.0e-10);
81  
82      }
83  
84      @Test
85      public void testCompareWithValladoField() {
86          doTestCompareWithValladoField(Binary64Field.getInstance());
87      }
88  
89      private <T extends CalculusFieldElement<T>> void doTestCompareWithValladoField(Field<T> field) {
90  
91          // Test extracted from Vallado D., Fundamentals of astrodynamics & applications, 4rd Edition, Example 7-4, p467.
92  
93          // Initialisation
94          final IodHerrickGibbs gibbs = new IodHerrickGibbs(398600.4418);
95  
96          // Observations vector (EME2000)
97          final FieldVector3D<T> r1 = new FieldVector3D<>(field, new Vector3D(3419.85564, 6019.82602, 2784.60022));
98          final FieldVector3D<T> r2 = new FieldVector3D<>(field, new Vector3D(2935.91195, 6326.18324, 2660.59584));
99          final FieldVector3D<T> r3 = new FieldVector3D<>(field, new Vector3D(2434.95202, 6597.38674, 2521.52311));
100 
101         // Epoch corresponding to the observation vector
102         FieldAbsoluteDate<T> t1 = FieldAbsoluteDate.getJ2000Epoch(field);
103         FieldAbsoluteDate<T> t2 = t1.shiftedBy(76.48);
104         FieldAbsoluteDate<T> t3 = t1.shiftedBy(153.04);
105 
106         // Reference result (from Vallado's example)
107         //
108         // IMPORTANT: Reference results in the 4th Edition of the book are not correct! The reference below comes from
109         // Vallado's validation test available in its Github repository (see test_hgibbs()):
110         // https://github.com/CelesTrak/fundamentals-of-astrodynamics/blob/main/software/python/tests/astro/iod/test_utils.py
111         final Vector3D referenceV2 = new Vector3D(-6.441557227511062, 3.777559606719521, -1.7205675602414345);
112 
113         // Herrick-Gibbs IOD
114         final FieldOrbit<T> orbit = gibbs.estimate(FramesFactory.getEME2000(), r1, t1, r2, t2, r3, t3);
115 
116         // Verify
117         Assertions.assertEquals(0.0, orbit.durationFrom(t2).getReal());
118         Assertions.assertEquals(r2.getNorm().getReal(), orbit.getPosition().getNorm().getReal(), 1.0e-10);
119         Assertions.assertEquals(referenceV2.getNorm(), orbit.getVelocity().getNorm().getReal(),  1.0e-10);
120 
121     }
122 
123     @Test
124     public void testNonDifferentDatesForObservations() {
125         // Initialization
126         final AbsoluteDate date = new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, TimeScalesFactory.getUTC());
127         final ObservableSatellite satellite = new ObservableSatellite(0);
128         final PV pv1 = new PV(date, Vector3D.PLUS_I, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
129         final PV pv2 = new PV(date.shiftedBy(1.0), Vector3D.PLUS_J, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
130         final PV pv3 = new PV(date, Vector3D.PLUS_K, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
131         // Action
132         try {
133             new IodHerrickGibbs(Constants.WGS84_EARTH_MU).estimate(FramesFactory.getEME2000(), pv1, pv2, pv3);
134         } catch (OrekitException oe) {
135             Assertions.assertEquals(OrekitMessages.NON_DIFFERENT_DATES_FOR_OBSERVATIONS, oe.getSpecifier());
136         }
137     }
138 
139     @Test
140     public void testNonDifferentDatesForObservationsField() {
141         doTestNonDifferentDatesForObservationsField(Binary64Field.getInstance());
142     }
143 
144     private <T extends CalculusFieldElement<T>> void doTestNonDifferentDatesForObservationsField(Field<T> field) {
145         // Initialization
146         final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
147         // Action
148         try {
149             new IodHerrickGibbs(Constants.WGS84_EARTH_MU).estimate(FramesFactory.getEME2000(), FieldVector3D.getPlusI(field), date,
150                                                                    FieldVector3D.getPlusJ(field), date,
151                                                                    FieldVector3D.getPlusK(field), date);
152         } catch (OrekitException oe) {
153             Assertions.assertEquals(OrekitMessages.NON_DIFFERENT_DATES_FOR_OBSERVATIONS, oe.getSpecifier());
154         }
155     }
156 
157     @Test
158     public void testNonCoplanarPoints() {
159         // Initialization
160         final AbsoluteDate date = new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, TimeScalesFactory.getUTC());
161         final ObservableSatellite satellite = new ObservableSatellite(0);
162         final PV pv1 = new PV(date, Vector3D.PLUS_I, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
163         final PV pv2 = new PV(date.shiftedBy(1.0), Vector3D.PLUS_J, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
164         final PV pv3 = new PV(date.shiftedBy(2.0), Vector3D.PLUS_K, Vector3D.ZERO, 1.0, 1.0, 1.0, satellite);
165         // Action
166         try {
167             new IodHerrickGibbs(Constants.WGS84_EARTH_MU).estimate(FramesFactory.getEME2000(), pv1, pv2, pv3);
168         } catch (OrekitException oe) {
169             Assertions.assertEquals(OrekitMessages.NON_COPLANAR_POINTS, oe.getSpecifier());
170         }
171     }
172 
173     @Test
174     public void testNonCoplanarPointsField() {
175         doTestNonCoplanarPointsField(Binary64Field.getInstance());
176     }
177 
178     private  <T extends CalculusFieldElement<T>> void doTestNonCoplanarPointsField(Field<T> field) {
179         // Initialization
180         final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
181         // Action
182         try {
183             new IodHerrickGibbs(Constants.WGS84_EARTH_MU).estimate(FramesFactory.getEME2000(), FieldVector3D.getPlusI(field), date,
184                                                                    FieldVector3D.getPlusJ(field), date.shiftedBy(1.0),
185                                                                    FieldVector3D.getPlusK(field), date.shiftedBy(2.0));
186         } catch (OrekitException oe) {
187             Assertions.assertEquals(OrekitMessages.NON_COPLANAR_POINTS, oe.getSpecifier());
188         }
189     }
190 
191     @BeforeAll
192     public static void init() {
193         Utils.setDataRoot("regular-data:potential:tides");
194     }
195 }