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.models.earth.troposphere;
18  
19  import org.hipparchus.analysis.differentiation.Gradient;
20  import org.hipparchus.util.FastMath;
21  import org.hipparchus.util.Precision;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.BeforeEach;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.Utils;
26  import org.orekit.bodies.FieldGeodeticPoint;
27  import org.orekit.bodies.GeodeticPoint;
28  import org.orekit.time.AbsoluteDate;
29  import org.orekit.time.FieldAbsoluteDate;
30  import org.orekit.time.TimeScalesFactory;
31  import org.orekit.utils.FieldTrackingCoordinates;
32  import org.orekit.utils.TrackingCoordinates;
33  
34  public abstract class AbstractMappingFunctionTest<T extends TroposphereMappingFunction> {
35  
36      protected abstract T buildMappingFunction();
37  
38      // default values for doTestMappingFactors
39      protected AbsoluteDate defaultDate;
40      protected GeodeticPoint defaultPoint;
41      protected TrackingCoordinates defaultTrackingCoordinates;
42  
43      @Test
44      public abstract void testMappingFactors();
45  
46      protected void doTestMappingFactors(final AbsoluteDate date, final GeodeticPoint point,
47                                          final TrackingCoordinates trackingCoordinates,
48                                          final double expectedHydro, final double expectedWet) {
49          final T model = buildMappingFunction();
50          final double[] computedMapping = model.mappingFactors(trackingCoordinates, point,
51                                                                date);
52          Assertions.assertEquals(expectedHydro, computedMapping[0], 1.0e-3);
53          Assertions.assertEquals(expectedWet,   computedMapping[1], 1.0e-3);
54      }
55  
56      @Test
57      public void testFixedHeight() {
58          final AbsoluteDate date = new AbsoluteDate();
59          final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(45.0), 350.0);
60          T model = buildMappingFunction();
61          double[] lastFactors = new double[] {
62              Double.MAX_VALUE,
63              Double.MAX_VALUE
64          };
65          // mapping functions shall decline with increasing elevation angle
66          for (double elev = 10d; elev < 90d; elev += 8d) {
67              final double[] factors = model.mappingFactors(new TrackingCoordinates(0.0, FastMath.toRadians(elev), 0.0),
68                                                            point,
69                                                            date);
70              Assertions.assertTrue(Precision.compareTo(factors[0], lastFactors[0], 1.0e-6) < 0);
71              Assertions.assertTrue(Precision.compareTo(factors[1], lastFactors[1], 1.0e-6) < 0);
72              lastFactors[0] = factors[0];
73              lastFactors[1] = factors[1];
74          }
75      }
76  
77      @Test
78      public abstract void testDerivatives();
79  
80      protected void doTestDerivatives(final double tolValue,
81                                       final double tolTimeDerivative,
82                                       final double tolAzimuthDerivative,
83                                       final double tolElevationDerivative,
84                                       final double tolRangeDerivative) {
85  
86          T model = buildMappingFunction();
87  
88          final AbsoluteDate dateD = new AbsoluteDate(2004, 11, 16, 10, 54, 17.0, TimeScalesFactory.getUTC());
89          final GeodeticPoint kuroishiD =
90              new GeodeticPoint(FastMath.toRadians(40.64236), FastMath.toRadians(140.59590), 59.0);
91          final TrackingCoordinates trackingCoordinatesD = new TrackingCoordinates(FastMath.toRadians(210.0),
92                                                                                   FastMath.toRadians(15.0),
93                                                                                   1.4e6);
94  
95          // gradients with respect to date, azimuth, elevation and range
96          final FieldAbsoluteDate<Gradient> dateG = new FieldAbsoluteDate<>(dateD, Gradient.variable(4, 0, 0.0));
97          final FieldGeodeticPoint<Gradient> kuroishiG =
98              new FieldGeodeticPoint<>(dateG.getField(), kuroishiD);
99          final FieldTrackingCoordinates<Gradient> trackingCoordinatesG =
100             new FieldTrackingCoordinates<>(Gradient.variable(4, 1, trackingCoordinatesD.getAzimuth()),
101                                            Gradient.variable(4, 2, trackingCoordinatesD.getElevation()),
102                                            Gradient.variable(4, 3, trackingCoordinatesD.getRange()));
103 
104         final Gradient[] finiteDifferences = mappingFactorsGradient(model,
105                                                                     trackingCoordinatesD, kuroishiD, dateD,
106                                                                     1000.0, 1.0e-2, 1.0e-2, 100.0);
107         final Gradient[] direct = model.mappingFactors(trackingCoordinatesG, kuroishiG, dateG);
108 
109         // values
110         Assertions.assertEquals(finiteDifferences[0].getValue(), direct[0].getValue(), tolValue);
111         Assertions.assertEquals(finiteDifferences[1].getValue(), direct[1].getValue(), tolValue);
112 
113         // derivatives with respect to date
114         Assertions.assertEquals(finiteDifferences[0].getGradient()[0], direct[0].getGradient()[0], tolTimeDerivative);
115         Assertions.assertEquals(finiteDifferences[1].getGradient()[0], direct[1].getGradient()[0], tolTimeDerivative);
116 
117         // derivatives with respect to azimuth
118         Assertions.assertEquals(finiteDifferences[0].getGradient()[1], direct[0].getGradient()[1], tolAzimuthDerivative);
119         Assertions.assertEquals(finiteDifferences[1].getGradient()[1], direct[1].getGradient()[1], tolAzimuthDerivative);
120 
121         // derivatives with respect to elevation
122         Assertions.assertEquals(finiteDifferences[0].getGradient()[2], direct[0].getGradient()[2], tolElevationDerivative);
123         Assertions.assertEquals(finiteDifferences[1].getGradient()[2], direct[1].getGradient()[2], tolElevationDerivative);
124 
125         // derivatives with respect to range
126         Assertions.assertEquals(finiteDifferences[0].getGradient()[3], direct[0].getGradient()[3], tolRangeDerivative);
127         Assertions.assertEquals(finiteDifferences[1].getGradient()[3], direct[1].getGradient()[3], tolRangeDerivative);
128 
129     }
130 
131     private Gradient[] mappingFactorsGradient(final T model,
132                                               final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
133                                               final AbsoluteDate date,
134                                               final double dt, final double da, final double de, final double dr) {
135         return new Gradient[] {
136             new Gradient(model.mappingFactors(trackingCoordinates, point, date)[0],
137                          mappingFactorsDerivative(model, trackingCoordinates, point, date, dt, 0)[0],
138                          mappingFactorsDerivative(model, trackingCoordinates, point, date, da, 1)[0],
139                          mappingFactorsDerivative(model, trackingCoordinates, point, date, de, 2)[0],
140                          mappingFactorsDerivative(model, trackingCoordinates, point, date, dr, 3)[0]),
141             new Gradient(model.mappingFactors(trackingCoordinates, point, date)[1],
142                          mappingFactorsDerivative(model, trackingCoordinates, point, date, dt, 0)[1],
143                          mappingFactorsDerivative(model, trackingCoordinates, point, date, da, 1)[1],
144                          mappingFactorsDerivative(model, trackingCoordinates, point, date, de, 2)[1],
145                          mappingFactorsDerivative(model, trackingCoordinates, point, date, dr, 3)[1])
146         };
147     }
148 
149     private double[] mappingFactorsDerivative(final T model,
150                                               final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
151                                               final AbsoluteDate date,
152                                               final double delta, final int index) {
153 
154         final double dt = index == 0 ? delta : 0.0;
155         final double da = index == 1 ? delta : 0.0;
156         final double de = index == 2 ? delta : 0.0;
157         final double dr = index == 3 ? delta : 0.0;
158         final double[] mM4h = shiftedMappingFactors(model, trackingCoordinates, point, date,
159                                                     -4 * dt, -4 * da, -4 * de, -4 * dr);
160         final double[] mM3h = shiftedMappingFactors(model, trackingCoordinates, point, date,
161                                                     -3 * dt, -3 * da, -3 * de, -3 * dr);
162         final double[] mM2h = shiftedMappingFactors(model, trackingCoordinates, point, date,
163                                                     -2 * dt, -2 * da, -2 * de, -2 * dr);
164         final double[] mM1h = shiftedMappingFactors(model, trackingCoordinates, point, date,
165                                                     -1 * dt, -1 * da, -1 * de, -1 * dr);
166         final double[] mP1h = shiftedMappingFactors(model, trackingCoordinates, point, date,
167                                                      1 * dt,  1 * da,  1 * de,  1 * dr);
168         final double[] mP2h = shiftedMappingFactors(model, trackingCoordinates, point, date,
169                                                      2 * dt,  2 * da,  2 * de,  2 * dr);
170         final double[] mP3h = shiftedMappingFactors(model, trackingCoordinates, point, date,
171                                                      3 * dt,  3 * da,  3 * de,  3 * dr);
172         final double[] mP4h = shiftedMappingFactors(model, trackingCoordinates, point, date,
173                                                      4 * dt,  4 * da,  4 * de,  4 * dr);
174 
175         return new double[] {
176             differential8(mM4h[0], mM3h[0], mM2h[0], mM1h[0], mP1h[0], mP2h[0], mP3h[0], mP4h[0], delta),
177             differential8(mM4h[1], mM3h[1], mM2h[1], mM1h[1], mP1h[1], mP2h[1], mP3h[1], mP4h[1], delta)
178         };
179 
180     }
181 
182     private double[] shiftedMappingFactors(final T model, final TrackingCoordinates trackingCoordinates,
183                                            final GeodeticPoint point, final AbsoluteDate date,
184                                            final double dt, final double da, final double de, final double dr) {
185         return model.mappingFactors(new TrackingCoordinates(trackingCoordinates.getAzimuth()   + da,
186                                                             trackingCoordinates.getElevation() + de,
187                                                             trackingCoordinates.getRange()     + dr),
188                                     point, date.shiftedBy(dt));
189     }
190 
191     private double differential8(final double fM4h, final double fM3h, final double fM2h, final double fM1h,
192                                  final double fP1h, final double fP2h, final double fP3h, final double fP4h,
193                                  final double h) {
194 
195         // eight-points finite differences
196         // the remaining error is -h^8/630 d^9f/dx^9 + O(h^10)
197         return (-3 * (fP4h - fM4h) + 32 * (fP3h - fM3h) - 168 * (fP2h - fM2h) + 672 * (fP1h - fM1h)) / (840 * h);
198 
199     }
200 
201     @BeforeEach
202     public void setUp() {
203         Utils.setDataRoot("regular-data");
204         defaultDate                = new AbsoluteDate(1994, 1, 1, TimeScalesFactory.getUTC());
205         defaultPoint               = new GeodeticPoint(FastMath.toRadians(48.0), FastMath.toRadians(0.20), 68.0);
206         defaultTrackingCoordinates = new TrackingCoordinates(0.0, FastMath.toRadians(5.0), 0.0);
207     }
208 
209 }