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.models.earth.ionosphere;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.Decimal64Field;
24  import org.hipparchus.util.FastMath;
25  import org.junit.Assert;
26  import org.junit.Before;
27  import org.junit.Test;
28  import org.orekit.Utils;
29  import org.orekit.bodies.FieldGeodeticPoint;
30  import org.orekit.bodies.GeodeticPoint;
31  import org.orekit.bodies.OneAxisEllipsoid;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.frames.TopocentricFrame;
34  import org.orekit.gnss.Frequency;
35  import org.orekit.orbits.CartesianOrbit;
36  import org.orekit.orbits.FieldCartesianOrbit;
37  import org.orekit.orbits.FieldOrbit;
38  import org.orekit.orbits.Orbit;
39  import org.orekit.propagation.FieldSpacecraftState;
40  import org.orekit.propagation.SpacecraftState;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.time.TimeScalesFactory;
44  import org.orekit.utils.Constants;
45  import org.orekit.utils.FieldPVCoordinates;
46  import org.orekit.utils.IERSConventions;
47  import org.orekit.utils.PVCoordinates;
48  
49  /**
50   * Reference values for the tests are from : "European Union (2016). European GNSS (Galileo)
51   * Open Service-Ionospheric Correction Algorithm for Galileo Single Frequency Users. 1.2."
52   */
53  public class NeQuickModelTest {
54  
55      private double[] medium;
56      private double[] high;
57  
58      @Before
59      public void setUp() {
60          Utils.setDataRoot("regular-data");
61          high = new double[] {
62              236.831641, -0.39362878, 0.00402826613
63          };
64          medium = new double[] {
65              121.129893, 0.351254133, 0.0134635348
66          };
67          
68      }
69  
70      @Test
71      public void testHighSolarActivity() {
72  
73          // Model
74          final NeQuickModel model = new NeQuickModel(high);
75  
76          // Geodetic points
77          final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(82.49), FastMath.toRadians(297.66), 78.11);
78          final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(54.29), FastMath.toRadians(8.23), 20281546.18);
79  
80          // Date
81          final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 0, 0, 0, TimeScalesFactory.getUTC());
82  
83          // STEC
84          final double stec = model.stec(date, recP, satP);
85          Assert.assertEquals(20.40, stec, 0.09);
86      }
87  
88      @Test
89      public void testFieldHighSolarActivity() {
90          doTestFieldHighSolarActivity(Decimal64Field.getInstance());
91      }
92  
93      private <T extends CalculusFieldElement<T>> void doTestFieldHighSolarActivity(final Field<T> field) {
94  
95          // Zero
96          final T zero = field.getZero();
97  
98          // Model
99          final NeQuickModel model = new NeQuickModel(high);
100 
101         // Geodetic points
102         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(82.49)),
103                         zero.add(FastMath.toRadians(297.66)), zero.add(78.11));
104         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(54.29)),
105                         zero.add(FastMath.toRadians(8.23)), zero.add(20281546.18));
106 
107         // Date
108         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 0, 0, 0, TimeScalesFactory.getUTC());
109 
110         // STEC
111         final T stec = model.stec(date, recP, satP);
112         Assert.assertEquals(20.40, stec.getReal(), 0.09);
113     }
114 
115     @Test
116     public void testMediumSolarActivity() {
117 
118         // Model
119         final NeQuickModel model = new NeQuickModel(medium);
120 
121         // Geodetic points
122         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
123         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
124 
125         // Date
126         final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
127 
128         // STEC
129         final double stec = model.stec(date, recP, satP);
130         Assert.assertEquals(6.96, stec, 0.05);
131     }    
132 
133     @Test
134     public void testFieldMediumSolarActivity() {
135         doTestFieldMediumSolarActivity(Decimal64Field.getInstance());
136     }
137 
138     private <T extends CalculusFieldElement<T>> void doTestFieldMediumSolarActivity(final Field<T> field) {
139 
140         // Zero
141         final T zero = field.getZero();
142 
143         // Model
144         final NeQuickModel model = new NeQuickModel(medium);
145 
146         // Geodetic points
147         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(-31.80)),
148                         zero.add(FastMath.toRadians(115.89)), zero.add(12.78));
149         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(-14.31)),
150                         zero.add(FastMath.toRadians(124.09)), zero.add(20100697.90));
151 
152         // Date
153         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
154 
155         // STEC
156         final T stec = model.stec(date, recP, satP);
157         Assert.assertEquals(6.96, stec.getReal(), 0.05);
158     }   
159 
160     @Test
161     public void testDelay() {
162 
163         // Model
164         final NeQuickModel model = new NeQuickModel(medium);
165 
166         // Geodetic points
167         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
168         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
169 
170         // Date
171         final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
172 
173         // Earth
174         final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
175                                                                 Constants.WGS84_EARTH_FLATTENING,
176                                                                 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
177         // Satellite position
178         final Vector3D satPosInITRF    = ellipsoid.transform(satP);
179         final Vector3D satPosInEME2000 = ellipsoid.getBodyFrame().getTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
180 
181         // Spacecraft state
182         final PVCoordinates   pv      = new PVCoordinates(satPosInEME2000, new Vector3D(1.0, 1.0, 1.0));
183         final Orbit           orbit   = new CartesianOrbit(pv, FramesFactory.getEME2000(), date, Constants.WGS84_EARTH_MU);
184         final SpacecraftState state   = new SpacecraftState(orbit);
185 
186         final double delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
187                                              Frequency.G01.getMHzFrequency() * 1.0E6, model.getParameters());
188        
189         // Verify
190         Assert.assertEquals(1.13, delay, 0.01);
191     }
192 
193     @Test
194     public void testFieldDelay() {
195         doTestFieldDelay(Decimal64Field.getInstance());
196     }
197 
198     private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
199 
200         // Zero and One
201         final T zero = field.getZero();
202         final T one  = field.getOne();        
203 
204         // Model
205         final NeQuickModel model = new NeQuickModel(medium);
206 
207         // Geodetic points
208         final double recLat = FastMath.toRadians(-31.80);
209         final double recLon = FastMath.toRadians(115.89);
210         final double recAlt = 12.78;
211         final GeodeticPoint         recP = new GeodeticPoint(recLat, recLon, recAlt);
212         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(-14.31)),
213                         zero.add(FastMath.toRadians(124.09)), zero.add(20100697.90));
214 
215         // Date
216         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
217 
218         // Earth
219         final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
220                                                                 Constants.WGS84_EARTH_FLATTENING,
221                                                                 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
222         // Satellite position
223         final FieldVector3D<T> satPosInITRF    = ellipsoid.transform(satP);
224         final FieldVector3D<T> satPosInEME2000 = ellipsoid.getBodyFrame().getTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
225 
226         // Spacecraft state
227         final FieldPVCoordinates<T>   pv      = new FieldPVCoordinates<>(satPosInEME2000, new FieldVector3D<>(one, one, one));
228         final FieldOrbit<T>           orbit   = new FieldCartesianOrbit<>(pv, FramesFactory.getEME2000(), date, zero.add(Constants.WGS84_EARTH_MU));
229         final FieldSpacecraftState<T> state   = new FieldSpacecraftState<>(orbit);
230 
231         final T delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
232                                         Frequency.G01.getMHzFrequency() * 1.0E6, model.getParameters(field));
233        
234         // Verify
235         Assert.assertEquals(1.13, delay.getReal(), 0.01);
236     }
237 
238 }