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.List;
26  import java.util.Locale;
27  
28  import org.hipparchus.util.FastMath;
29  import org.hipparchus.util.Precision;
30  import org.orekit.errors.OrekitException;
31  import org.orekit.errors.OrekitMessages;
32  import org.orekit.time.DateComponents;
33  import org.orekit.utils.Constants;
34  
35  /** Reader for the SHM gravity field format.
36   *
37   * <p> This format was used to describe the gravity field of EIGEN models
38   * published by the GFZ Potsdam up to 2003. It was then replaced by
39   * {@link ICGEMFormatReader ICGEM format}. The SHM format is described in
40   * <a href="http://op.gfz-potsdam.de/champ/docs_CHAMP/CH-FORMAT-REFLINKS.html"> Potsdam university
41   * website</a>.
42   *
43   * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
44   *  which will determine which reader to use with the selected gravity field file.</p>
45   *
46   * @see GravityFieldFactory
47   * @author Fabien Maussion
48   */
49  public class SHMFormatReader extends PotentialCoefficientsReader {
50  
51      /** First field labels. */
52      private static final String GRCOEF = "GRCOEF";
53  
54      /** Second field labels. */
55      private static final String GRCOF2 = "GRCOF2";
56  
57      /** Drift coefficients labels. */
58      private static final String GRDOTA = "GRDOTA";
59  
60      /** Reference date. */
61      private DateComponents referenceDate;
62  
63      /** Secular drift of the cosine coefficients. */
64      private final List<List<Double>> cDot;
65  
66      /** Secular drift of the sine coefficients. */
67      private final List<List<Double>> sDot;
68  
69      /** Simple constructor.
70       * @param supportedNames regular expression for supported files names
71       * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
72       */
73      public SHMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
74          super(supportedNames, missingCoefficientsAllowed);
75          referenceDate = null;
76          cDot = new ArrayList<List<Double>>();
77          sDot = new ArrayList<List<Double>>();
78      }
79  
80      /** {@inheritDoc} */
81      public void loadData(final InputStream input, final String name)
82          throws IOException, ParseException, OrekitException {
83  
84          // reset the indicator before loading any data
85          setReadComplete(false);
86          referenceDate = null;
87          cDot.clear();
88          sDot.clear();
89  
90          boolean    normalized = false;
91          TideSystem tideSystem = TideSystem.UNKNOWN;
92  
93          final BufferedReader r = new BufferedReader(new InputStreamReader(input, "UTF-8"));
94          boolean okEarth            = false;
95          boolean okSHM              = false;
96          boolean okCoeffs           = false;
97          double[][] c               = null;
98          double[][] s               = null;
99          String line = r.readLine();
100         if ((line != null) &&
101             "FIRST ".equals(line.substring(0, 6)) &&
102             "SHM    ".equals(line.substring(49, 56))) {
103             for (line = r.readLine(); line != null; line = r.readLine()) {
104                 if (line.length() >= 6) {
105                     final String[] tab = line.split("\\s+");
106 
107                     // read the earth values
108                     if ("EARTH".equals(tab[0])) {
109                         setMu(parseDouble(tab[1]));
110                         setAe(parseDouble(tab[2]));
111                         okEarth = true;
112                     }
113 
114                     // initialize the arrays
115                     if ("SHM".equals(tab[0])) {
116 
117                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
118                         final int order  = FastMath.min(getMaxParseOrder(), degree);
119                         c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
120                         s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
121                         final String lowerCaseLine = line.toLowerCase(Locale.US);
122                         normalized = lowerCaseLine.contains("fully normalized");
123                         if (lowerCaseLine.contains("exclusive permanent tide")) {
124                             tideSystem = TideSystem.TIDE_FREE;
125                         } else {
126                             tideSystem = TideSystem.UNKNOWN;
127                         }
128                         okSHM = true;
129                     }
130 
131                     // fill the arrays
132                     if (GRCOEF.equals(line.substring(0, 6)) || GRCOF2.equals(tab[0]) || GRDOTA.equals(tab[0])) {
133                         final int i = Integer.parseInt(tab[1]);
134                         final int j = Integer.parseInt(tab[2]);
135                         if (i < c.length && j < c[i].length) {
136                             if (GRDOTA.equals(tab[0])) {
137 
138                                 // store the secular drift coefficients
139                                 extendListOfLists(cDot, i, j, 0.0);
140                                 extendListOfLists(sDot, i, j, 0.0);
141                                 parseCoefficient(tab[3], cDot, i, j, "Cdot", name);
142                                 parseCoefficient(tab[4], sDot, i, j, "Sdot", name);
143 
144                                 // check the reference date (format yyyymmdd)
145                                 final DateComponents#DateComponents">DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
146                                                                                    Integer.parseInt(tab[7].substring(4, 6)),
147                                                                                    Integer.parseInt(tab[7].substring(6, 8)));
148                                 if (referenceDate == null) {
149                                     // first reference found, store it
150                                     referenceDate = localRef;
151                                 } else if (!referenceDate.equals(localRef)) {
152                                     throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
153                                                               referenceDate, localRef, name);
154                                 }
155 
156                             } else {
157 
158                                 // store the constant coefficients
159                                 parseCoefficient(tab[3], c, i, j, "C", name);
160                                 parseCoefficient(tab[4], s, i, j, "S", name);
161                                 okCoeffs = true;
162 
163                             }
164                         }
165                     }
166 
167                 }
168             }
169         }
170 
171         if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
172             // ensure at least the (0, 0) element is properly set
173             if (Precision.equals(c[0][0], 0.0, 0)) {
174                 c[0][0] = 1.0;
175             }
176         }
177 
178         if (!(okEarth && okSHM && okCoeffs)) {
179             String loaderName = getClass().getName();
180             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
181             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
182                                       name, loaderName);
183         }
184 
185         setRawCoefficients(normalized, c, s, name);
186         setTideSystem(tideSystem);
187         setReadComplete(true);
188 
189     }
190 
191     /** Get a provider for read spherical harmonics coefficients.
192      * <p>
193      * SHM fields do include time-dependent parts which are taken into account
194      * in the returned provider.
195      * </p>
196      * @param wantNormalized if true, the provider will provide normalized coefficients,
197      * otherwise it will provide un-normalized coefficients
198      * @param degree maximal degree
199      * @param order maximal order
200      * @return a new provider
201      * @since 6.0
202      */
203     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
204                                                      final int degree, final int order) {
205 
206         // get the constant part
207         RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
208 
209         if (!cDot.isEmpty()) {
210 
211             // add the secular trend layer
212             final double[][] cArray = toArray(cDot);
213             final double[][] sArray = toArray(sDot);
214             rescale(1.0 / Constants.JULIAN_YEAR, true, cArray, sArray, wantNormalized, cArray, sArray);
215             provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArray, sArray);
216 
217         }
218 
219         return provider;
220 
221     }
222 
223 }