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.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.GNSSOrbitalElements;
34  import org.orekit.propagation.analytical.gnss.data.QZSSAlmanac;
35  import org.orekit.propagation.analytical.gnss.data.QZSSNavigationMessage;
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 QZSSPropagatorTest {
46  
47      private static QZSSAlmanac almanac;
48  
49      @BeforeClass
50      public static void setUpBeforeClass() {
51          Utils.setDataRoot("gnss");
52  
53          // Almanac for satellite 193 for May 27th 2019 (q201914.alm)
54          almanac = new QZSSAlmanac();
55          almanac.setPRN(193);
56          almanac.setWeek(7);
57          almanac.setTime(348160.0);
58          almanac.setSqrtA(6493.145996);
59          almanac.setE(7.579761505E-02);
60          almanac.setI0(0.7201680272);
61          almanac.setOmega0(-1.643310999);
62          almanac.setOmegaDot(-3.005839491E-09);
63          almanac.setPa(-1.561775201);
64          almanac.setM0(-4.050903957E-01);
65          almanac.setAf0(-2.965927124E-04);
66          almanac.setAf1(7.275957614E-12);
67          almanac.setHealth(0);
68          almanac.setDate(new GNSSDate(almanac.getWeek(), 1000.0 * almanac.getTime(), SatelliteSystem.QZSS).getDate());
69  
70      }
71  
72      @Test
73      public void testQZSSCycle() {
74          // Builds the QZSS propagator from the almanac
75          final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
76          // Propagate at the QZSS date and one QZSS cycle later
77          final AbsoluteDate date0 = almanac.getDate();
78          final Vector3D p0 = propagator.propagateInEcef(date0).getPosition();
79          final double bdtCycleDuration = almanac.getCycleDuration();
80          final AbsoluteDate date1 = date0.shiftedBy(bdtCycleDuration);
81          final Vector3D p1 = propagator.propagateInEcef(date1).getPosition();
82  
83          // Checks
84          Assert.assertEquals(0., p0.distance(p1), 0.);
85      }
86  
87      @Test
88      public void testFrames() {
89          // Builds the QZSS propagator from the almanac
90          final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
91          Assert.assertEquals("EME2000", propagator.getFrame().getName());
92          Assert.assertEquals(3.986005e+14, almanac.getMu(), 1.0e6);
93          // Defines some date
94          final AbsoluteDate date = new AbsoluteDate(2016, 3, 3, 12, 0, 0., TimeScalesFactory.getUTC());
95          // Get PVCoordinates at the date in the ECEF
96          final PVCoordinates pv0 = propagator.propagateInEcef(date);
97          // Get PVCoordinates at the date in the ECEF
98          final PVCoordinates pv1 = propagator.getPVCoordinates(date, propagator.getECEF());
99  
100         // Checks
101         Assert.assertEquals(0., pv0.getPosition().distance(pv1.getPosition()), 3.3e-8);
102         Assert.assertEquals(0., pv0.getVelocity().distance(pv1.getVelocity()), 3.9e-12);
103     }
104 
105     @Test
106     public void testNoReset() {
107         try {
108             GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
109             propagator.resetInitialState(propagator.getInitialState());
110             Assert.fail("an exception should have been thrown");
111         } catch (OrekitException oe) {
112             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
113         }
114         try {
115             GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
116             propagator.resetIntermediateState(propagator.getInitialState(), true);
117             Assert.fail("an exception should have been thrown");
118         } catch (OrekitException oe) {
119             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
120         }
121     }
122 
123     @Test
124     public void testDerivativesConsistency() {
125 
126         final Frame eme2000 = FramesFactory.getEME2000();
127         double errorP = 0;
128         double errorV = 0;
129         double errorA = 0;
130         GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
131         GNSSOrbitalElements elements = propagator.getOrbitalElements();
132         AbsoluteDate t0 = new GNSSDate(elements.getWeek(), 0.001 * elements.getTime(), SatelliteSystem.QZSS).getDate();
133         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 600) {
134             final AbsoluteDate central = t0.shiftedBy(dt);
135             final PVCoordinates pv = propagator.getPVCoordinates(central, eme2000);
136             final double h = 10.0;
137             List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
138             for (int i = -3; i <= 3; ++i) {
139                 sample.add(propagator.getPVCoordinates(central.shiftedBy(i * h), eme2000));
140             }
141             final PVCoordinates interpolated =
142                             TimeStampedPVCoordinates.interpolate(central,
143                                                                  CartesianDerivativesFilter.USE_P,
144                                                                  sample);
145             errorP = FastMath.max(errorP, Vector3D.distance(pv.getPosition(), interpolated.getPosition()));
146             errorV = FastMath.max(errorV, Vector3D.distance(pv.getVelocity(), interpolated.getVelocity()));
147             errorA = FastMath.max(errorA, Vector3D.distance(pv.getAcceleration(), interpolated.getAcceleration()));
148         }
149 
150         Assert.assertEquals(0.0, errorP, 3.8e-9);
151         Assert.assertEquals(0.0, errorV, 8.4e-8);
152         Assert.assertEquals(0.0, errorA, 2.1e-8);
153 
154     }
155 
156     @Test
157     public void testPosition() {
158         // Initial QZSS orbital elements (Ref: IGS)
159         final QZSSNavigationMessage qoe = new QZSSNavigationMessage();
160         qoe.setPRN(195);
161         qoe.setWeek(21);
162         qoe.setTime(226800.0);
163         qoe.setSqrtA(6493.226968765259);
164         qoe.setE(0.07426900835707784);
165         qoe.setDeltaN(4.796628370253418E-10);
166         qoe.setI0(0.7116940567084221);
167         qoe.setIDot(4.835915721014987E-10);
168         qoe.setOmega0(0.6210371871830609);
169         qoe.setOmegaDot(-8.38963517626603E-10);
170         qoe.setPa(-1.5781555771543598);
171         qoe.setM0(1.077008903618136);
172         qoe.setCuc(-8.8568776845932E-6);
173         qoe.setCus(1.794286072254181E-5);
174         qoe.setCrc(-344.03125);
175         qoe.setCrs(-305.6875);
176         qoe.setCic(1.2032687664031982E-6);
177         qoe.setCis(-2.6728957891464233E-6);
178         qoe.setDate(new GNSSDate(qoe.getWeek(), 1000.0 * qoe.getTime(), SatelliteSystem.QZSS).getDate());
179         // Date of the QZSS orbital elements
180         final AbsoluteDate target = qoe.getDate();
181         // Build the QZSS propagator
182         final GNSSPropagator propagator = new GNSSPropagatorBuilder(qoe).build();
183         // Compute the PV coordinates at the date of the QZSS orbital elements
184         final PVCoordinates pv = propagator.getPVCoordinates(target, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
185         // Computed position
186         final Vector3D computedPos = pv.getPosition();
187         // Expected position (reference from QZSS sp3 file qzu20693_00.sp3)
188         final Vector3D expectedPos = new Vector3D(-35047225.493, 18739632.916, -9522204.569);
189         Assert.assertEquals(0., Vector3D.distance(expectedPos, computedPos), 0.7);
190     }
191 
192     @Test
193     public void testIssue544() {
194         // Builds the QZSSPropagator from the almanac
195         final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanac).build();
196         // In order to test the issue, we voluntary set a Double.NaN value in the date.
197         final AbsoluteDate date0 = new AbsoluteDate(2010, 5, 7, 7, 50, Double.NaN, TimeScalesFactory.getUTC());
198         final PVCoordinates pv0 = propagator.propagateInEcef(date0);
199         // Verify that an infinite loop did not occur
200         Assert.assertEquals(Vector3D.NaN, pv0.getPosition());
201         Assert.assertEquals(Vector3D.NaN, pv0.getVelocity());
202         
203     }
204 }