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