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.List;
27 import java.util.regex.Matcher;
28 import java.util.regex.Pattern;
29
30 import org.hipparchus.util.FastMath;
31 import org.hipparchus.util.Precision;
32 import org.orekit.annotation.DefaultDataContext;
33 import org.orekit.data.DataContext;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.errors.OrekitMessages;
36 import org.orekit.errors.OrekitParseException;
37 import org.orekit.time.AbsoluteDate;
38 import org.orekit.time.DateComponents;
39 import org.orekit.time.TimeScale;
40 import org.orekit.utils.Constants;
41
42 /** Reader for the GRGS gravity field format.
43 *
44 * <p> This format was used to describe various gravity fields at GRGS (Toulouse).
45 *
46 * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
47 * which will determine which reader to use with the selected gravity field file.</p>
48 *
49 * @see GravityFields
50 * @author Luc Maisonobe
51 */
52 public class GRGSFormatReader extends PotentialCoefficientsReader {
53
54 /** Patterns for lines (the last pattern is repeated for all data lines). */
55 private static final Pattern[] LINES;
56
57 /** Reference date. */
58 private AbsoluteDate referenceDate;
59
60 /** Secular drift of the cosine coefficients. */
61 private final List<List<Double>> cDot;
62
63 /** Secular drift of the sine coefficients. */
64 private final List<List<Double>> sDot;
65
66 static {
67
68 // sub-patterns
69 final String real = "[-+]?\\d?\\.\\d+[eEdD][-+]\\d\\d";
70 final String sep = ")\\s*(";
71
72 // regular expression for header lines
73 final String[] header = {
74 "^\\s*FIELD - .*$",
75 "^\\s+AE\\s+1/F\\s+GM\\s+OMEGA\\s*$",
76 "^\\s*(" + real + sep + real + sep + real + sep + real + ")\\s*$",
77 "^\\s*REFERENCE\\s+DATE\\s+:\\s+(\\d+)\\.0+\\s*$",
78 "^\\s*MAXIMAL\\s+DEGREE\\s+:\\s+(\\d+)\\s.*$",
79 "^\\s*L\\s+M\\s+DOT\\s+CBAR\\s+SBAR\\s+SIGMA C\\s+SIGMA S(\\s+LIB)?\\s*$"
80 };
81
82 // regular expression for data lines
83 final String data = "^([ 0-9]{3})([ 0-9]{3})( |DOT)\\s*(" +
84 real + sep + real + sep + real + sep + real +
85 ")(\\s+[0-9]+)?\\s*$";
86
87 // compile the regular expressions
88 LINES = new Pattern[header.length + 1];
89 for (int i = 0; i < header.length; ++i) {
90 LINES[i] = Pattern.compile(header[i]);
91 }
92 LINES[LINES.length - 1] = Pattern.compile(data);
93
94 }
95
96 /** Simple constructor.
97 *
98 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
99 *
100 * @param supportedNames regular expression for supported files names
101 * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
102 * @see #GRGSFormatReader(String, boolean, TimeScale)
103 */
104 @DefaultDataContext
105 public GRGSFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
106 this(supportedNames, missingCoefficientsAllowed,
107 DataContext.getDefault().getTimeScales().getTT());
108 }
109
110 /**
111 * Simple constructor.
112 *
113 * @param supportedNames regular expression for supported files names
114 * @param missingCoefficientsAllowed if true, allows missing coefficients in the input
115 * data
116 * @param timeScale to use when parsing dates.
117 * @since 10.1
118 */
119 public GRGSFormatReader(final String supportedNames,
120 final boolean missingCoefficientsAllowed,
121 final TimeScale timeScale) {
122 super(supportedNames, missingCoefficientsAllowed, timeScale);
123 referenceDate = null;
124 cDot = new ArrayList<>();
125 sDot = new ArrayList<>();
126 }
127
128 /** {@inheritDoc} */
129 public void loadData(final InputStream input, final String name)
130 throws IOException, ParseException, OrekitException {
131
132 // reset the indicator before loading any data
133 setReadComplete(false);
134 referenceDate = null;
135 cDot.clear();
136 sDot.clear();
137
138 // FIELD - GRIM5, VERSION : C1, november 1999
139 // AE 1/F GM OMEGA
140 //0.63781364600000E+070.29825765000000E+030.39860044150000E+150.72921150000000E-04
141 //REFERENCE DATE : 1997.00
142 //MAXIMAL DEGREE : 120 Sigmas calibration factor : .5000E+01 (applied)
143 //L M DOT CBAR SBAR SIGMA C SIGMA S
144 // 2 0DOT 0.13637590952454E-10 0.00000000000000E+00 .143968E-11 .000000E+00
145 // 3 0DOT 0.28175700027753E-11 0.00000000000000E+00 .496704E-12 .000000E+00
146 // 4 0DOT 0.12249148508277E-10 0.00000000000000E+00 .129977E-11 .000000E+00
147 // 0 0 .99999999988600E+00 .00000000000000E+00 .153900E-09 .000000E+00
148 // 2 0 -0.48416511550920E-03 0.00000000000000E+00 .204904E-10 .000000E+00
149
150 int lineNumber = 0;
151 double[][] c = null;
152 double[][] s = null;
153 try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
154 for (String line = r.readLine(); line != null; line = r.readLine()) {
155
156 ++lineNumber;
157
158 // match current header or data line
159 final Matcher matcher = LINES[FastMath.min(LINES.length, lineNumber) - 1].matcher(line);
160 if (!matcher.matches()) {
161 throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
162 lineNumber, name, line);
163 }
164
165 if (lineNumber == 3) {
166 // header line defining ae, 1/f, GM and Omega
167 setAe(parseDouble(matcher.group(1)));
168 setMu(parseDouble(matcher.group(3)));
169 } else if (lineNumber == 4) {
170 // header line containing the reference date
171 referenceDate = toDate(
172 new DateComponents(Integer.parseInt(matcher.group(1)), 1, 1));
173 } else if (lineNumber == 5) {
174 // header line defining max degree
175 final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(matcher.group(1)));
176 final int order = FastMath.min(getMaxParseOrder(), degree);
177 c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
178 s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
179 } else if (lineNumber > 6) {
180 // data line
181 final int i = Integer.parseInt(matcher.group(1).trim());
182 final int j = Integer.parseInt(matcher.group(2).trim());
183 if (i < c.length && j < c[i].length) {
184 if ("DOT".equals(matcher.group(3).trim())) {
185
186 // store the secular drift coefficients
187 extendListOfLists(cDot, i, j, 0.0);
188 extendListOfLists(sDot, i, j, 0.0);
189 parseCoefficient(matcher.group(4), cDot, i, j, "Cdot", name);
190 parseCoefficient(matcher.group(5), sDot, i, j, "Sdot", name);
191
192 } else {
193
194 // store the constant coefficients
195 parseCoefficient(matcher.group(4), c, i, j, "C", name);
196 parseCoefficient(matcher.group(5), s, i, j, "S", name);
197
198 }
199 }
200 }
201
202 }
203 }
204
205 if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
206 // ensure at least the (0, 0) element is properly set
207 if (Precision.equals(c[0][0], 0.0, 0)) {
208 c[0][0] = 1.0;
209 }
210 }
211
212 setRawCoefficients(true, c, s, name);
213 setTideSystem(TideSystem.UNKNOWN);
214 setReadComplete(true);
215
216 }
217
218 /** Get a provider for read spherical harmonics coefficients.
219 * <p>
220 * GRGS fields may include time-dependent parts which are taken into account
221 * in the returned provider.
222 * </p>
223 * @param wantNormalized if true, the provider will provide normalized coefficients,
224 * otherwise it will provide un-normalized coefficients
225 * @param degree maximal degree
226 * @param order maximal order
227 * @return a new provider
228 * @since 6.0
229 */
230 public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
231 final int degree, final int order) {
232
233 // get the constant part
234 RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
235
236 if (!cDot.isEmpty()) {
237
238 // add the secular trend layer
239 final double[][] cArray = toArray(cDot);
240 final double[][] sArray = toArray(sDot);
241 rescale(1.0 / Constants.JULIAN_YEAR, true, cArray, sArray, wantNormalized, cArray, sArray);
242 provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArray, sArray);
243
244 }
245
246 return provider;
247
248 }
249
250 }