1   /* Copyright 2002-2021 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.ArrayList;
26  import java.util.HashMap;
27  import java.util.List;
28  import java.util.Map;
29  import java.util.regex.Pattern;
30  
31  import org.hipparchus.util.FastMath;
32  import org.hipparchus.util.Precision;
33  import org.orekit.annotation.DefaultDataContext;
34  import org.orekit.data.DataContext;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.errors.OrekitMessages;
37  import org.orekit.errors.OrekitParseException;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.time.DateComponents;
40  import org.orekit.time.TimeScale;
41  import org.orekit.utils.Constants;
42  
43  /** Reader for the ICGEM gravity field format.
44   *
45   * <p>This format is used to describe the gravity field of EIGEN models
46   * published by the GFZ Potsdam since 2004. It is described in Franz
47   * Barthelmes and Christoph F&ouml;rste paper: "the ICGEM-format".
48   * The 2006-02-28 version of this paper can be found <a
49   * href="http://op.gfz-potsdam.de/grace/results/grav/g005_ICGEM-Format.pdf">here</a>
50   * and the 2011-06-07 version of this paper can be found <a
51   * href="http://icgem.gfz-potsdam.de/ICGEM-Format-2011.pdf">here</a>.
52   * These versions differ in time-dependent coefficients, which are linear-only prior
53   * to 2011 (up to eigen-5 model) and have also harmonic effects after that date
54   * (starting with eigen-6 model). Both versions are supported by the class.</p>
55   * <p>
56   * This reader uses relaxed check on the gravity constant key so any key ending
57   * in gravity_constant is accepted and not only earth_gravity_constant as specified
58   * in the previous documents. This allows to read also non Earth gravity fields
59   * as found in <a href="http://icgem.gfz-potsdam.de/tom_celestial">ICGEM
60   * - Gravity Field Models of other Celestial Bodies</a> page to be read.
61   * </p>
62   * <p>
63   * In order to simplify implementation, some design choices have been made: the
64   * reference date and the periods of harmonic pulsations are stored globally and
65   * not on a per-coefficient basis. This has the following implications:
66   * </p>
67   * <ul>
68   *   <li>
69   *     all time-stamped coefficients must share the same reference date, otherwise
70   *     an error will be triggered during parsing,
71   *   </li>
72   *   <li>
73   *     in order to avoid too large memory and CPU consumption, only a few different
74   *     periods should appear in the file,
75   *   </li>
76   *   <li>
77   *     only one occurrence of each coefficient may appear in the file, otherwise
78   *     an error will be triggered during parsing. Multiple occurrences with different
79   *     time stamps are forbidden (both because they correspond to a duplicated entry
80   *     and because they define two different reference dates as per previous design
81   *     choice).
82   *   </li>
83   * </ul>
84   *
85   * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
86   *  which will determine which reader to use with the selected gravity field file.</p>
87   *
88   * @see GravityFields
89   * @author Luc Maisonobe
90   */
91  public class ICGEMFormatReader extends PotentialCoefficientsReader {
92  
93      /** Product type. */
94      private static final String PRODUCT_TYPE            = "product_type";
95  
96      /** Gravity field product type. */
97      private static final String GRAVITY_FIELD           = "gravity_field";
98  
99      /** Gravity constant marker. */
100     private static final String GRAVITY_CONSTANT        = "gravity_constant";
101 
102     /** Reference radius. */
103     private static final String REFERENCE_RADIUS        = "radius";
104 
105     /** Max degree. */
106     private static final String MAX_DEGREE              = "max_degree";
107 
108     /** Tide system indicator. */
109     private static final String TIDE_SYSTEM_INDICATOR   = "tide_system";
110 
111     /** Indicator value for zero-tide system. */
112     private static final String ZERO_TIDE               = "zero_tide";
113 
114     /** Indicator value for tide-free system. */
115     private static final String TIDE_FREE               = "tide_free";
116 
117     /** Indicator value for unknown tide system. */
118     private static final String TIDE_UNKNOWN            = "unknown";
119 
120     /** Normalization indicator. */
121     private static final String NORMALIZATION_INDICATOR = "norm";
122 
123     /** Indicator value for normalized coefficients. */
124     private static final String NORMALIZED              = "fully_normalized";
125 
126     /** Indicator value for un-normalized coefficients. */
127     private static final String UNNORMALIZED            = "unnormalized";
128 
129     /** End of header marker. */
130     private static final String END_OF_HEADER           = "end_of_head";
131 
132     /** Gravity field coefficient. */
133     private static final String GFC                     = "gfc";
134 
135     /** Time stamped gravity field coefficient. */
136     private static final String GFCT                    = "gfct";
137 
138     /** Gravity field coefficient first time derivative. */
139     private static final String DOT                     = "dot";
140 
141     /** Gravity field coefficient trend. */
142     private static final String TRND                    = "trnd";
143 
144     /** Gravity field coefficient sine amplitude. */
145     private static final String ASIN                    = "asin";
146 
147     /** Gravity field coefficient cosine amplitude. */
148     private static final String ACOS                    = "acos";
149 
150     /** Pattern for delimiting regular expressions. */
151     private static final Pattern SEPARATOR = Pattern.compile("\\s+");
152 
153     /** Indicator for normalized coefficients. */
154     private boolean normalized;
155 
156     /** Reference date. */
157     private AbsoluteDate referenceDate;
158 
159     /** Secular trend of the cosine coefficients. */
160     private final List<List<Double>> cTrend;
161 
162     /** Secular trend of the sine coefficients. */
163     private final List<List<Double>> sTrend;
164 
165     /** Cosine part of the cosine coefficients pulsation. */
166     private final Map<Double, List<List<Double>>> cCos;
167 
168     /** Sine part of the cosine coefficients pulsation. */
169     private final Map<Double, List<List<Double>>> cSin;
170 
171     /** Cosine part of the sine coefficients pulsation. */
172     private final Map<Double, List<List<Double>>> sCos;
173 
174     /** Sine part of the sine coefficients pulsation. */
175     private final Map<Double, List<List<Double>>> sSin;
176 
177     /** Simple constructor.
178      *
179      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
180      *
181      * @param supportedNames regular expression for supported files names
182      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
183      * @see #ICGEMFormatReader(String, boolean, TimeScale)
184      */
185     @DefaultDataContext
186     public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
187         this(supportedNames, missingCoefficientsAllowed,
188                 DataContext.getDefault().getTimeScales().getTT());
189     }
190 
191     /**
192      * Simple constructor.
193      *
194      * @param supportedNames             regular expression for supported files names
195      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input
196      *                                   data
197      * @param timeScale                  to use when parsing dates.
198      * @since 10.1
199      */
200     public ICGEMFormatReader(final String supportedNames,
201                              final boolean missingCoefficientsAllowed,
202                              final TimeScale timeScale) {
203         super(supportedNames, missingCoefficientsAllowed, timeScale);
204         referenceDate = null;
205         cTrend = new ArrayList<>();
206         sTrend = new ArrayList<>();
207         cCos   = new HashMap<>();
208         cSin   = new HashMap<>();
209         sCos   = new HashMap<>();
210         sSin   = new HashMap<>();
211     }
212 
213     /** {@inheritDoc} */
214     public void loadData(final InputStream input, final String name)
215         throws IOException, ParseException, OrekitException {
216 
217         // reset the indicator before loading any data
218         setReadComplete(false);
219         referenceDate = null;
220         cTrend.clear();
221         sTrend.clear();
222         cCos.clear();
223         cSin.clear();
224         sCos.clear();
225         sSin.clear();
226 
227         // by default, the field is normalized with unknown tide system
228         // (will be overridden later if non-default)
229         normalized = true;
230         TideSystem tideSystem = TideSystem.UNKNOWN;
231 
232         boolean inHeader = true;
233         double[][] c     = null;
234         double[][] s     = null;
235         boolean okCoeffs = false;
236         int lineNumber   = 0;
237         String line      = null;
238         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
239             for (line = r.readLine(); line != null; line = r.readLine()) {
240                 ++lineNumber;
241                 line = line.trim();
242                 if (line.length() == 0) {
243                     continue;
244                 }
245                 final String[] tab = SEPARATOR.split(line);
246                 if (inHeader) {
247                     if (tab.length == 2 && PRODUCT_TYPE.equals(tab[0])) {
248                         if (!GRAVITY_FIELD.equals(tab[1])) {
249                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
250                                                            lineNumber, name, line);
251                         }
252                     } else if (tab.length == 2 && tab[0].endsWith(GRAVITY_CONSTANT)) {
253                         setMu(parseDouble(tab[1]));
254                     } else if (tab.length == 2 && REFERENCE_RADIUS.equals(tab[0])) {
255                         setAe(parseDouble(tab[1]));
256                     } else if (tab.length == 2 && MAX_DEGREE.equals(tab[0])) {
257 
258                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
259                         final int order  = FastMath.min(getMaxParseOrder(), degree);
260                         c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
261                         s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
262 
263                     } else if (tab.length == 2 && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
264                         if (ZERO_TIDE.equals(tab[1])) {
265                             tideSystem = TideSystem.ZERO_TIDE;
266                         } else if (TIDE_FREE.equals(tab[1])) {
267                             tideSystem = TideSystem.TIDE_FREE;
268                         } else if (TIDE_UNKNOWN.equals(tab[1])) {
269                             tideSystem = TideSystem.UNKNOWN;
270                         } else {
271                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
272                                                            lineNumber, name, line);
273                         }
274                     } else if (tab.length == 2 && NORMALIZATION_INDICATOR.equals(tab[0])) {
275                         if (NORMALIZED.equals(tab[1])) {
276                             normalized = true;
277                         } else if (UNNORMALIZED.equals(tab[1])) {
278                             normalized = false;
279                         } else {
280                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
281                                                            lineNumber, name, line);
282                         }
283                     } else if (tab.length == 2 && END_OF_HEADER.equals(tab[0])) {
284                         inHeader = false;
285                     }
286                 } else {
287                     if (tab.length == 7 && GFC.equals(tab[0]) || tab.length == 8 && GFCT.equals(tab[0])) {
288 
289                         final int i = Integer.parseInt(tab[1]);
290                         final int j = Integer.parseInt(tab[2]);
291                         if (i < c.length && j < c[i].length) {
292 
293                             parseCoefficient(tab[3], c, i, j, "C", name);
294                             parseCoefficient(tab[4], s, i, j, "S", name);
295                             okCoeffs = true;
296 
297                             if (tab.length == 8) {
298                                 // check the reference date (format yyyymmdd)
299                                 final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
300                                                                                    Integer.parseInt(tab[7].substring(4, 6)),
301                                                                                    Integer.parseInt(tab[7].substring(6, 8)));
302                                 if (referenceDate == null) {
303                                     // first reference found, store it
304                                     referenceDate = toDate(localRef);
305                                 } else if (!referenceDate.equals(toDate(localRef))) {
306                                     final AbsoluteDate localDate = toDate(localRef);
307                                     throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
308                                                               referenceDate, localDate, name,
309                                                               localDate.durationFrom(referenceDate));
310                                 }
311                             }
312 
313                         }
314                     } else if (tab.length == 7 && (DOT.equals(tab[0]) || TRND.equals(tab[0]))) {
315 
316                         final int i = Integer.parseInt(tab[1]);
317                         final int j = Integer.parseInt(tab[2]);
318                         if (i < c.length && j < c[i].length) {
319 
320                             // store the secular trend coefficients
321                             extendListOfLists(cTrend, i, j, 0.0);
322                             extendListOfLists(sTrend, i, j, 0.0);
323                             parseCoefficient(tab[3], cTrend, i, j, "Ctrend", name);
324                             parseCoefficient(tab[4], sTrend, i, j, "Strend", name);
325 
326                         }
327 
328                     } else if (tab.length == 8 && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
329 
330                         final int i = Integer.parseInt(tab[1]);
331                         final int j = Integer.parseInt(tab[2]);
332                         if (i < c.length && j < c[i].length) {
333 
334                             // extract arrays corresponding to period
335                             final Double period = Double.valueOf(tab[7]);
336                             if (!cCos.containsKey(period)) {
337                                 cCos.put(period, new ArrayList<>());
338                                 cSin.put(period, new ArrayList<>());
339                                 sCos.put(period, new ArrayList<>());
340                                 sSin.put(period, new ArrayList<>());
341                             }
342                             final List<List<Double>> cCosPeriod = cCos.get(period);
343                             final List<List<Double>> cSinPeriod = cSin.get(period);
344                             final List<List<Double>> sCosPeriod = sCos.get(period);
345                             final List<List<Double>> sSinPeriod = sSin.get(period);
346 
347                             // store the pulsation coefficients
348                             extendListOfLists(cCosPeriod, i, j, 0.0);
349                             extendListOfLists(cSinPeriod, i, j, 0.0);
350                             extendListOfLists(sCosPeriod, i, j, 0.0);
351                             extendListOfLists(sSinPeriod, i, j, 0.0);
352                             if (ACOS.equals(tab[0])) {
353                                 parseCoefficient(tab[3], cCosPeriod, i, j, "Ccos", name);
354                                 parseCoefficient(tab[4], sCosPeriod, i, j, "SCos", name);
355                             } else {
356                                 parseCoefficient(tab[3], cSinPeriod, i, j, "Csin", name);
357                                 parseCoefficient(tab[4], sSinPeriod, i, j, "Ssin", name);
358                             }
359 
360                         }
361 
362                     } else {
363                         throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
364                                                        lineNumber, name, line);
365                     }
366                 }
367 
368             }
369         } catch (NumberFormatException nfe) {
370             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
371                                       lineNumber, name, line);
372         }
373 
374 
375         if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
376             // ensure at least the (0, 0) element is properly set
377             if (Precision.equals(c[0][0], 0.0, 0)) {
378                 c[0][0] = 1.0;
379             }
380         }
381 
382         if (Double.isNaN(getAe()) || Double.isNaN(getMu()) || !okCoeffs) {
383             String loaderName = getClass().getName();
384             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
385             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
386                                       name, loaderName);
387         }
388 
389         setRawCoefficients(normalized, c, s, name);
390         setTideSystem(tideSystem);
391         setReadComplete(true);
392 
393     }
394 
395     /** Get a provider for read spherical harmonics coefficients.
396      * <p>
397      * ICGEM fields do include time-dependent parts which are taken into account
398      * in the returned provider.
399      * </p>
400      * @param wantNormalized if true, the provider will provide normalized coefficients,
401      * otherwise it will provide un-normalized coefficients
402      * @param degree maximal degree
403      * @param order maximal order
404      * @return a new provider
405      * @since 6.0
406      */
407     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
408                                                      final int degree, final int order) {
409 
410         RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
411         if (cTrend.isEmpty() && cCos.isEmpty()) {
412             // there are no time-dependent coefficients
413             return provider;
414         }
415 
416         if (!cTrend.isEmpty()) {
417 
418             // add the secular trend layer
419             final double[][] cArrayTrend = toArray(cTrend);
420             final double[][] sArrayTrend = toArray(sTrend);
421             rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
422             provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);
423 
424         }
425 
426         for (final Map.Entry<Double, List<List<Double>>> entry : cCos.entrySet()) {
427 
428             final double period = entry.getKey();
429 
430             // add the pulsating layer for the current period
431             final double[][] cArrayCos = toArray(cCos.get(period));
432             final double[][] sArrayCos = toArray(sCos.get(period));
433             final double[][] cArraySin = toArray(cSin.get(period));
434             final double[][] sArraySin = toArray(sSin.get(period));
435             rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
436             rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
437             provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
438                                                        cArrayCos, cArraySin, sArrayCos, sArraySin);
439 
440         }
441 
442         return provider;
443 
444     }
445 
446 }