1   /* Copyright 2002-2024 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 }