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.forces.gravity.potential;
18  
19  import org.hipparchus.util.FastMath;
20  import org.hipparchus.util.Precision;
21  import org.junit.jupiter.api.Assertions;
22  import org.junit.jupiter.api.Test;
23  import org.orekit.Utils;
24  import org.orekit.errors.OrekitException;
25  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics;
26  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
27  import org.orekit.time.AbsoluteDate;
28  import org.orekit.utils.Constants;
29  
30  public class EGMFormatReaderTest {
31  
32      @Test
33      void testReadNormalizedEgm96() {
34          Utils.setDataRoot("potential");
35          GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96_to5.ascii", true));
36          NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
37          NormalizedSphericalHarmonics harmonics = provider.onDate(AbsoluteDate.FUTURE_INFINITY);
38          Assertions.assertEquals(TideSystem.TIDE_FREE, provider.getTideSystem());
39          Assertions.assertEquals( 0.957254173792E-06, harmonics.getNormalizedCnm(3, 0), 1.0e-15);
40          Assertions.assertEquals( 0.174971983203E-06, harmonics.getNormalizedCnm(5, 5), 1.0e-15);
41          Assertions.assertEquals( 0.0,                harmonics.getNormalizedSnm(4, 0), 1.0e-15);
42          Assertions.assertEquals( 0.308853169333E-06, harmonics.getNormalizedSnm(4, 4), 1.0e-15);
43          Assertions.assertEquals(-0.295301647654E-06, harmonics.getNormalizedCnm(5, 4), 1.0e-15);
44      }
45  
46      @Test
47      void testReadNormalizedEgm2008() {
48          Utils.setDataRoot("potential");
49          GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("EGM2008_to2190_TideFree.ascii", true));
50          NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
51          NormalizedSphericalHarmonics harmonics = provider.onDate(AbsoluteDate.FUTURE_INFINITY);
52          Assertions.assertEquals(TideSystem.TIDE_FREE, provider.getTideSystem());
53          Assertions.assertEquals( 0.957161207093473E-06, harmonics.getNormalizedCnm(3, 0), 1.0e-15);
54          Assertions.assertEquals( 0.174811795496002E-06, harmonics.getNormalizedCnm(5, 5), 1.0e-15);
55          Assertions.assertEquals( 0.0,                harmonics.getNormalizedSnm(4, 0), 1.0e-15);
56          Assertions.assertEquals( 0.308803882149194E-06, harmonics.getNormalizedSnm(4, 4), 1.0e-15);
57          Assertions.assertEquals(-0.295328761175629E-06, harmonics.getNormalizedCnm(5, 4), 1.0e-15);
58      }
59  
60      @Test
61      void testReadUnnormalized() {
62          Utils.setDataRoot("potential");
63          GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96_to5.ascii", true));
64          UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
65          UnnormalizedSphericalHarmonics harmonics = provider.onDate(AbsoluteDate.FUTURE_INFINITY);
66          Assertions.assertEquals(TideSystem.TIDE_FREE, provider.getTideSystem());
67          int maxUlps = 1;
68          checkValue(harmonics.getUnnormalizedCnm(3, 0), 3, 0, 0.957254173792E-06, maxUlps);
69          checkValue(harmonics.getUnnormalizedCnm(5, 5), 5, 5, 0.174971983203E-06, maxUlps);
70          checkValue(harmonics.getUnnormalizedSnm(4, 0), 4, 0, 0.0,                maxUlps);
71          checkValue(harmonics.getUnnormalizedSnm(4, 4), 4, 4, 0.308853169333E-06, maxUlps);
72  
73          double a = (-0.295301647654E-06);
74          double b = 9*8*7*6*5*4*3*2;
75          double c = 2*11/b;
76          double result = a*FastMath.sqrt(c);
77  
78          Assertions.assertEquals(result, harmonics.getUnnormalizedCnm(5, 4), 1.0e-20);
79  
80          a = -0.188560802735E-06;
81          b = 8*7*6*5*4*3*2;
82          c=2*9/b;
83          result = a*FastMath.sqrt(c);
84          Assertions.assertEquals(result, harmonics.getUnnormalizedCnm(4, 4), 1.0e-20);
85  
86          Assertions.assertEquals(1.0826266835531513e-3, -harmonics.getUnnormalizedCnm(2, 0), 1.0e-20);
87  
88      }
89  
90      @Test
91      void testReadLimits() {
92          Utils.setDataRoot("potential");
93          GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96_to5.ascii", true));
94          UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(3, 2);
95          UnnormalizedSphericalHarmonics harmonics = provider.onDate(null);
96          try {
97              harmonics.getUnnormalizedCnm(3, 3);
98              Assertions.fail("an exception should have been thrown");
99          } catch (OrekitException oe) {
100             // expected
101         } catch (Exception e) {
102             Assertions.fail("wrong exception caught: " + e.getLocalizedMessage());
103         }
104         try {
105             harmonics.getUnnormalizedCnm(4, 2);
106             Assertions.fail("an exception should have been thrown");
107         } catch (OrekitException oe) {
108             // expected
109         } catch (Exception e) {
110             Assertions.fail("wrong exception caught: " + e.getLocalizedMessage());
111         }
112         harmonics.getUnnormalizedCnm(3, 2);
113         Assertions.assertEquals(3, provider.getMaxDegree());
114         Assertions.assertEquals(2, provider.getMaxOrder());
115     }
116 
117     @Test
118     void testCorruptedFile1() {
119         Assertions.assertThrows(OrekitException.class, () -> {
120             Utils.setDataRoot("potential");
121             GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("corrupted-1-egm96_to5", false));
122             GravityFieldFactory.getUnnormalizedProvider(5, 5);
123         });
124     }
125 
126     @Test
127     void testCorruptedFile2() {
128         Assertions.assertThrows(OrekitException.class, () -> {
129             Utils.setDataRoot("potential");
130             GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("corrupted-2-egm96_to5", false));
131             GravityFieldFactory.getUnnormalizedProvider(5, 5);
132         });
133     }
134 
135     @Test
136     void testCorruptedFile3() {
137         Assertions.assertThrows(OrekitException.class, () -> {
138             Utils.setDataRoot("potential");
139             GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("corrupted-3-egm96_to5", false));
140             GravityFieldFactory.getUnnormalizedProvider(5, 5);
141         });
142     }
143 
144     @Test
145     void testZeroTidePattern1() {
146         Utils.setDataRoot("potential");
147         GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("dummy_egm2008", true));
148         Assertions.assertEquals(TideSystem.ZERO_TIDE,
149                             GravityFieldFactory.getUnnormalizedProvider(5, 5).getTideSystem());
150     }
151 
152     @Test
153     void testZeroTidePattern2() {
154         Utils.setDataRoot("potential");
155         GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("dummy_zerotide", true));
156         Assertions.assertEquals(TideSystem.ZERO_TIDE,
157                             GravityFieldFactory.getUnnormalizedProvider(5, 5).getTideSystem());
158     }
159 
160     @Test
161     void testEGM2008TideFreePattern() {
162         Utils.setDataRoot("potential");
163         GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("dummy_egm2008_tidefree", true));
164         Assertions.assertEquals(TideSystem.TIDE_FREE,
165                             GravityFieldFactory.getUnnormalizedProvider(5, 5).getTideSystem());
166     }
167 
168     @Test
169     void testWgs84CoefficientOverride()
170         {
171         final double epsilon = Precision.EPSILON;
172 
173         Utils.setDataRoot("potential");
174         EGMFormatReader egm96Reader = new EGMFormatReader("egm96_to5.ascii", true);
175         GravityFieldFactory.addPotentialCoefficientsReader(egm96Reader);
176         GravityFieldFactory.getNormalizedProvider(5, 5);
177         Assertions.assertEquals(Constants.EGM96_EARTH_EQUATORIAL_RADIUS, egm96Reader.getAe(), epsilon);
178         Assertions.assertEquals(Constants.EGM96_EARTH_MU, egm96Reader.getMu(), epsilon);
179 
180         Utils.setDataRoot("potential");
181         EGMFormatReader wgs84Egm96Reader = new EGMFormatReader("egm96_to5.ascii", true, true);
182         GravityFieldFactory.addPotentialCoefficientsReader(wgs84Egm96Reader);
183         GravityFieldFactory.getNormalizedProvider(5, 5);
184         Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, wgs84Egm96Reader.getAe(), epsilon);
185         Assertions.assertEquals(Constants.WGS84_EARTH_MU, wgs84Egm96Reader.getMu(), epsilon);
186 
187         Utils.setDataRoot("potential");
188         EGMFormatReader egm2008Reader = new EGMFormatReader("EGM2008_to2190_TideFree.ascii", true);
189         GravityFieldFactory.addPotentialCoefficientsReader(egm2008Reader);
190         GravityFieldFactory.getNormalizedProvider(5, 5);
191         Assertions.assertEquals(Constants.EGM96_EARTH_EQUATORIAL_RADIUS, egm2008Reader.getAe(), epsilon);
192         Assertions.assertEquals(Constants.EGM96_EARTH_MU, egm2008Reader.getMu(), epsilon);
193 
194         Utils.setDataRoot("potential");
195         EGMFormatReader wgs84Egm2008Reader = new EGMFormatReader("EGM2008_to2190_TideFree.ascii", true, true);
196         GravityFieldFactory.addPotentialCoefficientsReader(wgs84Egm2008Reader);
197         GravityFieldFactory.getNormalizedProvider(5, 5);
198         Assertions.assertEquals(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, wgs84Egm2008Reader.getAe(), epsilon);
199         Assertions.assertEquals(Constants.WGS84_EARTH_MU, wgs84Egm2008Reader.getMu(), epsilon);
200 
201     }
202 
203     private void checkValue(final double value, final int n, final int m,
204                             final double constant, final int maxUlps)
205         {
206         double factor = GravityFieldFactory.getUnnormalizationFactors(n, m)[n][m];
207         double normalized = factor * constant;
208         double epsilon = maxUlps * FastMath.ulp(normalized);
209         Assertions.assertEquals(normalized, value, epsilon);
210     }
211 
212 }