1   /* Copyright 2022-2025 Thales Alenia Space
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.iturp834;
18  
19  import org.hipparchus.analysis.differentiation.Gradient;
20  import org.hipparchus.util.FastMath;
21  import org.junit.jupiter.api.Assertions;
22  import org.junit.jupiter.api.Test;
23  import org.orekit.bodies.FieldGeodeticPoint;
24  import org.orekit.bodies.GeodeticPoint;
25  
26  public abstract class AbstractGridTest<T extends AbstractGrid> {
27  
28      @Test
29      public abstract void testMetadata();
30  
31      protected void doTestMetadata(final T grid,
32                                    final double expectedMinLatDeg, final double expectedMaxLatDeg, final int expectedNbLat,
33                                    final double expectedMinLonDeg, final double expectedMaxLonDeg, final int expectedNbLon) {
34  
35         Assertions.assertEquals(expectedNbLat,     grid.getLatitudeAxis().size());
36         Assertions.assertEquals(expectedMinLatDeg, FastMath.toDegrees(grid.getLatitudeAxis().node(0)), 1.0e-10);
37         Assertions.assertEquals(expectedMaxLatDeg, FastMath.toDegrees(grid.getLatitudeAxis().node(120)), 1.0e-10);
38  
39         Assertions.assertEquals(expectedNbLon,     grid.getLongitudeAxis().size());
40         Assertions.assertEquals(expectedMinLonDeg, FastMath.toDegrees(grid.getLongitudeAxis().node(0)), 1.0e-10);
41         Assertions.assertEquals(expectedMaxLonDeg, FastMath.toDegrees(grid.getLongitudeAxis().node(240)), 1.0e-10);
42  
43      }
44  
45      @Test
46      public abstract void testMinMax();
47  
48      protected void doTestMinMax(final T grid, final double expectedMin, final double expectedMax, double tol) {
49  
50          double min = Double.POSITIVE_INFINITY;
51          double max = Double.NEGATIVE_INFINITY;
52          for (int i = 0; i < grid.getLongitudeAxis().size(); ++i) {
53              final double lon = grid.getLongitudeAxis().node(i);
54              for (int j = 0; j < grid.getLatitudeAxis().size(); ++j) {
55                  final double lat = grid.getLatitudeAxis().node(j);
56                  final double hCell = grid.getCell(new GeodeticPoint(lat, lon, 0), 0.0).evaluate();
57                  min = FastMath.min(min, hCell);
58                  max = FastMath.max(max, hCell);
59              }
60          }
61          Assertions.assertEquals(expectedMin, min, tol);
62          Assertions.assertEquals(expectedMax, max, tol);
63  
64      }
65  
66      @Test
67      public abstract void testValue();
68  
69      protected void doTestValue(final T grid, final GeodeticPoint point, final double soy,
70                                 final double expectedValue, final double tol) {
71         Assertions.assertEquals(expectedValue, grid.getCell(point, soy).evaluate(), tol);
72      }
73  
74      @Test
75      public abstract void testGradient();
76  
77      protected void doTestGradient(final T grid, final GeodeticPoint point, final double soy,
78                                    final double tolVal, final double tolDer) {
79  
80          final Gradient                     soyG   = Gradient.variable(4, 0, soy);
81          final FieldGeodeticPoint<Gradient> pointG =
82              new FieldGeodeticPoint<>(Gradient.variable(4, 1, point.getLatitude()),
83                                       Gradient.variable(4, 2, point.getLongitude()),
84                                       Gradient.variable(4, 3, point.getAltitude()));
85          final Gradient gradient = grid.getCell(pointG, soyG).evaluate();
86  
87          Assertions.assertEquals(grid.getCell(point, soy).evaluate(), gradient.getValue(), tolVal);
88          Assertions.assertEquals(cellDerivative(grid, point, soy, 60.0,   0), gradient.getPartialDerivative(0), tolDer);
89          Assertions.assertEquals(cellDerivative(grid, point, soy, 0.0001, 1), gradient.getPartialDerivative(1), tolDer);
90          Assertions.assertEquals(cellDerivative(grid, point, soy, 0.0001, 2), gradient.getPartialDerivative(2), tolDer);
91          Assertions.assertEquals(cellDerivative(grid, point, soy, 0.0001, 3), gradient.getPartialDerivative(3), tolDer);
92  
93      }
94  
95      private double cellDerivative(final T grid,
96                                    final GeodeticPoint point, final double soy,
97                                    final double delta, final int index) {
98  
99          final double dt   = index == 0 ? delta : 0.0;
100         final double dlat = index == 1 ? delta : 0.0;
101         final double dlon = index == 2 ? delta : 0.0;
102         final double dh   = index == 3 ? delta : 0.0;
103         final double fM4h = shiftedCell(grid, point, soy, -4 * dt, -4 * dlat, -4 * dlon, -4 * dh);
104         final double fM3h = shiftedCell(grid, point, soy, -3 * dt, -3 * dlat, -3 * dlon, -3 * dh);
105         final double fM2h = shiftedCell(grid, point, soy, -2 * dt, -2 * dlat, -2 * dlon, -2 * dh);
106         final double fM1h = shiftedCell(grid, point, soy, -1 * dt, -1 * dlat, -1 * dlon, -1 * dh);
107         final double fP1h = shiftedCell(grid, point, soy,  1 * dt,  1 * dlat,  1 * dlon,  1 * dh);
108         final double fP2h = shiftedCell(grid, point, soy,  2 * dt,  2 * dlat,  2 * dlon,  2 * dh);
109         final double fP3h = shiftedCell(grid, point, soy,  3 * dt,  3 * dlat,  3 * dlon,  3 * dh);
110         final double fP4h = shiftedCell(grid, point, soy,  4 * dt,  4 * dlat,  4 * dlon,  4 * dh);
111 
112         // eight-points finite differences, the remaining error is -h⁸/630 d⁹f/dx⁹ + O(h¹⁰)
113         return (-3 * (fP4h - fM4h) + 32 * (fP3h - fM3h) - 168 * (fP2h - fM2h) + 672 * (fP1h - fM1h)) /
114                (840 * delta);
115 
116     }
117 
118     private double shiftedCell(final T grid,
119                                final GeodeticPoint point, final double soy,
120                                final double dt, final double dlat, final double dlon, final double dh) {
121         return grid.getCell(new GeodeticPoint(point.getLatitude()  + dlat,
122                                               point.getLongitude() + dlon,
123                                               point.getAltitude()  + dh),
124                             soy + dt).evaluate();
125     }
126 
127 }