1 /* Copyright 2002-2017 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 colums 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 column of the GMST tide multiplier (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 * @exception OrekitException if τ has been configured by a previous call
366 * to {@link #withDoodson(int, int)}
367 * @see #withDoodson(int, int)
368 */
369 public PoissonSeriesParser withGamma(final int column) throws OrekitException {
370
371 // check we don't try to have both τ and γ configured at the same time
372 if (firstDoodson > 0 && column > 0) {
373 throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
374 }
375
376 // update the fields pattern to expect 1 integer at the right index
377 final String[] newFieldsPatterns = fieldsPatterns.clone();
378 setPatterns(newFieldsPatterns, gamma, 1, UNKNOWN_TYPE_PATTERN);
379 setPatterns(newFieldsPatterns, column, 1, INTEGER_TYPE_PATTERN);
380
381 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
382 column, firstDoodson, doodson, firstDelaunay,
383 firstPlanetary, sinCosColumns, sinCosFactors);
384
385 }
386
387 /** Set up columns for Doodson multiplers and Doodson number.
388 * @param firstMultiplierColumn column of the first Doodson multiplier which
389 * corresponds to τ (counting from 1)
390 * @param numberColumn column of the Doodson number (counting from 1)
391 * @return a new parser, with updated columns settings
392 * @exception OrekitException if γ has been configured by a previous call
393 * to {@link #withGamma(int)}
394 * @see #withGamma(int)
395 * @see #withFirstDelaunay(int)
396 */
397 public PoissonSeriesParser withDoodson(final int firstMultiplierColumn, final int numberColumn)
398 throws OrekitException {
399
400 // check we don't try to have both τ and γ configured at the same time
401 if (gamma > 0 && firstMultiplierColumn > 0) {
402 throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
403 }
404
405 final String[] newFieldsPatterns = fieldsPatterns.clone();
406
407 // update the fields pattern to expect 6 integers at the right indices
408 setPatterns(newFieldsPatterns, firstDoodson, 6, UNKNOWN_TYPE_PATTERN);
409 setPatterns(newFieldsPatterns, firstMultiplierColumn, 6, INTEGER_TYPE_PATTERN);
410
411 // update the fields pattern to expect 1 Doodson number at the right index
412 setPatterns(newFieldsPatterns, doodson, 1, UNKNOWN_TYPE_PATTERN);
413 setPatterns(newFieldsPatterns, numberColumn, 1, DOODSON_TYPE_PATTERN);
414
415 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
416 gamma, firstMultiplierColumn, numberColumn, firstDelaunay,
417 firstPlanetary, sinCosColumns, sinCosFactors);
418
419 }
420
421 /** Set up first column of Delaunay multiplier.
422 * @param firstColumn column of the first Delaunay multiplier (counting from 1)
423 * @return a new parser, with updated columns settings
424 */
425 public PoissonSeriesParser withFirstDelaunay(final int firstColumn) {
426
427 // update the fields pattern to expect 5 integers at the right indices
428 final String[] newFieldsPatterns = fieldsPatterns.clone();
429 setPatterns(newFieldsPatterns, firstDelaunay, 5, UNKNOWN_TYPE_PATTERN);
430 setPatterns(newFieldsPatterns, firstColumn, 5, INTEGER_TYPE_PATTERN);
431
432 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
433 gamma, firstDoodson, doodson, firstColumn,
434 firstPlanetary, sinCosColumns, sinCosFactors);
435
436 }
437
438 /** Set up first column of planetary multiplier.
439 * @param firstColumn column of the first planetary multiplier (counting from 1)
440 * @return a new parser, with updated columns settings
441 */
442 public PoissonSeriesParser withFirstPlanetary(final int firstColumn) {
443
444 // update the fields pattern to expect 9 integers at the right indices
445 final String[] newFieldsPatterns = fieldsPatterns.clone();
446 setPatterns(newFieldsPatterns, firstPlanetary, 9, UNKNOWN_TYPE_PATTERN);
447 setPatterns(newFieldsPatterns, firstColumn, 9, INTEGER_TYPE_PATTERN);
448
449 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
450 gamma, firstDoodson, doodson, firstDelaunay,
451 firstColumn, sinCosColumns, sinCosFactors);
452
453 }
454
455 /** Set up columns of the sine and cosine coefficients.
456 * @param degree degree to set up
457 * @param sinColumn column of the sine coefficient for t<sup>degree</sup> counting from 1
458 * (may be -1 if there are no sine coefficients)
459 * @param sinFactor multiplicative factor for the sine coefficient
460 * @param cosColumn column of the cosine coefficient for t<sup>degree</sup> counting from 1
461 * (may be -1 if there are no cosine coefficients)
462 * @param cosFactor multiplicative factor for the cosine coefficient
463 * @return a new parser, with updated columns settings
464 */
465 public PoissonSeriesParser withSinCos(final int degree,
466 final int sinColumn, final double sinFactor,
467 final int cosColumn, final double cosFactor) {
468
469 // update the sin/cos columns array
470 final int maxDegree = FastMath.max(degree, sinCosColumns.length / 2 - 1);
471 final int[] newSinCosColumns = new int[2 * (maxDegree + 1)];
472 Arrays.fill(newSinCosColumns, -1);
473 System.arraycopy(sinCosColumns, 0, newSinCosColumns, 0, sinCosColumns.length);
474 newSinCosColumns[2 * degree] = sinColumn;
475 newSinCosColumns[2 * degree + 1] = cosColumn;
476
477 final double[] newSinCosFactors = new double[2 * (maxDegree + 1)];
478 Arrays.fill(newSinCosFactors, Double.NaN);
479 System.arraycopy(sinCosFactors, 0, newSinCosFactors, 0, sinCosFactors.length);
480 newSinCosFactors[2 * degree] = sinFactor;
481 newSinCosFactors[2 * degree + 1] = cosFactor;
482
483 // update the fields pattern to expect real numbers at the right indices
484 final String[] newFieldsPatterns = fieldsPatterns.clone();
485 if (2 * degree < sinCosColumns.length) {
486 setPatterns(newFieldsPatterns, sinCosColumns[2 * degree], 1, UNKNOWN_TYPE_PATTERN);
487 }
488 setPatterns(newFieldsPatterns, sinColumn, 1, REAL_TYPE_PATTERN);
489 if (2 * degree + 1 < sinCosColumns.length) {
490 setPatterns(newFieldsPatterns, sinCosColumns[2 * degree + 1], 1, UNKNOWN_TYPE_PATTERN);
491 }
492 setPatterns(newFieldsPatterns, cosColumn, 1, REAL_TYPE_PATTERN);
493
494 return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
495 gamma, firstDoodson, doodson, firstDelaunay,
496 firstPlanetary, newSinCosColumns, newSinCosFactors);
497
498 }
499
500 /** Parse a stream.
501 * @param stream stream containing the IERS table
502 * @param name name of the resource file (for error messages only)
503 * @return parsed Poisson series
504 * @exception OrekitException if stream is null or the table cannot be parsed
505 */
506 public PoissonSeries parse(final InputStream stream, final String name) throws OrekitException {
507
508 if (stream == null) {
509 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
510 }
511
512 // the degrees section header should read something like:
513 // j = 0 Nb of terms = 1306
514 // or something like:
515 // j = 0 Number of terms = 1037
516 final Pattern degreeSectionHeaderPattern =
517 Pattern.compile("^\\p{Space}*j\\p{Space}*=\\p{Space}*(\\p{Digit}+)" +
518 "[\\p{Alpha}\\p{Space}]+=\\p{Space}*(\\p{Digit}+)\\p{Space}*$");
519
520 // regular lines are simply a space separated list of numbers
521 final StringBuilder builder = new StringBuilder("^\\p{Space}*");
522 for (int i = 0; i < fieldsPatterns.length; ++i) {
523 builder.append("(");
524 builder.append(fieldsPatterns[i]);
525 builder.append(")");
526 builder.append((i < fieldsPatterns.length - 1) ? "\\p{Space}+" : "\\p{Space}*$");
527 }
528 final Pattern regularLinePattern = Pattern.compile(builder.toString());
529
530 try {
531
532 // setup the reader
533 final BufferedReader reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"));
534 int lineNumber = 0;
535 int expectedIndex = -1;
536 int nTerms = -1;
537 int count = 0;
538 int degree = 0;
539
540 // prepare the container for the parsed data
541 PolynomialNutation polynomial;
542 if (polynomialParser == null) {
543 // we don't expect any polynomial, we directly set the zero polynomial
544 polynomial = new PolynomialNutation(new double[0]);
545 } else {
546 // the dedicated parser will fill in the polynomial later
547 polynomial = null;
548 }
549 final Map<Long, SeriesTerm> series = new HashMap<Long, SeriesTerm>();
550
551 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
552
553 // replace unicode minus sign ('−') by regular hyphen ('-') for parsing
554 // such unicode characters occur in tables that are copy-pasted from PDF files
555 line = line.replace('\u2212', '-');
556 ++lineNumber;
557
558 final Matcher regularMatcher = regularLinePattern.matcher(line);
559 if (regularMatcher.matches()) {
560 // we have found a regular data line
561
562 if (expectedIndex > 0) {
563 // we are in a file were terms are numbered, we check the index
564 if (Integer.parseInt(regularMatcher.group(1)) != expectedIndex) {
565 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
566 lineNumber, name, regularMatcher.group());
567 }
568 }
569
570 // get the Doodson multipliers as well as the Doodson number
571 final int cTau = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson));
572 final int cS = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 1));
573 final int cH = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 2));
574 final int cP = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 3));
575 final int cNprime = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 4));
576 final int cPs = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 5));
577 final int nDoodson = (doodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(doodson).replaceAll("[.,]", ""));
578
579 // get the tide multipler
580 int cGamma = (gamma < 0) ? 0 : Integer.parseInt(regularMatcher.group(gamma));
581
582 // get the Delaunay multipliers
583 int cL = Integer.parseInt(regularMatcher.group(firstDelaunay));
584 int cLPrime = Integer.parseInt(regularMatcher.group(firstDelaunay + 1));
585 int cF = Integer.parseInt(regularMatcher.group(firstDelaunay + 2));
586 int cD = Integer.parseInt(regularMatcher.group(firstDelaunay + 3));
587 int cOmega = Integer.parseInt(regularMatcher.group(firstDelaunay + 4));
588
589 // get the planetary multipliers
590 final int cMe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary));
591 final int cVe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 1));
592 final int cE = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 2));
593 final int cMa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 3));
594 final int cJu = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 4));
595 final int cSa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 5));
596 final int cUr = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 6));
597 final int cNe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 7));
598 final int cPa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 8));
599
600 if (nDoodson > 0) {
601
602 // set up the traditional parameters corresponding to the Doodson arguments
603 cGamma = cTau;
604 cL = -cL;
605 cLPrime = -cLPrime;
606 cF = -cF;
607 cD = -cD;
608 cOmega = -cOmega;
609
610 // check Doodson number, Doodson multiplers and Delaunay multipliers consistency
611 if (nDoodson != doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs) ||
612 nDoodson != delaunayToDoodsonNumber(cGamma, cL, cLPrime, cF, cD, cOmega)) {
613 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
614 lineNumber, name, regularMatcher.group());
615 }
616
617 }
618
619 final long key = NutationCodec.encode(cGamma, cL, cLPrime, cF, cD, cOmega,
620 cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
621
622 // retrieved the term, or build it if it's the first time it is encountered in the file
623 final SeriesTerm term;
624 if (series.containsKey(key)) {
625 // the term was already known, from another degree
626 term = series.get(key);
627 } else {
628 // the term is a new one
629 term = SeriesTerm.buildTerm(cGamma, cL, cLPrime, cF, cD, cOmega,
630 cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
631 }
632
633 boolean nonZero = false;
634 for (int d = 0; d < sinCosColumns.length / 2; ++d) {
635 final double sinCoeff =
636 parseCoefficient(regularMatcher, sinCosColumns[2 * d], sinCosFactors[2 * d]);
637 final double cosCoeff =
638 parseCoefficient(regularMatcher, sinCosColumns[2 * d + 1], sinCosFactors[2 * d + 1]);
639 if (!Precision.equals(sinCoeff, 0.0, 0) || !Precision.equals(cosCoeff, 0.0, 0)) {
640 nonZero = true;
641 term.add(0, degree + d, sinCoeff, cosCoeff);
642 ++count;
643 }
644 }
645 if (nonZero) {
646 series.put(key, term);
647 }
648
649 if (expectedIndex > 0) {
650 // we are in a file were terms are numbered
651 // we must update the expected value for next term
652 ++expectedIndex;
653 }
654
655 } else {
656
657 final Matcher headerMatcher = degreeSectionHeaderPattern.matcher(line);
658 if (headerMatcher.matches()) {
659
660 // we have found a degree section header
661 final int nextDegree = Integer.parseInt(headerMatcher.group(1));
662 if ((nextDegree != degree + 1) && (degree != 0 || nextDegree != 0)) {
663 throw new OrekitException(OrekitMessages.MISSING_SERIE_J_IN_FILE,
664 degree + 1, name, lineNumber);
665 }
666
667 if (nextDegree == 0) {
668 // in IERS files split in sections, all terms are numbered
669 // we can check the indices
670 expectedIndex = 1;
671 }
672
673 if (nextDegree > 0 && count != nTerms) {
674 // the previous degree does not have the expected number of terms
675 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
676 }
677
678 // remember the number of terms the upcoming sublist should have
679 nTerms = Integer.parseInt(headerMatcher.group(2));
680 count = 0;
681 degree = nextDegree;
682
683 } else if (polynomial == null) {
684 // look for the polynomial part
685 final double[] coefficients = polynomialParser.parse(line);
686 if (coefficients != null) {
687 polynomial = new PolynomialNutation(coefficients);
688 }
689 }
690
691 }
692
693 }
694
695 if (polynomial == null || series.isEmpty()) {
696 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
697 }
698
699 if (nTerms > 0 && count != nTerms) {
700 // the last degree does not have the expected number of terms
701 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
702 }
703
704 // build the series
705 return new PoissonSeries(polynomial, series);
706
707 } catch (IOException ioe) {
708 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
709 }
710
711 }
712
713 /** Parse a scaled coefficient.
714 * @param matcher line matcher holding the coefficient
715 * @param group group number of the coefficient, or -1 if line does not contain coefficient
716 * @param scale scaling factor to apply
717 * @return scaled factor, or 0.0 if group is -1
718 */
719 private double parseCoefficient(final Matcher matcher, final int group, final double scale) {
720 if (group < 0) {
721 return 0.0;
722 } else {
723 return scale * Double.parseDouble(matcher.group(group));
724 }
725 }
726
727 /** Compute Doodson number from Delaunay multipliers.
728 * @param cGamma coefficient for γ = GMST + π tide parameter
729 * @param cL coefficient for mean anomaly of the Moon
730 * @param cLPrime coefficient for mean anomaly of the Sun
731 * @param cF coefficient for L - Ω where L is the mean longitude of the Moon
732 * @param cD coefficient for mean elongation of the Moon from the Sun
733 * @param cOmega coefficient for mean longitude of the ascending node of the Moon
734 * @return computed Doodson number
735 */
736 private int delaunayToDoodsonNumber(final int cGamma,
737 final int cL, final int cLPrime, final int cF,
738 final int cD, final int cOmega) {
739
740 // reconstruct Doodson multipliers from gamma and Delaunay multipliers
741 final int cTau = cGamma;
742 final int cS = cGamma + (cL + cF + cD);
743 final int cH = cLPrime - cD;
744 final int cP = -cL;
745 final int cNprime = cF - cOmega;
746 final int cPs = -cLPrime;
747
748 return doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs);
749
750 }
751
752 /** Compute Doodson number from Doodson multipliers.
753 * @param cTau coefficient for mean lunar time
754 * @param cS coefficient for mean longitude of the Moon
755 * @param cH coefficient for mean longitude of the Sun
756 * @param cP coefficient for longitude of Moon mean perigee
757 * @param cNprime negative of the longitude of the Moon's mean ascending node on the ecliptic
758 * @param cPs coefficient for longitude of Sun mean perigee
759 * @return computed Doodson number
760 */
761 private int doodsonToDoodsonNumber(final int cTau,
762 final int cS, final int cH, final int cP,
763 final int cNprime, final int cPs) {
764
765 return ((((cTau * 10 + (cS + 5)) * 10 + (cH + 5)) * 10 + (cP + 5)) * 10 + (cNprime + 5)) * 10 + (cPs + 5);
766
767 }
768
769 }