1   /* Copyright 2002-2021 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.BeforeClass;
26  import org.junit.Test;
27  import org.orekit.Utils;
28  import org.orekit.errors.OrekitException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.frames.Frame;
31  import org.orekit.frames.FramesFactory;
32  import org.orekit.gnss.SatelliteSystem;
33  import org.orekit.propagation.analytical.gnss.data.BeidouAlmanac;
34  import org.orekit.propagation.analytical.gnss.data.BeidouNavigationMessage;
35  import org.orekit.propagation.analytical.gnss.data.GNSSOrbitalElements;
36  import org.orekit.time.AbsoluteDate;
37  import org.orekit.time.GNSSDate;
38  import org.orekit.time.TimeScalesFactory;
39  import org.orekit.utils.CartesianDerivativesFilter;
40  import org.orekit.utils.Constants;
41  import org.orekit.utils.IERSConventions;
42  import org.orekit.utils.PVCoordinates;
43  import org.orekit.utils.TimeStampedPVCoordinates;
44  
45  public class BeidouPropagatorTest {
46  
47      private static BeidouAlmanac almanac;
48  
49      @BeforeClass
50      public static void setUpBeforeClass() {
51          Utils.setDataRoot("gnss");
52  
53          // Almanac for satellite 18 for May 28th 2019
54          almanac = new BeidouAlmanac();
55          almanac.setPRN(18);
56          almanac.setWeek(694);
57          almanac.setTime(4096.0);
58          almanac.setSqrtA(6493.3076);
59          almanac.setE(0.00482368);
60          almanac.setI0(0.0, -0.01365602);
61          almanac.setOmega0(1.40069711);
62          almanac.setOmegaDot(-2.11437379e-9);
63          almanac.setPa(3.11461541);
64          almanac.setM0(-2.53029382);
65          almanac.setAf0(0.0001096725);
66          almanac.setAf1(7.27596e-12);
67          almanac.setHealth(0);
68          almanac.setDate(new GNSSDate(almanac.getWeek(), 1000.0 * almanac.getTime(), SatelliteSystem.BEIDOU).getDate());
69      }
70  
71      @Test
72      public void testBeidouCycle() {
73          // Builds the BeiDou propagator from the almanac
74          final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
75          // Intermediate verification
76          Assert.assertEquals(18,           almanac.getPRN());
77          Assert.assertEquals(0,            almanac.getHealth());
78          Assert.assertEquals(0.0001096725, almanac.getAf0(), 1.0e-15);
79          Assert.assertEquals(7.27596e-12,  almanac.getAf1(), 1.0e-15);
80          // Propagate at the BeiDou date and one BeiDou cycle later
81          final AbsoluteDate date0 = almanac.getDate();
82          final Vector3D p0 = propagator.propagateInEcef(date0).getPosition();
83          final double bdtCycleDuration = almanac.getCycleDuration();
84          final AbsoluteDate date1 = date0.shiftedBy(bdtCycleDuration);
85          final Vector3D p1 = propagator.propagateInEcef(date1).getPosition();
86  
87          // Checks
88          Assert.assertEquals(0., p0.distance(p1), 0.);
89      }
90  
91      @Test
92      public void testFrames() {
93          // Builds the BeiDou propagator from the almanac
94          final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
95          Assert.assertEquals("EME2000", propagator.getFrame().getName());
96          Assert.assertEquals(3.986004418e+14, almanac.getMu(), 1.0e6);
97          // Defines some date
98          final AbsoluteDate date = new AbsoluteDate(2016, 3, 3, 12, 0, 0., TimeScalesFactory.getUTC());
99          // Get PVCoordinates at the date in the ECEF
100         final PVCoordinates pv0 = propagator.propagateInEcef(date);
101         // Get PVCoordinates at the date in the ECEF
102         final PVCoordinates pv1 = propagator.getPVCoordinates(date, propagator.getECEF());
103 
104         // Checks
105         Assert.assertEquals(0., pv0.getPosition().distance(pv1.getPosition()), 3.3e-8);
106         Assert.assertEquals(0., pv0.getVelocity().distance(pv1.getVelocity()), 3.9e-12);
107     }
108 
109     @Test
110     public void testNoReset() {
111         try {
112             GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
113             propagator.resetInitialState(propagator.getInitialState());
114             Assert.fail("an exception should have been thrown");
115         } catch (OrekitException oe) {
116             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
117         }
118         try {
119             GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
120             propagator.resetIntermediateState(propagator.getInitialState(), true);
121             Assert.fail("an exception should have been thrown");
122         } catch (OrekitException oe) {
123             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
124         }
125     }
126 
127     @Test
128     public void testDerivativesConsistency() {
129 
130         final Frame eme2000 = FramesFactory.getEME2000();
131         double errorP = 0;
132         double errorV = 0;
133         double errorA = 0;
134         GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
135         GNSSOrbitalElements elements = propagator.getOrbitalElements();
136         AbsoluteDate t0 = new GNSSDate(elements.getWeek(), 0.001 * elements.getTime(), SatelliteSystem.BEIDOU).getDate();
137         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 600) {
138             final AbsoluteDate central = t0.shiftedBy(dt);
139             final PVCoordinates pv = propagator.getPVCoordinates(central, eme2000);
140             final double h = 10.0;
141             List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
142             for (int i = -3; i <= 3; ++i) {
143                 sample.add(propagator.getPVCoordinates(central.shiftedBy(i * h), eme2000));
144             }
145             final PVCoordinates interpolated =
146                             TimeStampedPVCoordinates.interpolate(central,
147                                                                  CartesianDerivativesFilter.USE_P,
148                                                                  sample);
149             errorP = FastMath.max(errorP, Vector3D.distance(pv.getPosition(), interpolated.getPosition()));
150             errorV = FastMath.max(errorV, Vector3D.distance(pv.getVelocity(), interpolated.getVelocity()));
151             errorA = FastMath.max(errorA, Vector3D.distance(pv.getAcceleration(), interpolated.getAcceleration()));
152         }
153 
154         Assert.assertEquals(0.0, errorP, 3.8e-9);
155         Assert.assertEquals(0.0, errorV, 8.0e-8);
156         Assert.assertEquals(0.0, errorA, 2.0e-8);
157 
158     }
159 
160     @Test
161     public void testPosition() {
162         // Initial BeiDou orbital elements (Ref: IGS)
163         final BeidouNavigationMessage boe = new BeidouNavigationMessage();
164         boe.setPRN(7);
165         boe.setWeek(713);
166         boe.setTime(284400.0);
167         boe.setSqrtA(6492.84515953064);
168         boe.setE(0.00728036486543715);
169         boe.setDeltaN(2.1815194404696853E-9);
170         boe.setI0(0.9065628903946735);
171         boe.setIDot(0.0);
172         boe.setOmega0(-0.6647664535282437);
173         boe.setOmegaDot(-3.136916379444212E-9);
174         boe.setPa(-2.6584351442773304);
175         boe.setM0(0.9614806010234702);
176         boe.setCuc(7.306225597858429E-6);
177         boe.setCus(-6.314832717180252E-6);
178         boe.setCrc(406.96875);
179         boe.setCrs(225.9375);
180         boe.setCic(-7.450580596923828E-9);
181         boe.setCis(-1.4062970876693726E-7);
182         boe.setDate(new GNSSDate(boe.getWeek(), 1000.0 * boe.getTime(), SatelliteSystem.BEIDOU).getDate());
183         // Date of the BeiDou orbital elements (GPStime - BDTtime = 14s)
184         final AbsoluteDate target = boe.getDate().shiftedBy(-14.0);
185         // Build the BeiDou propagator
186         final GNSSPropagator propagator = new GNSSPropagatorBuilder(boe).build();
187         // Compute the PV coordinates at the date of the BeiDou orbital elements
188         final PVCoordinates pv = propagator.getPVCoordinates(target, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
189         // Computed position
190         final Vector3D computedPos = pv.getPosition();
191         // Expected position (reference from sp3 file WUM0MGXULA_20192470700_01D_05M_ORB.SP33)
192         final Vector3D expectedPos = new Vector3D(-10260690.520, 24061180.795, -32837341.074);
193         Assert.assertEquals(0., Vector3D.distance(expectedPos, computedPos), 3.1);
194     }
195 
196     @Test
197     public void testIssue544() {
198         // Builds the BeidouPropagator from the almanac
199         final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
200         // In order to test the issue, we voluntary set a Double.NaN value in the date.
201         final AbsoluteDate date0 = new AbsoluteDate(2010, 5, 7, 7, 50, Double.NaN, TimeScalesFactory.getUTC());
202         final PVCoordinates pv0 = propagator.propagateInEcef(date0);
203         // Verify that an infinite loop did not occur
204         Assert.assertEquals(Vector3D.NaN, pv0.getPosition());
205         Assert.assertEquals(Vector3D.NaN, pv0.getVelocity());
206         
207     }
208 
209 }