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  package org.orekit.propagation.analytical.gnss;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.FastMath;
24  import org.junit.Assert;
25  import org.junit.Before;
26  import org.junit.BeforeClass;
27  import org.junit.Test;
28  import org.orekit.Utils;
29  import org.orekit.errors.OrekitException;
30  import org.orekit.errors.OrekitMessages;
31  import org.orekit.frames.Frame;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.gnss.SatelliteSystem;
34  import org.orekit.propagation.analytical.gnss.data.GNSSOrbitalElements;
35  import org.orekit.propagation.analytical.gnss.data.GalileoAlmanac;
36  import org.orekit.propagation.analytical.gnss.data.GalileoNavigationMessage;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.GNSSDate;
39  import org.orekit.time.TimeScalesFactory;
40  import org.orekit.utils.CartesianDerivativesFilter;
41  import org.orekit.utils.Constants;
42  import org.orekit.utils.IERSConventions;
43  import org.orekit.utils.PVCoordinates;
44  import org.orekit.utils.TimeStampedPVCoordinates;
45  
46  public class GalileoPropagatorTest {
47  
48      private GalileoNavigationMessage goe;
49  
50      @Before
51      public void setUp() {
52          goe = new GalileoNavigationMessage();
53          goe.setPRN(4);
54          goe.setWeek(1024);
55          goe.setTime(293400.0);
56          goe.setSqrtA(5440.602949142456);
57          goe.setDeltaN(3.7394414770330066E-9);
58          goe.setE(2.4088891223073006E-4);
59          goe.setI0(0.9531656087278083);
60          goe.setIDot(-2.36081262303612E-10);
61          goe.setOmega0(-0.36639513583951266);
62          goe.setOmegaDot(-5.7695260382035525E-9);
63          goe.setPa(-1.6870064194345724);
64          goe.setM0(-0.38716557650888);
65          goe.setCuc(-8.903443813323975E-7);
66          goe.setCus(6.61797821521759E-6);
67          goe.setCrc(194.0625);
68          goe.setCrs(-18.78125);
69          goe.setCic(3.166496753692627E-8);
70          goe.setCis(-1.862645149230957E-8);
71          goe.setDate(new GNSSDate(goe.getWeek(), 1000. * goe.getTime(), SatelliteSystem.GALILEO).getDate());
72      }
73      
74      @BeforeClass
75      public static void setUpBeforeClass() {
76          Utils.setDataRoot("gnss");
77      }
78  
79      @Test
80      public void testGalileoCycle() {
81          // Reference for the almanac: 2019-05-28T09:40:01.0Z
82          final GalileoAlmanac almanac = new GalileoAlmanac();
83          almanac.setPRN(1);
84          almanac.setWeek(1024);
85          almanac.setTime(293400.0);
86          almanac.setDeltaSqrtA(0.013671875);
87          almanac.setE(0.000152587890625);
88          almanac.setDeltaInc(0.003356933593);
89          almanac.setIOD(4);
90          almanac.setOmega0(0.2739257812499857891);
91          almanac.setOmegaDot(-1.74622982740407E-9);
92          almanac.setPa(0.7363586425);
93          almanac.setM0(0.27276611328124);
94          almanac.setAf0(-0.0006141662597);
95          almanac.setAf1(-7.275957614183E-12);
96          almanac.setHealthE1(0);
97          almanac.setHealthE5a(0);
98          almanac.setHealthE5b(0);
99          almanac.setDate(new GNSSDate(almanac.getWeek(), 1000.0 * almanac.getTime(), SatelliteSystem.GALILEO).getDate());
100 
101         // Intermediate verification
102         Assert.assertEquals(1,                   almanac.getPRN());
103         Assert.assertEquals(1024,                almanac.getWeek());
104         Assert.assertEquals(4,                   almanac.getIOD());
105         Assert.assertEquals(0,                   almanac.getHealthE1());
106         Assert.assertEquals(0,                   almanac.getHealthE5a());
107         Assert.assertEquals(0,                   almanac.getHealthE5b());
108         Assert.assertEquals(-0.0006141662597,    almanac.getAf0(), 1.0e-15);
109         Assert.assertEquals(-7.275957614183E-12, almanac.getAf1(), 1.0e-15);
110 
111         // Builds the GalileoPropagator from the almanac
112         GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
113         // Propagate at the Galileo date and one Galileo cycle later
114         final AbsoluteDate date0 = almanac.getDate();
115         final Vector3D p0 = propagator.propagateInEcef(date0).getPosition();
116         final double galCycleDuration = almanac.getCycleDuration();
117         final AbsoluteDate date1 = date0.shiftedBy(galCycleDuration);
118         final Vector3D p1 = propagator.propagateInEcef(date1).getPosition();
119 
120         // Checks
121         Assert.assertEquals(0., p0.distance(p1), 0.);
122     }
123 
124     @Test
125     public void testFrames() {
126         // Builds the GalileoPropagator from the ephemeris
127         GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
128         Assert.assertEquals("EME2000", propagator.getFrame().getName());
129         Assert.assertEquals(3.986004418e+14, goe.getMu(), 1.0e6);
130         // Defines some date
131         final AbsoluteDate date = new AbsoluteDate(2016, 3, 3, 12, 0, 0., TimeScalesFactory.getUTC());
132         // Get PVCoordinates at the date in the ECEF
133         final PVCoordinates pv0 = propagator.propagateInEcef(date);
134         // Get PVCoordinates at the date in the ECEF
135         final PVCoordinates pv1 = propagator.getPVCoordinates(date, propagator.getECEF());
136 
137         // Checks
138         Assert.assertEquals(0., pv0.getPosition().distance(pv1.getPosition()), 2.4e-8);
139         Assert.assertEquals(0., pv0.getVelocity().distance(pv1.getVelocity()), 2.8e-12);
140     }
141 
142     @Test
143     public void testNoReset() {
144         try {
145             GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
146             propagator.resetInitialState(propagator.getInitialState());
147             Assert.fail("an exception should have been thrown");
148         } catch (OrekitException oe) {
149             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
150         }
151         try {
152             GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
153             propagator.resetIntermediateState(propagator.getInitialState(), true);
154             Assert.fail("an exception should have been thrown");
155         } catch (OrekitException oe) {
156             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
157         }
158     }
159 
160     @Test
161     public void testDerivativesConsistency() {
162 
163         final Frame eme2000 = FramesFactory.getEME2000();
164         double errorP = 0;
165         double errorV = 0;
166         double errorA = 0;
167         GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
168         GNSSOrbitalElements elements = propagator.getOrbitalElements();
169         AbsoluteDate t0 = new GNSSDate(elements.getWeek(), 0.001 * elements.getTime(), SatelliteSystem.GALILEO).getDate();
170         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 600) {
171             final AbsoluteDate central = t0.shiftedBy(dt);
172             final PVCoordinates pv = propagator.getPVCoordinates(central, eme2000);
173             final double h = 10.0;
174             List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
175             for (int i = -3; i <= 3; ++i) {
176                 sample.add(propagator.getPVCoordinates(central.shiftedBy(i * h), eme2000));
177             }
178             final PVCoordinates interpolated =
179                             TimeStampedPVCoordinates.interpolate(central,
180                                                                  CartesianDerivativesFilter.USE_P,
181                                                                  sample);
182             errorP = FastMath.max(errorP, Vector3D.distance(pv.getPosition(), interpolated.getPosition()));
183             errorV = FastMath.max(errorV, Vector3D.distance(pv.getVelocity(), interpolated.getVelocity()));
184             errorA = FastMath.max(errorA, Vector3D.distance(pv.getAcceleration(), interpolated.getAcceleration()));
185         }
186         Assert.assertEquals(0.0, errorP, 1.5e-11);
187         Assert.assertEquals(0.0, errorV, 2.2e-7);
188         Assert.assertEquals(0.0, errorA, 4.9e-8);
189 
190     }
191 
192     @Test
193     public void testPosition() {
194         // Date of the Galileo orbital elements, 10 April 2019 at 09:30:00 UTC
195         final AbsoluteDate target = goe.getDate();
196         // Build the Galileo propagator
197         final GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
198         // Compute the PV coordinates at the date of the Galileo orbital elements
199         final PVCoordinates pv = propagator.getPVCoordinates(target, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
200         // Computed position
201         final Vector3D computedPos = pv.getPosition();
202         // Expected position (reference from IGS file WUM0MGXULA_20191010500_01D_15M_ORB.sp3)
203         final Vector3D expectedPos = new Vector3D(10487480.721, 17867448.753, -21131462.002);
204         Assert.assertEquals(0., Vector3D.distance(expectedPos, computedPos), 2.1);
205     }
206 
207     @Test
208     public void testIssue544() {
209         // Builds the GalileoPropagator from the almanac
210         final GNSSPropagator propagator = new GNSSPropagatorBuilder(goe).build();
211         // In order to test the issue, we voluntary set a Double.NaN value in the date.
212         final AbsoluteDate date0 = new AbsoluteDate(2010, 5, 7, 7, 50, Double.NaN, TimeScalesFactory.getUTC());
213         final PVCoordinates pv0 = propagator.propagateInEcef(date0);
214         // Verify that an infinite loop did not occur
215         Assert.assertEquals(Vector3D.NaN, pv0.getPosition());
216         Assert.assertEquals(Vector3D.NaN, pv0.getVelocity()); 
217     }
218 
219 }