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