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