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.data;
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.util.Arrays;
25 import java.util.HashMap;
26 import java.util.Map;
27 import java.util.regex.Matcher;
28 import java.util.regex.Pattern;
29
30 import org.hipparchus.exception.DummyLocalizable;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.Precision;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitMessages;
35
36 /**
37 * Parser for {@link PoissonSeries Poisson series} files.
38 * <p>
39 * A Poisson series is composed of a time polynomial part and a non-polynomial
40 * part which consist in summation series. The {@link SeriesTerm series terms}
41 * are harmonic functions (combination of sines and cosines) of polynomial
42 * <em>arguments</em>. The polynomial arguments are combinations of luni-solar or
43 * planetary {@link BodiesElements elements}.
44 * </p>
45 * <p>
46 * The Poisson series files from IERS have various formats, with or without
47 * polynomial part, with or without planetary components, with or without
48 * period column, with terms of increasing degrees either in dedicated columns
49 * or in successive sections of the file ... This class attempts to read all the
50 * commonly found formats, by specifying the columns of interest.
51 * </p>
52 * <p>
53 * The handling of increasing degrees terms (i.e. sin, cos, t sin, t cos, t^2 sin,
54 * t^2 cos ...) is done as follows.
55 * </p>
56 * <ul>
57 * <li>user must specify pairs of columns to be extracted at each line,
58 * in increasing degree order</li>
59 * <li>negative columns indices correspond to inexistent values that will be
60 * replaced by 0.0)</li>
61 * <li>file may provide section headers to specify a degree, which is added
62 * to the current column degree</li>
63 * </ul>
64 * <p>
65 * A file from an old convention, like table 5.1 in IERS conventions 1996, uses
66 * separate columns for degree 0 and degree 1, and uses only sine for nutation in
67 * longitude and cosine for nutation in obliquity. It reads as follows:
68 * </p>
69 * <pre>
70 * ∆ψ = Σ (Ai+A'it) sin(ARGUMENT), ∆ε = Σ (Bi+B'it) cos(ARGUMENT)
71 *
72 * MULTIPLIERS OF PERIOD LONGITUDE OBLIQUITY
73 * l l' F D Om days Ai A'i Bi B'i
74 *
75 * 0 0 0 0 1 -6798.4 -171996 -174.2 92025 8.9
76 * 0 0 2 -2 2 182.6 -13187 -1.6 5736 -3.1
77 * 0 0 2 0 2 13.7 -2274 -0.2 977 -0.5
78 * 0 0 0 0 2 -3399.2 2062 0.2 -895 0.5
79 * </pre>
80 * <p>
81 * In order to parse the nutation in longitude from the previous table, the
82 * following settings should be used:
83 * </p>
84 * <ul>
85 * <li>totalColumns = 10 (see {@link #PoissonSeriesParser(int)})</li>
86 * <li>firstDelaunay = 1 (see {@link #withFirstDelaunay(int)})</li>
87 * <li>no calls to {@link #withFirstPlanetary(int)} as there are no planetary columns in this table</li>
88 * <li>sinCosColumns = 7, -1 for degree 0 for Ai (see {@link #withSinCos(int, int, double, int, double)})</li>
89 * <li>sinCosColumns = 8, -1 for degree 1 for A'i (see {@link #withSinCos(int, int, double, int, double)})</li>
90 * </ul>
91 * <p>
92 * In order to parse the nutation in obliquity from the previous table, the
93 * following settings should be used:
94 * </p>
95 * <ul>
96 * <li>totalColumns = 10 (see {@link #PoissonSeriesParser(int)})</li>
97 * <li>firstDelaunay = 1 (see {@link #withFirstDelaunay(int)})</li>
98 * <li>no calls to {@link #withFirstPlanetary(int)} as there are no planetary columns in this table</li>
99 * <li>sinCosColumns = -1, 9 for degree 0 for Bi (see {@link #withSinCos(int, int, double, int, double)})</li>
100 * <li>sinCosColumns = -1, 10 for degree 1 for B'i (see {@link #withSinCos(int, int, double, int, double)})</li>
101 * </ul>
102 * <p>
103 * A file from a recent convention, like table 5.3a in IERS conventions 2010, uses
104 * only two columns for sin and cos, and separate degrees in successive sections with
105 * dedicated headers. It reads as follows:
106 * </p>
107 * <pre>
108 * ---------------------------------------------------------------------------------------------------
109 *
110 * (unit microarcsecond; cut-off: 0.1 microarcsecond)
111 * (ARG being for various combination of the fundamental arguments of the nutation theory)
112 *
113 * Sum_i[A_i * sin(ARG) + A"_i * cos(ARG)]
114 *
115 * + Sum_i[A'_i * sin(ARG) + A"'_i * cos(ARG)] * t (see Chapter 5, Eq. (35))
116 *
117 * The Table below provides the values for A_i and A"_i (j=0) and then A'_i and A"'_i (j=1)
118 *
119 * The expressions for the fundamental arguments appearing in columns 4 to 8 (luni-solar part)
120 * and in columns 9 to 17 (planetary part) are those of the IERS Conventions 2003
121 *
122 * ----------------------------------------------------------------------------------------------------------
123 * j = 0 Number of terms = 1320
124 * ----------------------------------------------------------------------------------------------------------
125 * i A_i A"_i l l' F D Om L_Me L_Ve L_E L_Ma L_J L_Sa L_U L_Ne p_A
126 * ----------------------------------------------------------------------------------------------------------
127 * 1 -17206424.18 3338.60 0 0 0 0 1 0 0 0 0 0 0 0 0 0
128 * 2 -1317091.22 -1369.60 0 0 2 -2 2 0 0 0 0 0 0 0 0 0
129 * 3 -227641.81 279.60 0 0 2 0 2 0 0 0 0 0 0 0 0 0
130 * 4 207455.40 -69.80 0 0 0 0 2 0 0 0 0 0 0 0 0 0
131 * 5 147587.70 1181.70 0 1 0 0 0 0 0 0 0 0 0 0 0 0
132 *
133 * ...
134 *
135 * 1319 -0.10 0.00 0 0 0 0 0 1 0 -3 0 0 0 0 0 -2
136 * 1320 -0.10 0.00 0 0 0 0 0 0 0 1 0 1 -2 0 0 0
137 *
138 * --------------------------------------------------------------------------------------------------------------
139 * j = 1 Number of terms = 38
140 * --------------------------------------------------------------------------------------------------------------
141 * i A'_i A"'_i l l' F D Om L_Me L_Ve L_E L_Ma L_J L_Sa L_U L_Ne p_A
142 * --------------------------------------------------------------------------------------------------------------
143 * 1321 -17418.82 2.89 0 0 0 0 1 0 0 0 0 0 0 0 0 0
144 * 1322 -363.71 -1.50 0 1 0 0 0 0 0 0 0 0 0 0 0 0
145 * 1323 -163.84 1.20 0 0 2 -2 2 0 0 0 0 0 0 0 0 0
146 * 1324 122.74 0.20 0 1 2 -2 2 0 0 0 0 0 0 0 0 0
147 * </pre>
148 * <p>
149 * In order to parse the nutation in longitude from the previous table, the
150 * following settings should be used:
151 * </p>
152 * <ul>
153 * <li>totalColumns = 17 (see {@link #PoissonSeriesParser(int)})</li>
154 * <li>firstDelaunay = 4 (see {@link #withFirstDelaunay(int)})</li>
155 * <li>firstPlanetary = 9 (see {@link #withFirstPlanetary(int)})</li>
156 * <li>sinCosColumns = 2,3 (we specify only degree 0, so when we read
157 * section j = 0 we read degree 0, when we read section j = 1 we read
158 * degree 1, see {@link #withSinCos(int, int, double, int, double)} ...)</li>
159 * </ul>
160 * <p>
161 * A file from a recent convention, like table 6.5a in IERS conventions 2010, contains
162 * both Doodson arguments (τ, s, h, p, N', ps), Doodson numbers and Delaunay parameters.
163 * In this case, the coefficients for the Delaunay parameters must be <em>subtracted</em>
164 * from the τ = GMST + π tide parameter, so the signs in the files must be reversed
165 * in order to match the Doodson arguments and Doodson numbers. This is done automatically
166 * (and consistency is checked) only when the {@link #withDoodson(int, int)} method is
167 * called at parser configuration time. Some other files use the γ = GMST + π tide parameter
168 * rather than Doodson τ argument and the coefficients for the Delaunay parameters must be
169 * <em>added</em> to the γ parameter, so no sign reversal is performed. In order to avoid
170 * ambiguity as the two cases are incompatible with each other, trying to add a configuration
171 * for τ by calling {@link #withDoodson(int, int)} and to also add a configuration for γ by
172 * calling {@link #withGamma(int)} triggers an exception.
173 * </p>
174 * <p>The table 6.5a file also contains a column for the waves names (the Darwin's symbol)
175 * which may be empty, so it must be identified explicitly by calling {@link
176 * #withOptionalColumn(int)}. The 6.5a table reads as follows:
177 * </p>
178 * <pre>
179 * The in-phase (ip) amplitudes (A₁ δkfR Hf) and the out-of-phase (op) amplitudes (A₁ δkfI Hf)
180 * of the corrections for frequency dependence of k₂₁⁽⁰⁾, taking the nominal value k₂₁ for the
181 * diurnal tides as (0.29830 − i 0.00144). Units: 10⁻¹² . The entries for δkfR and δkfI are in
182 * units of 10⁻⁵. Multipliers of the Doodson arguments identifying the tidal terms are given,
183 * as also those of the Delaunay variables characterizing the nutations produced by these
184 * terms.
185 *
186 * Name deg/hr Doodson τ s h p N' ps l l' F D Ω δkfR δkfI Amp. Amp.
187 * No. /10−5 /10−5 (ip) (op)
188 * 2Q₁ 12.85429 125,755 1 -3 0 2 0 0 2 0 2 0 2 -29 3 -0.1 0.0
189 * σ₁ 12.92714 127,555 1 -3 2 0 0 0 0 0 2 2 2 -30 3 -0.1 0.0
190 * 13.39645 135,645 1 -2 0 1 -1 0 1 0 2 0 1 -45 5 -0.1 0.0
191 * Q₁ 13.39866 135,655 1 -2 0 1 0 0 1 0 2 0 2 -46 5 -0.7 0.1
192 * ρ₁ 13.47151 137,455 1 -2 2 -1 0 0 -1 0 2 2 2 -49 5 -0.1 0.0
193 * </pre>
194 * <ul>
195 * <li>totalColumns = 18 (see {@link #PoissonSeriesParser(int)})</li>
196 * <li>optionalColumn = 1 (see {@link #withOptionalColumn(int)})</li>
197 * <li>firstDoodson, Doodson number = 4, 3 (see {@link #withDoodson(int, int)})</li>
198 * <li>firstDelaunay = 10 (see {@link #withFirstDelaunay(int)})</li>
199 * <li>sinCosColumns = 17, 18, see {@link #withSinCos(int, int, double, int, double)} ...)</li>
200 * </ul>
201 * <p>
202 * Our parsing algorithm involves adding the section degree from the "j = 0, 1, 2 ..." header
203 * to the column degree. A side effect of this algorithm is that it is theoretically possible
204 * to mix both formats and have for example degree two term appear as degree 2 column in section
205 * j=0 and as degree 1 column in section j=1 and as degree 0 column in section j=2. This case
206 * is not expected to be encountered in practice. The real files use either several columns
207 * <em>or</em> several sections, but not both at the same time.
208 * </p>
209 *
210 * @author Luc Maisonobe
211 * @see SeriesTerm
212 * @see PolynomialNutation
213 * @since 6.1
214 */
215 public class PoissonSeriesParser {
216
217 /** Default pattern for fields with unknown type (non-space characters). */
218 private static final String UNKNOWN_TYPE_PATTERN = "\\S+";
219
220 /** Pattern for optional fields (either nothing or non-space characters). */
221 private static final String OPTIONAL_FIELD_PATTERN = "\\S*";
222
223 /** Pattern for fields with integer type. */
224 private static final String INTEGER_TYPE_PATTERN = "[-+]?\\p{Digit}+";
225
226 /** Pattern for fields with real type. */
227 private static final String REAL_TYPE_PATTERN = "[-+]?(?:(?:\\p{Digit}+(?:\\.\\p{Digit}*)?)|(?:\\.\\p{Digit}+))(?:[eE][-+]?\\p{Digit}+)?";
228
229 /** Pattern for fields with Doodson number. */
230 private static final String DOODSON_TYPE_PATTERN = "\\p{Digit}{2,3}[.,]\\p{Digit}{3}";
231
232 /** Pattern for String replacement. */
233 private static final Pattern PATTERN = Pattern.compile("[.,]");
234
235 /** Parser for the polynomial part. */
236 private final PolynomialParser polynomialParser;
237
238 /** Fields patterns. */
239 private final String[] fieldsPatterns;
240
241 /** Optional column (counting from 1). */
242 private final int optional;
243
244 /** Column of the γ = GMST + π tide multiplier (counting from 1). */
245 private final int gamma;
246
247 /** Column of the first Doodson multiplier (counting from 1). */
248 private final int firstDoodson;
249
250 /** Column of the Doodson number (counting from 1). */
251 private final int doodson;
252
253 /** Column of the first Delaunay multiplier (counting from 1). */
254 private final int firstDelaunay;
255
256 /** Column of the first planetary multiplier (counting from 1). */
257 private final int firstPlanetary;
258
259 /** columns of the sine and cosine coefficients for successive degrees.
260 * <p>
261 * The ordering is: sin, cos, t sin, t cos, t^2 sin, t^2 cos ...
262 * </p>
263 */
264 private final int[] sinCosColumns;
265
266 /** Multiplicative factors to use for various columns. */
267 private final double[] sinCosFactors;
268
269 /** Build a parser for a Poisson series from an IERS table file.
270 * @param polynomialParser polynomial parser to use
271 * @param fieldsPatterns patterns for fields
272 * @param optional optional column
273 * @param gamma column of the GMST tide multiplier
274 * @param firstDoodson column of the first Doodson multiplier
275 * @param doodson column of the Doodson number
276 * @param firstDelaunay column of the first Delaunay multiplier
277 * @param firstPlanetary column of the first planetary multiplier
278 * @param sinCosColumns columns of the sine and cosine coefficients
279 * @param factors multiplicative factors to use for various columns
280 */
281 private PoissonSeriesParser(final PolynomialParser polynomialParser,
282 final String[] fieldsPatterns, final int optional,
283 final int gamma, final int firstDoodson,
284 final int doodson, final int firstDelaunay,
285 final int firstPlanetary, final int[] sinCosColumns,
286 final double[] factors) {
287 this.polynomialParser = polynomialParser;
288 this.fieldsPatterns = fieldsPatterns;
289 this.optional = optional;
290 this.gamma = gamma;
291 this.firstDoodson = firstDoodson;
292 this.doodson = doodson;
293 this.firstDelaunay = firstDelaunay;
294 this.firstPlanetary = firstPlanetary;
295 this.sinCosColumns = sinCosColumns;
296 this.sinCosFactors = factors;
297 }
298
299 /** Build a parser for a Poisson series from an IERS table file.
300 * @param totalColumns total number of columns in the non-polynomial sections
301 */
302 public PoissonSeriesParser(final int totalColumns) {
303 this(null, createInitialFieldsPattern(totalColumns), -1,
304 -1, -1, -1, -1, -1, new int[0], new double[0]);
305 }
306
307 /** Create an array with only non-space fields patterns.
308 * @param totalColumns total number of columns
309 * @return a new fields pattern array
310 */
311 private static String[] createInitialFieldsPattern(final int totalColumns) {
312 final String[] patterns = new String[totalColumns];
313 setPatterns(patterns, 1, totalColumns, UNKNOWN_TYPE_PATTERN);
314 return patterns;
315 }
316
317 /** Set fields patterns.
318 * @param array fields pattern array to modify
319 * @param first first column to set (counting from 1), do nothing if non-positive
320 * @param count number of columns to set
321 * @param pattern pattern to use
322 */
323 private static void setPatterns(final String[] array, final int first, final int count,
324 final String pattern) {
325 if (first > 0) {
326 Arrays.fill(array, first - 1, first - 1 + count, pattern);
327 }
328 }
329
330 /** Set up polynomial part parsing.
331 * @param freeVariable name of the free variable in the polynomial part
332 * @param unit default unit for polynomial, if not explicit within the file
333 * @return a new parser, with polynomial parser updated
334 */
335 public PoissonSeriesParser withPolynomialPart(final char freeVariable, final PolynomialParser.Unit unit) {
336 return new PoissonSeriesParser(new PolynomialParser(freeVariable, unit), fieldsPatterns, optional,
337 gamma, firstDoodson, doodson, firstDelaunay,
338 firstPlanetary, sinCosColumns, sinCosFactors);
339 }
340
341 /** Set up optional column.
342 * <p>
343 * Optional columns typically appears in tides-related files, as some waves have
344 * specific names (χ₁, M₂, ...) and other waves don't have names and hence are
345 * replaced by spaces in the corresponding file line.
346 * </p>
347 * <p>
348 * At most one column may be optional.
349 * </p>
350 * @param column optional column (counting from 1)
351 * @return a new parser, with updated columns settings
352 */
353 public PoissonSeriesParser withOptionalColumn(final int column) {
354
355 // update the fields pattern to expect 1 optional field at the right index
356 final String[] newFieldsPatterns = fieldsPatterns.clone();
357 setPatterns(newFieldsPatterns, optional, 1, UNKNOWN_TYPE_PATTERN);
358 setPatterns(newFieldsPatterns, column, 1, OPTIONAL_FIELD_PATTERN);
359
360 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, column,
361 gamma, firstDoodson, doodson, firstDelaunay,
362 firstPlanetary, sinCosColumns, sinCosFactors);
363
364 }
365
366 /** Set up column of GMST tide multiplier.
367 * @param column column of the GMST tide multiplier (counting from 1)
368 * @return a new parser, with updated columns settings
369 * @see #withDoodson(int, int)
370 */
371 public PoissonSeriesParser withGamma(final int column) {
372
373 // check we don't try to have both τ and γ configured at the same time
374 if (firstDoodson > 0 && column > 0) {
375 throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
376 }
377
378 // update the fields pattern to expect 1 integer at the right index
379 final String[] newFieldsPatterns = fieldsPatterns.clone();
380 setPatterns(newFieldsPatterns, gamma, 1, UNKNOWN_TYPE_PATTERN);
381 setPatterns(newFieldsPatterns, column, 1, INTEGER_TYPE_PATTERN);
382
383 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
384 column, firstDoodson, doodson, firstDelaunay,
385 firstPlanetary, sinCosColumns, sinCosFactors);
386
387 }
388
389 /** Set up columns for Doodson multipliers and Doodson number.
390 * @param firstMultiplierColumn column of the first Doodson multiplier which
391 * corresponds to τ (counting from 1)
392 * @param numberColumn column of the Doodson number (counting from 1)
393 * @return a new parser, with updated columns settings
394 * @see #withGamma(int)
395 * @see #withFirstDelaunay(int)
396 */
397 public PoissonSeriesParser withDoodson(final int firstMultiplierColumn, final int numberColumn) {
398
399 // check we don't try to have both τ and γ configured at the same time
400 if (gamma > 0 && firstMultiplierColumn > 0) {
401 throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
402 }
403
404 final String[] newFieldsPatterns = fieldsPatterns.clone();
405
406 // update the fields pattern to expect 6 integers at the right indices
407 setPatterns(newFieldsPatterns, firstDoodson, 6, UNKNOWN_TYPE_PATTERN);
408 setPatterns(newFieldsPatterns, firstMultiplierColumn, 6, INTEGER_TYPE_PATTERN);
409
410 // update the fields pattern to expect 1 Doodson number at the right index
411 setPatterns(newFieldsPatterns, doodson, 1, UNKNOWN_TYPE_PATTERN);
412 setPatterns(newFieldsPatterns, numberColumn, 1, DOODSON_TYPE_PATTERN);
413
414 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
415 gamma, firstMultiplierColumn, numberColumn, firstDelaunay,
416 firstPlanetary, sinCosColumns, sinCosFactors);
417
418 }
419
420 /** Set up first column of Delaunay multiplier.
421 * @param firstColumn column of the first Delaunay multiplier (counting from 1)
422 * @return a new parser, with updated columns settings
423 */
424 public PoissonSeriesParser withFirstDelaunay(final int firstColumn) {
425
426 // update the fields pattern to expect 5 integers at the right indices
427 final String[] newFieldsPatterns = fieldsPatterns.clone();
428 setPatterns(newFieldsPatterns, firstDelaunay, 5, UNKNOWN_TYPE_PATTERN);
429 setPatterns(newFieldsPatterns, firstColumn, 5, INTEGER_TYPE_PATTERN);
430
431 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
432 gamma, firstDoodson, doodson, firstColumn,
433 firstPlanetary, sinCosColumns, sinCosFactors);
434
435 }
436
437 /** Set up first column of planetary multiplier.
438 * @param firstColumn column of the first planetary multiplier (counting from 1)
439 * @return a new parser, with updated columns settings
440 */
441 public PoissonSeriesParser withFirstPlanetary(final int firstColumn) {
442
443 // update the fields pattern to expect 9 integers at the right indices
444 final String[] newFieldsPatterns = fieldsPatterns.clone();
445 setPatterns(newFieldsPatterns, firstPlanetary, 9, UNKNOWN_TYPE_PATTERN);
446 setPatterns(newFieldsPatterns, firstColumn, 9, INTEGER_TYPE_PATTERN);
447
448 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
449 gamma, firstDoodson, doodson, firstDelaunay,
450 firstColumn, sinCosColumns, sinCosFactors);
451
452 }
453
454 /** Set up columns of the sine and cosine coefficients.
455 * @param degree degree to set up
456 * @param sinColumn column of the sine coefficient for t<sup>degree</sup> counting from 1
457 * (may be -1 if there are no sine coefficients)
458 * @param sinFactor multiplicative factor for the sine coefficient
459 * @param cosColumn column of the cosine coefficient for t<sup>degree</sup> counting from 1
460 * (may be -1 if there are no cosine coefficients)
461 * @param cosFactor multiplicative factor for the cosine coefficient
462 * @return a new parser, with updated columns settings
463 */
464 public PoissonSeriesParser withSinCos(final int degree,
465 final int sinColumn, final double sinFactor,
466 final int cosColumn, final double cosFactor) {
467
468 // update the sin/cos columns array
469 final int maxDegree = FastMath.max(degree, sinCosColumns.length / 2 - 1);
470 final int[] newSinCosColumns = new int[2 * (maxDegree + 1)];
471 Arrays.fill(newSinCosColumns, -1);
472 System.arraycopy(sinCosColumns, 0, newSinCosColumns, 0, sinCosColumns.length);
473 newSinCosColumns[2 * degree] = sinColumn;
474 newSinCosColumns[2 * degree + 1] = cosColumn;
475
476 final double[] newSinCosFactors = new double[2 * (maxDegree + 1)];
477 Arrays.fill(newSinCosFactors, Double.NaN);
478 System.arraycopy(sinCosFactors, 0, newSinCosFactors, 0, sinCosFactors.length);
479 newSinCosFactors[2 * degree] = sinFactor;
480 newSinCosFactors[2 * degree + 1] = cosFactor;
481
482 // update the fields pattern to expect real numbers at the right indices
483 final String[] newFieldsPatterns = fieldsPatterns.clone();
484 if (2 * degree < sinCosColumns.length) {
485 setPatterns(newFieldsPatterns, sinCosColumns[2 * degree], 1, UNKNOWN_TYPE_PATTERN);
486 }
487 setPatterns(newFieldsPatterns, sinColumn, 1, REAL_TYPE_PATTERN);
488 if (2 * degree + 1 < sinCosColumns.length) {
489 setPatterns(newFieldsPatterns, sinCosColumns[2 * degree + 1], 1, UNKNOWN_TYPE_PATTERN);
490 }
491 setPatterns(newFieldsPatterns, cosColumn, 1, REAL_TYPE_PATTERN);
492
493 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
494 gamma, firstDoodson, doodson, firstDelaunay,
495 firstPlanetary, newSinCosColumns, newSinCosFactors);
496
497 }
498
499 /** Parse a stream.
500 * @param stream stream containing the IERS table
501 * @param name name of the resource file (for error messages only)
502 * @return parsed Poisson series
503 */
504 public PoissonSeries parse(final InputStream stream, final String name) {
505
506 if (stream == null) {
507 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
508 }
509
510 // the degrees section header should read something like:
511 // j = 0 Nb of terms = 1306
512 // or something like:
513 // j = 0 Number of terms = 1037
514 final Pattern degreeSectionHeaderPattern =
515 Pattern.compile("^\\p{Space}*j\\p{Space}*=\\p{Space}*(\\p{Digit}+)" +
516 "[\\p{Alpha}\\p{Space}]+=\\p{Space}*(\\p{Digit}+)\\p{Space}*$");
517
518 // regular lines are simply a space separated list of numbers
519 final StringBuilder builder = new StringBuilder("^\\p{Space}*");
520 for (int i = 0; i < fieldsPatterns.length; ++i) {
521 builder.append("(");
522 builder.append(fieldsPatterns[i]);
523 builder.append(")");
524 builder.append((i < fieldsPatterns.length - 1) ? "\\p{Space}+" : "\\p{Space}*$");
525 }
526 final Pattern regularLinePattern = Pattern.compile(builder.toString());
527
528 // setup the reader
529 try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {
530
531 int lineNumber = 0;
532 int expectedIndex = -1;
533 int nTerms = -1;
534 int count = 0;
535 int degree = 0;
536
537 // prepare the container for the parsed data
538 PolynomialNutation polynomial;
539 if (polynomialParser == null) {
540 // we don't expect any polynomial, we directly set the zero polynomial
541 polynomial = new PolynomialNutation(new double[0]);
542 } else {
543 // the dedicated parser will fill in the polynomial later
544 polynomial = null;
545 }
546 final Map<Long, SeriesTerm> series = new HashMap<Long, SeriesTerm>();
547
548 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
549
550 // replace unicode minus sign ('−') by regular hyphen ('-') for parsing
551 // such unicode characters occur in tables that are copy-pasted from PDF files
552 line = line.replace('\u2212', '-');
553 ++lineNumber;
554
555 final Matcher regularMatcher = regularLinePattern.matcher(line);
556 if (regularMatcher.matches()) {
557 // we have found a regular data line
558
559 if (expectedIndex > 0) {
560 // we are in a file were terms are numbered, we check the index
561 if (Integer.parseInt(regularMatcher.group(1)) != expectedIndex) {
562 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
563 lineNumber, name, regularMatcher.group());
564 }
565 }
566
567 // get the Doodson multipliers as well as the Doodson number
568 final int cTau = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson));
569 final int cS = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 1));
570 final int cH = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 2));
571 final int cP = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 3));
572 final int cNprime = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 4));
573 final int cPs = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 5));
574 final int nDoodson = (doodson < 0) ? 0 : Integer.parseInt(PATTERN.matcher(regularMatcher.group(doodson)).replaceAll(""));
575
576 // get the tide multiplier
577 int cGamma = (gamma < 0) ? 0 : Integer.parseInt(regularMatcher.group(gamma));
578
579 // get the Delaunay multipliers
580 int cL = Integer.parseInt(regularMatcher.group(firstDelaunay));
581 int cLPrime = Integer.parseInt(regularMatcher.group(firstDelaunay + 1));
582 int cF = Integer.parseInt(regularMatcher.group(firstDelaunay + 2));
583 int cD = Integer.parseInt(regularMatcher.group(firstDelaunay + 3));
584 int cOmega = Integer.parseInt(regularMatcher.group(firstDelaunay + 4));
585
586 // get the planetary multipliers
587 final int cMe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary));
588 final int cVe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 1));
589 final int cE = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 2));
590 final int cMa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 3));
591 final int cJu = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 4));
592 final int cSa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 5));
593 final int cUr = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 6));
594 final int cNe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 7));
595 final int cPa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 8));
596
597 if (nDoodson > 0) {
598
599 // set up the traditional parameters corresponding to the Doodson arguments
600 cGamma = cTau;
601 cL = -cL;
602 cLPrime = -cLPrime;
603 cF = -cF;
604 cD = -cD;
605 cOmega = -cOmega;
606
607 // check Doodson number, Doodson multipliers and Delaunay multipliers consistency
608 if (nDoodson != doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs) ||
609 nDoodson != delaunayToDoodsonNumber(cGamma, cL, cLPrime, cF, cD, cOmega)) {
610 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
611 lineNumber, name, regularMatcher.group());
612 }
613
614 }
615
616 final long key = NutationCodec.encode(cGamma, cL, cLPrime, cF, cD, cOmega,
617 cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
618
619 // retrieved the term, or build it if it's the first time it is encountered in the file
620 final SeriesTerm term;
621 if (series.containsKey(key)) {
622 // the term was already known, from another degree
623 term = series.get(key);
624 } else {
625 // the term is a new one
626 term = SeriesTerm.buildTerm(cGamma, cL, cLPrime, cF, cD, cOmega,
627 cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
628 }
629
630 boolean nonZero = false;
631 for (int d = 0; d < sinCosColumns.length / 2; ++d) {
632 final double sinCoeff =
633 parseCoefficient(regularMatcher, sinCosColumns[2 * d], sinCosFactors[2 * d]);
634 final double cosCoeff =
635 parseCoefficient(regularMatcher, sinCosColumns[2 * d + 1], sinCosFactors[2 * d + 1]);
636 if (!Precision.equals(sinCoeff, 0.0, 0) || !Precision.equals(cosCoeff, 0.0, 0)) {
637 nonZero = true;
638 term.add(0, degree + d, sinCoeff, cosCoeff);
639 ++count;
640 }
641 }
642 if (nonZero) {
643 series.put(key, term);
644 }
645
646 if (expectedIndex > 0) {
647 // we are in a file were terms are numbered
648 // we must update the expected value for next term
649 ++expectedIndex;
650 }
651
652 } else {
653
654 final Matcher headerMatcher = degreeSectionHeaderPattern.matcher(line);
655 if (headerMatcher.matches()) {
656
657 // we have found a degree section header
658 final int nextDegree = Integer.parseInt(headerMatcher.group(1));
659 if (nextDegree != degree + 1 && (degree != 0 || nextDegree != 0)) {
660 throw new OrekitException(OrekitMessages.MISSING_SERIE_J_IN_FILE,
661 degree + 1, name, lineNumber);
662 }
663
664 if (nextDegree == 0) {
665 // in IERS files split in sections, all terms are numbered
666 // we can check the indices
667 expectedIndex = 1;
668 }
669
670 if (nextDegree > 0 && count != nTerms) {
671 // the previous degree does not have the expected number of terms
672 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
673 }
674
675 // remember the number of terms the upcoming sublist should have
676 nTerms = Integer.parseInt(headerMatcher.group(2));
677 count = 0;
678 degree = nextDegree;
679
680 } else if (polynomial == null) {
681 // look for the polynomial part
682 final double[] coefficients = polynomialParser.parse(line);
683 if (coefficients != null) {
684 polynomial = new PolynomialNutation(coefficients);
685 }
686 }
687
688 }
689
690 }
691
692 if (polynomial == null || series.isEmpty()) {
693 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
694 }
695
696 if (nTerms > 0 && count != nTerms) {
697 // the last degree does not have the expected number of terms
698 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
699 }
700
701 // build the series
702 return new PoissonSeries(polynomial, series);
703
704 } catch (IOException ioe) {
705 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
706 }
707
708 }
709
710 /** Parse a scaled coefficient.
711 * @param matcher line matcher holding the coefficient
712 * @param group group number of the coefficient, or -1 if line does not contain coefficient
713 * @param scale scaling factor to apply
714 * @return scaled factor, or 0.0 if group is -1
715 */
716 private double parseCoefficient(final Matcher matcher, final int group, final double scale) {
717 if (group < 0) {
718 return 0.0;
719 } else {
720 return scale * Double.parseDouble(matcher.group(group));
721 }
722 }
723
724 /** Compute Doodson number from Delaunay multipliers.
725 * @param cGamma coefficient for γ = GMST + π tide parameter
726 * @param cL coefficient for mean anomaly of the Moon
727 * @param cLPrime coefficient for mean anomaly of the Sun
728 * @param cF coefficient for L - Ω where L is the mean longitude of the Moon
729 * @param cD coefficient for mean elongation of the Moon from the Sun
730 * @param cOmega coefficient for mean longitude of the ascending node of the Moon
731 * @return computed Doodson number
732 */
733 private int delaunayToDoodsonNumber(final int cGamma,
734 final int cL, final int cLPrime, final int cF,
735 final int cD, final int cOmega) {
736
737 // reconstruct Doodson multipliers from gamma and Delaunay multipliers
738 final int cTau = cGamma;
739 final int cS = cGamma + (cL + cF + cD);
740 final int cH = cLPrime - cD;
741 final int cP = -cL;
742 final int cNprime = cF - cOmega;
743 final int cPs = -cLPrime;
744
745 return doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs);
746
747 }
748
749 /** Compute Doodson number from Doodson multipliers.
750 * @param cTau coefficient for mean lunar time
751 * @param cS coefficient for mean longitude of the Moon
752 * @param cH coefficient for mean longitude of the Sun
753 * @param cP coefficient for longitude of Moon mean perigee
754 * @param cNprime negative of the longitude of the Moon's mean ascending node on the ecliptic
755 * @param cPs coefficient for longitude of Sun mean perigee
756 * @return computed Doodson number
757 */
758 private int doodsonToDoodsonNumber(final int cTau,
759 final int cS, final int cH, final int cP,
760 final int cNprime, final int cPs) {
761
762 return ((((cTau * 10 + (cS + 5)) * 10 + (cH + 5)) * 10 + (cP + 5)) * 10 + (cNprime + 5)) * 10 + (cPs + 5);
763
764 }
765
766 }