1   /* Copyright 2002-2023 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 java.io.BufferedReader;
20  import java.io.IOException;
21  import java.io.InputStream;
22  import java.io.InputStreamReader;
23  import java.nio.charset.StandardCharsets;
24  import java.text.ParseException;
25  import java.util.Arrays;
26  import java.util.Locale;
27  import java.util.regex.Pattern;
28  
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.Precision;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.errors.OrekitMessages;
33  import org.orekit.utils.Constants;
34  
35  /**This reader is adapted to the EGM Format.
36   *
37   * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
38   *  which will determine which reader to use with the selected gravity field file.</p>
39   *
40   * @see GravityFields
41   * @author Fabien Maussion
42   */
43  public class EGMFormatReader extends PotentialCoefficientsReader {
44  
45      /** Pattern for delimiting regular expressions. */
46      private static final Pattern SEPARATOR = Pattern.compile("\\s+");
47  
48      /** Start degree and order for coefficients container. */
49      private static final int START_DEGREE_ORDER = 120;
50  
51      /** Flag for using WGS84 values for equatorial radius and central attraction coefficient. */
52      private final boolean useWgs84Coefficients;
53  
54      /** Simple constructor.
55       * @param supportedNames regular expression for supported files names
56       * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
57       */
58      public EGMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
59          this(supportedNames, missingCoefficientsAllowed, false);
60      }
61  
62      /**
63       * Simple constructor that allows overriding 'standard' EGM96 ae and mu with
64       * WGS84 variants.
65       *
66       * @param supportedNames regular expression for supported files names
67       * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
68       * @param useWgs84Coefficients if true, the WGS84 values will be used for equatorial radius
69       * and central attraction coefficient
70       */
71      public EGMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed,
72                             final boolean useWgs84Coefficients) {
73          super(supportedNames, missingCoefficientsAllowed, null);
74          this.useWgs84Coefficients = useWgs84Coefficients;
75      }
76  
77  
78      /** {@inheritDoc} */
79      public void loadData(final InputStream input, final String name)
80          throws IOException, ParseException, OrekitException {
81  
82          // reset the indicator before loading any data
83          setReadComplete(false);
84  
85          // both EGM96 and EGM2008 use the same values for ae and mu
86          // if a new EGM model changes them, we should have some selection logic
87          // based on file name (a better way would be to have the data in the
88          // file...)
89          if (this.useWgs84Coefficients) {
90              setAe(Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
91              setMu(Constants.WGS84_EARTH_MU);
92          } else {
93              setAe(Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
94              setMu(Constants.EGM96_EARTH_MU);
95          }
96  
97          final String lowerCaseName = name.toLowerCase(Locale.US);
98          if (lowerCaseName.contains("2008") || lowerCaseName.contains("zerotide")) {
99              setTideSystem(TideSystem.ZERO_TIDE);
100         } else {
101             setTideSystem(TideSystem.TIDE_FREE);
102         }
103 
104         Container container = new Container(START_DEGREE_ORDER, START_DEGREE_ORDER,
105                                             missingCoefficientsAllowed() ? 0.0 : Double.NaN);
106         boolean okFields = true;
107         int       maxDegree  = -1;
108         int       maxOrder   = -1;
109         int lineNumber = 0;
110         String line = null;
111         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
112             for (line = r.readLine(); okFields && line != null; line = r.readLine()) {
113                 lineNumber++;
114                 if (line.length() >= 15) {
115 
116                     // get the fields defining the current potential terms
117                     final String[] tab = SEPARATOR.split(line.trim());
118                     if (tab.length != 6) {
119                         okFields = false;
120                     }
121 
122                     final int i = Integer.parseInt(tab[0]);
123                     final int j = Integer.parseInt(tab[1]);
124                     if (i < 0 || j < 0) {
125                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
126                                                   lineNumber, name, line);
127                     }
128 
129                     if (i <= getMaxParseDegree() && j <= getMaxParseOrder()) {
130 
131                         while (!container.flattener.withinRange(i, j)) {
132                             // we need to resize the container
133                             container = container.resize(container.flattener.getDegree() * 2,
134                                                          container.flattener.getOrder() * 2);
135                         }
136 
137                         parseCoefficient(tab[2], container.flattener, container.c, i, j, "C", name);
138                         parseCoefficient(tab[3], container.flattener, container.s, i, j, "S", name);
139                         maxDegree = FastMath.max(maxDegree, i);
140                         maxOrder  = FastMath.max(maxOrder,  j);
141 
142                     }
143 
144                 }
145             }
146         } catch (NumberFormatException nfe) {
147             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
148                                       lineNumber, name, line);
149         }
150 
151         if (missingCoefficientsAllowed() && getMaxParseDegree() > 0 && getMaxParseOrder() > 0) {
152             // ensure at least the (0, 0) element is properly set
153             if (Precision.equals(container.c[container.flattener.index(0, 0)], 0.0, 0)) {
154                 container.c[container.flattener.index(0, 0)] = 1.0;
155             }
156         }
157 
158         if (!(okFields && maxDegree >= 0)) {
159             String loaderName = getClass().getName();
160             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
161             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
162                                       name, loaderName);
163         }
164 
165         container = container.resize(maxDegree, maxOrder);
166         setRawCoefficients(true, container.flattener, container.c, container.s, name);
167         setReadComplete(true);
168 
169     }
170 
171     /** Get a provider for read spherical harmonics coefficients.
172      * <p>
173      * EGM fields don't include time-dependent parts, so this method returns
174      * directly a constant provider.
175      * </p>
176      * @param wantNormalized if true, the provider will provide normalized coefficients,
177      * otherwise it will provide un-normalized coefficients
178      * @param degree maximal degree
179      * @param order maximal order
180      * @return a new provider
181      * @since 6.0
182      */
183     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
184                                                      final int degree, final int order) {
185         return getBaseProvider(wantNormalized, degree, order);
186     }
187 
188     /** Temporary container for reading coefficients.
189      * @since 11.1
190      */
191     private static class Container {
192 
193         /** Converter from triangular to flat form. */
194         private final Flattener flattener;
195 
196         /** Cosine coefficients. */
197         private final double[] c;
198 
199         /** Sine coefficients. */
200         private final double[] s;
201 
202         /** Initial value for coefficients. */
203         private final double initialValue;
204 
205         /** Build a container with given degree and order.
206          * @param degree degree of the container
207          * @param order order of the container
208          * @param initialValue initial value for coefficients
209          */
210         Container(final int degree, final int order, final double initialValue) {
211             this.flattener    = new Flattener(degree, order);
212             this.c            = new double[flattener.arraySize()];
213             this.s            = new double[flattener.arraySize()];
214             this.initialValue = initialValue;
215             Arrays.fill(c, initialValue);
216             Arrays.fill(s, initialValue);
217         }
218 
219         /** Build a resized container.
220          * @param degree degree of the resized container
221          * @param order order of the resized container
222          * @return resized container
223          */
224         Container resize(final int degree, final int order) {
225             final Container resized = new Container(degree, order, initialValue);
226             for (int n = 0; n <= degree; ++n) {
227                 for (int m = 0; m <= FastMath.min(n, order); ++m) {
228                     if (flattener.withinRange(n, m)) {
229                         final int rIndex = resized.flattener.index(n, m);
230                         final int index  = flattener.index(n, m);
231                         resized.c[rIndex] = c[index];
232                         resized.s[rIndex] = s[index];
233                     }
234                 }
235             }
236             return resized;
237         }
238 
239     }
240 
241 }