PoissonSeriesParser.java
/* Copyright 2002-2024 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.data;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.nio.charset.StandardCharsets;
import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import org.hipparchus.exception.DummyLocalizable;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
/**
* Parser for {@link PoissonSeries Poisson series} files.
* <p>
* A Poisson series is composed of a time polynomial part and a non-polynomial
* part which consist in summation series. The {@link SeriesTerm series terms}
* are harmonic functions (combination of sines and cosines) of polynomial
* <em>arguments</em>. The polynomial arguments are combinations of luni-solar or
* planetary {@link BodiesElements elements}.
* </p>
* <p>
* The Poisson series files from IERS have various formats, with or without
* polynomial part, with or without planetary components, with or without
* period column, with terms of increasing degrees either in dedicated columns
* or in successive sections of the file ... This class attempts to read all the
* commonly found formats, by specifying the columns of interest.
* </p>
* <p>
* The handling of increasing degrees terms (i.e. sin, cos, t sin, t cos, t^2 sin,
* t^2 cos ...) is done as follows.
* </p>
* <ul>
* <li>user must specify pairs of columns to be extracted at each line,
* in increasing degree order</li>
* <li>negative columns indices correspond to inexistent values that will be
* replaced by 0.0)</li>
* <li>file may provide section headers to specify a degree, which is added
* to the current column degree</li>
* </ul>
* <p>
* A file from an old convention, like table 5.1 in IERS conventions 1996, uses
* separate columns for degree 0 and degree 1, and uses only sine for nutation in
* longitude and cosine for nutation in obliquity. It reads as follows:
* </p>
* <pre>
* ∆ψ = Σ (Ai+A'it) sin(ARGUMENT), ∆ε = Σ (Bi+B'it) cos(ARGUMENT)
*
* MULTIPLIERS OF PERIOD LONGITUDE OBLIQUITY
* l l' F D Om days Ai A'i Bi B'i
*
* 0 0 0 0 1 -6798.4 -171996 -174.2 92025 8.9
* 0 0 2 -2 2 182.6 -13187 -1.6 5736 -3.1
* 0 0 2 0 2 13.7 -2274 -0.2 977 -0.5
* 0 0 0 0 2 -3399.2 2062 0.2 -895 0.5
* </pre>
* <p>
* In order to parse the nutation in longitude from the previous table, the
* following settings should be used:
* </p>
* <ul>
* <li>totalColumns = 10 (see {@link #PoissonSeriesParser(int)})</li>
* <li>firstDelaunay = 1 (see {@link #withFirstDelaunay(int)})</li>
* <li>no calls to {@link #withFirstPlanetary(int)} as there are no planetary columns in this table</li>
* <li>sinCosColumns = 7, -1 for degree 0 for Ai (see {@link #withSinCos(int, int, double, int, double)})</li>
* <li>sinCosColumns = 8, -1 for degree 1 for A'i (see {@link #withSinCos(int, int, double, int, double)})</li>
* </ul>
* <p>
* In order to parse the nutation in obliquity from the previous table, the
* following settings should be used:
* </p>
* <ul>
* <li>totalColumns = 10 (see {@link #PoissonSeriesParser(int)})</li>
* <li>firstDelaunay = 1 (see {@link #withFirstDelaunay(int)})</li>
* <li>no calls to {@link #withFirstPlanetary(int)} as there are no planetary columns in this table</li>
* <li>sinCosColumns = -1, 9 for degree 0 for Bi (see {@link #withSinCos(int, int, double, int, double)})</li>
* <li>sinCosColumns = -1, 10 for degree 1 for B'i (see {@link #withSinCos(int, int, double, int, double)})</li>
* </ul>
* <p>
* A file from a recent convention, like table 5.3a in IERS conventions 2010, uses
* only two columns for sin and cos, and separate degrees in successive sections with
* dedicated headers. It reads as follows:
* </p>
* <pre>
* ---------------------------------------------------------------------------------------------------
*
* (unit microarcsecond; cut-off: 0.1 microarcsecond)
* (ARG being for various combination of the fundamental arguments of the nutation theory)
*
* Sum_i[A_i * sin(ARG) + A"_i * cos(ARG)]
*
* + Sum_i[A'_i * sin(ARG) + A"'_i * cos(ARG)] * t (see Chapter 5, Eq. (35))
*
* The Table below provides the values for A_i and A"_i (j=0) and then A'_i and A"'_i (j=1)
*
* The expressions for the fundamental arguments appearing in columns 4 to 8 (luni-solar part)
* and in columns 9 to 17 (planetary part) are those of the IERS Conventions 2003
*
* ----------------------------------------------------------------------------------------------------------
* j = 0 Number of terms = 1320
* ----------------------------------------------------------------------------------------------------------
* 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
* ----------------------------------------------------------------------------------------------------------
* 1 -17206424.18 3338.60 0 0 0 0 1 0 0 0 0 0 0 0 0 0
* 2 -1317091.22 -1369.60 0 0 2 -2 2 0 0 0 0 0 0 0 0 0
* 3 -227641.81 279.60 0 0 2 0 2 0 0 0 0 0 0 0 0 0
* 4 207455.40 -69.80 0 0 0 0 2 0 0 0 0 0 0 0 0 0
* 5 147587.70 1181.70 0 1 0 0 0 0 0 0 0 0 0 0 0 0
*
* ...
*
* 1319 -0.10 0.00 0 0 0 0 0 1 0 -3 0 0 0 0 0 -2
* 1320 -0.10 0.00 0 0 0 0 0 0 0 1 0 1 -2 0 0 0
*
* --------------------------------------------------------------------------------------------------------------
* j = 1 Number of terms = 38
* --------------------------------------------------------------------------------------------------------------
* 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
* --------------------------------------------------------------------------------------------------------------
* 1321 -17418.82 2.89 0 0 0 0 1 0 0 0 0 0 0 0 0 0
* 1322 -363.71 -1.50 0 1 0 0 0 0 0 0 0 0 0 0 0 0
* 1323 -163.84 1.20 0 0 2 -2 2 0 0 0 0 0 0 0 0 0
* 1324 122.74 0.20 0 1 2 -2 2 0 0 0 0 0 0 0 0 0
* </pre>
* <p>
* In order to parse the nutation in longitude from the previous table, the
* following settings should be used:
* </p>
* <ul>
* <li>totalColumns = 17 (see {@link #PoissonSeriesParser(int)})</li>
* <li>firstDelaunay = 4 (see {@link #withFirstDelaunay(int)})</li>
* <li>firstPlanetary = 9 (see {@link #withFirstPlanetary(int)})</li>
* <li>sinCosColumns = 2,3 (we specify only degree 0, so when we read
* section j = 0 we read degree 0, when we read section j = 1 we read
* degree 1, see {@link #withSinCos(int, int, double, int, double)} ...)</li>
* </ul>
* <p>
* A file from a recent convention, like table 6.5a in IERS conventions 2010, contains
* both Doodson arguments (τ, s, h, p, N', ps), Doodson numbers and Delaunay parameters.
* In this case, the coefficients for the Delaunay parameters must be <em>subtracted</em>
* from the τ = GMST + π tide parameter, so the signs in the files must be reversed
* in order to match the Doodson arguments and Doodson numbers. This is done automatically
* (and consistency is checked) only when the {@link #withDoodson(int, int)} method is
* called at parser configuration time. Some other files use the γ = GMST + π tide parameter
* rather than Doodson τ argument and the coefficients for the Delaunay parameters must be
* <em>added</em> to the γ parameter, so no sign reversal is performed. In order to avoid
* ambiguity as the two cases are incompatible with each other, trying to add a configuration
* for τ by calling {@link #withDoodson(int, int)} and to also add a configuration for γ by
* calling {@link #withGamma(int)} triggers an exception.
* </p>
* <p>The table 6.5a file also contains a column for the waves names (the Darwin's symbol)
* which may be empty, so it must be identified explicitly by calling {@link
* #withOptionalColumn(int)}. The 6.5a table reads as follows:
* </p>
* <pre>
* The in-phase (ip) amplitudes (A₁ δkfR Hf) and the out-of-phase (op) amplitudes (A₁ δkfI Hf)
* of the corrections for frequency dependence of k₂₁⁽⁰⁾, taking the nominal value k₂₁ for the
* diurnal tides as (0.29830 − i 0.00144). Units: 10⁻¹² . The entries for δkfR and δkfI are in
* units of 10⁻⁵. Multipliers of the Doodson arguments identifying the tidal terms are given,
* as also those of the Delaunay variables characterizing the nutations produced by these
* terms.
*
* Name deg/hr Doodson τ s h p N' ps l l' F D Ω δkfR δkfI Amp. Amp.
* No. /10−5 /10−5 (ip) (op)
* 2Q₁ 12.85429 125,755 1 -3 0 2 0 0 2 0 2 0 2 -29 3 -0.1 0.0
* σ₁ 12.92714 127,555 1 -3 2 0 0 0 0 0 2 2 2 -30 3 -0.1 0.0
* 13.39645 135,645 1 -2 0 1 -1 0 1 0 2 0 1 -45 5 -0.1 0.0
* Q₁ 13.39866 135,655 1 -2 0 1 0 0 1 0 2 0 2 -46 5 -0.7 0.1
* ρ₁ 13.47151 137,455 1 -2 2 -1 0 0 -1 0 2 2 2 -49 5 -0.1 0.0
* </pre>
* <ul>
* <li>totalColumns = 18 (see {@link #PoissonSeriesParser(int)})</li>
* <li>optionalColumn = 1 (see {@link #withOptionalColumn(int)})</li>
* <li>firstDoodson, Doodson number = 4, 3 (see {@link #withDoodson(int, int)})</li>
* <li>firstDelaunay = 10 (see {@link #withFirstDelaunay(int)})</li>
* <li>sinCosColumns = 17, 18, see {@link #withSinCos(int, int, double, int, double)} ...)</li>
* </ul>
* <p>
* Our parsing algorithm involves adding the section degree from the "j = 0, 1, 2 ..." header
* to the column degree. A side effect of this algorithm is that it is theoretically possible
* to mix both formats and have for example degree two term appear as degree 2 column in section
* j=0 and as degree 1 column in section j=1 and as degree 0 column in section j=2. This case
* is not expected to be encountered in practice. The real files use either several columns
* <em>or</em> several sections, but not both at the same time.
* </p>
*
* @author Luc Maisonobe
* @see SeriesTerm
* @see PolynomialNutation
* @since 6.1
*/
public class PoissonSeriesParser {
/** Default pattern for fields with unknown type (non-space characters). */
private static final String UNKNOWN_TYPE_PATTERN = "\\S+";
/** Pattern for optional fields (either nothing or non-space characters). */
private static final String OPTIONAL_FIELD_PATTERN = "\\S*";
/** Pattern for fields with integer type. */
private static final String INTEGER_TYPE_PATTERN = "[-+]?\\p{Digit}+";
/** Pattern for fields with real type. */
private static final String REAL_TYPE_PATTERN = "[-+]?(?:(?:\\p{Digit}+(?:\\.\\p{Digit}*)?)|(?:\\.\\p{Digit}+))(?:[eE][-+]?\\p{Digit}+)?";
/** Pattern for fields with Doodson number. */
private static final String DOODSON_TYPE_PATTERN = "\\p{Digit}{2,3}[.,]\\p{Digit}{3}";
/** Pattern for String replacement. */
private static final Pattern PATTERN = Pattern.compile("[.,]");
/** Parser for the polynomial part. */
private final PolynomialParser polynomialParser;
/** Fields patterns. */
private final String[] fieldsPatterns;
/** Optional column (counting from 1). */
private final int optional;
/** Column of the γ = GMST + π tide multiplier (counting from 1). */
private final int gamma;
/** Column of the first Doodson multiplier (counting from 1). */
private final int firstDoodson;
/** Column of the Doodson number (counting from 1). */
private final int doodson;
/** Column of the first Delaunay multiplier (counting from 1). */
private final int firstDelaunay;
/** Column of the first planetary multiplier (counting from 1). */
private final int firstPlanetary;
/** columns of the sine and cosine coefficients for successive degrees.
* <p>
* The ordering is: sin, cos, t sin, t cos, t^2 sin, t^2 cos ...
* </p>
*/
private final int[] sinCosColumns;
/** Multiplicative factors to use for various columns. */
private final double[] sinCosFactors;
/** Build a parser for a Poisson series from an IERS table file.
* @param polynomialParser polynomial parser to use
* @param fieldsPatterns patterns for fields
* @param optional optional column
* @param gamma column of the GMST tide multiplier
* @param firstDoodson column of the first Doodson multiplier
* @param doodson column of the Doodson number
* @param firstDelaunay column of the first Delaunay multiplier
* @param firstPlanetary column of the first planetary multiplier
* @param sinCosColumns columns of the sine and cosine coefficients
* @param factors multiplicative factors to use for various columns
*/
private PoissonSeriesParser(final PolynomialParser polynomialParser,
final String[] fieldsPatterns, final int optional,
final int gamma, final int firstDoodson,
final int doodson, final int firstDelaunay,
final int firstPlanetary, final int[] sinCosColumns,
final double[] factors) {
this.polynomialParser = polynomialParser;
this.fieldsPatterns = fieldsPatterns;
this.optional = optional;
this.gamma = gamma;
this.firstDoodson = firstDoodson;
this.doodson = doodson;
this.firstDelaunay = firstDelaunay;
this.firstPlanetary = firstPlanetary;
this.sinCosColumns = sinCosColumns;
this.sinCosFactors = factors;
}
/** Build a parser for a Poisson series from an IERS table file.
* @param totalColumns total number of columns in the non-polynomial sections
*/
public PoissonSeriesParser(final int totalColumns) {
this(null, createInitialFieldsPattern(totalColumns), -1,
-1, -1, -1, -1, -1, new int[0], new double[0]);
}
/** Create an array with only non-space fields patterns.
* @param totalColumns total number of columns
* @return a new fields pattern array
*/
private static String[] createInitialFieldsPattern(final int totalColumns) {
final String[] patterns = new String[totalColumns];
setPatterns(patterns, 1, totalColumns, UNKNOWN_TYPE_PATTERN);
return patterns;
}
/** Set fields patterns.
* @param array fields pattern array to modify
* @param first first column to set (counting from 1), do nothing if non-positive
* @param count number of columns to set
* @param pattern pattern to use
*/
private static void setPatterns(final String[] array, final int first, final int count,
final String pattern) {
if (first > 0) {
Arrays.fill(array, first - 1, first - 1 + count, pattern);
}
}
/** Set up polynomial part parsing.
* @param freeVariable name of the free variable in the polynomial part
* @param unit default unit for polynomial, if not explicit within the file
* @return a new parser, with polynomial parser updated
*/
public PoissonSeriesParser withPolynomialPart(final char freeVariable, final PolynomialParser.Unit unit) {
return new PoissonSeriesParser(new PolynomialParser(freeVariable, unit), fieldsPatterns, optional,
gamma, firstDoodson, doodson, firstDelaunay,
firstPlanetary, sinCosColumns, sinCosFactors);
}
/** Set up optional column.
* <p>
* Optional columns typically appears in tides-related files, as some waves have
* specific names (χ₁, M₂, ...) and other waves don't have names and hence are
* replaced by spaces in the corresponding file line.
* </p>
* <p>
* At most one column may be optional.
* </p>
* @param column optional column (counting from 1)
* @return a new parser, with updated columns settings
*/
public PoissonSeriesParser withOptionalColumn(final int column) {
// update the fields pattern to expect 1 optional field at the right index
final String[] newFieldsPatterns = fieldsPatterns.clone();
setPatterns(newFieldsPatterns, optional, 1, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, column, 1, OPTIONAL_FIELD_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, column,
gamma, firstDoodson, doodson, firstDelaunay,
firstPlanetary, sinCosColumns, sinCosFactors);
}
/** Set up column of GMST tide multiplier.
* @param column column of the GMST tide multiplier (counting from 1)
* @return a new parser, with updated columns settings
* @see #withDoodson(int, int)
*/
public PoissonSeriesParser withGamma(final int column) {
// check we don't try to have both τ and γ configured at the same time
if (firstDoodson > 0 && column > 0) {
throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
}
// update the fields pattern to expect 1 integer at the right index
final String[] newFieldsPatterns = fieldsPatterns.clone();
setPatterns(newFieldsPatterns, gamma, 1, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, column, 1, INTEGER_TYPE_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
column, firstDoodson, doodson, firstDelaunay,
firstPlanetary, sinCosColumns, sinCosFactors);
}
/** Set up columns for Doodson multipliers and Doodson number.
* @param firstMultiplierColumn column of the first Doodson multiplier which
* corresponds to τ (counting from 1)
* @param numberColumn column of the Doodson number (counting from 1)
* @return a new parser, with updated columns settings
* @see #withGamma(int)
* @see #withFirstDelaunay(int)
*/
public PoissonSeriesParser withDoodson(final int firstMultiplierColumn, final int numberColumn) {
// check we don't try to have both τ and γ configured at the same time
if (gamma > 0 && firstMultiplierColumn > 0) {
throw new OrekitException(OrekitMessages.CANNOT_PARSE_BOTH_TAU_AND_GAMMA);
}
final String[] newFieldsPatterns = fieldsPatterns.clone();
// update the fields pattern to expect 6 integers at the right indices
setPatterns(newFieldsPatterns, firstDoodson, 6, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, firstMultiplierColumn, 6, INTEGER_TYPE_PATTERN);
// update the fields pattern to expect 1 Doodson number at the right index
setPatterns(newFieldsPatterns, doodson, 1, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, numberColumn, 1, DOODSON_TYPE_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
gamma, firstMultiplierColumn, numberColumn, firstDelaunay,
firstPlanetary, sinCosColumns, sinCosFactors);
}
/** Set up first column of Delaunay multiplier.
* @param firstColumn column of the first Delaunay multiplier (counting from 1)
* @return a new parser, with updated columns settings
*/
public PoissonSeriesParser withFirstDelaunay(final int firstColumn) {
// update the fields pattern to expect 5 integers at the right indices
final String[] newFieldsPatterns = fieldsPatterns.clone();
setPatterns(newFieldsPatterns, firstDelaunay, 5, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, firstColumn, 5, INTEGER_TYPE_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
gamma, firstDoodson, doodson, firstColumn,
firstPlanetary, sinCosColumns, sinCosFactors);
}
/** Set up first column of planetary multiplier.
* @param firstColumn column of the first planetary multiplier (counting from 1)
* @return a new parser, with updated columns settings
*/
public PoissonSeriesParser withFirstPlanetary(final int firstColumn) {
// update the fields pattern to expect 9 integers at the right indices
final String[] newFieldsPatterns = fieldsPatterns.clone();
setPatterns(newFieldsPatterns, firstPlanetary, 9, UNKNOWN_TYPE_PATTERN);
setPatterns(newFieldsPatterns, firstColumn, 9, INTEGER_TYPE_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
gamma, firstDoodson, doodson, firstDelaunay,
firstColumn, sinCosColumns, sinCosFactors);
}
/** Set up columns of the sine and cosine coefficients.
* @param degree degree to set up
* @param sinColumn column of the sine coefficient for t<sup>degree</sup> counting from 1
* (may be -1 if there are no sine coefficients)
* @param sinFactor multiplicative factor for the sine coefficient
* @param cosColumn column of the cosine coefficient for t<sup>degree</sup> counting from 1
* (may be -1 if there are no cosine coefficients)
* @param cosFactor multiplicative factor for the cosine coefficient
* @return a new parser, with updated columns settings
*/
public PoissonSeriesParser withSinCos(final int degree,
final int sinColumn, final double sinFactor,
final int cosColumn, final double cosFactor) {
// update the sin/cos columns array
final int maxDegree = FastMath.max(degree, sinCosColumns.length / 2 - 1);
final int[] newSinCosColumns = new int[2 * (maxDegree + 1)];
Arrays.fill(newSinCosColumns, -1);
System.arraycopy(sinCosColumns, 0, newSinCosColumns, 0, sinCosColumns.length);
newSinCosColumns[2 * degree] = sinColumn;
newSinCosColumns[2 * degree + 1] = cosColumn;
final double[] newSinCosFactors = new double[2 * (maxDegree + 1)];
Arrays.fill(newSinCosFactors, Double.NaN);
System.arraycopy(sinCosFactors, 0, newSinCosFactors, 0, sinCosFactors.length);
newSinCosFactors[2 * degree] = sinFactor;
newSinCosFactors[2 * degree + 1] = cosFactor;
// update the fields pattern to expect real numbers at the right indices
final String[] newFieldsPatterns = fieldsPatterns.clone();
if (2 * degree < sinCosColumns.length) {
setPatterns(newFieldsPatterns, sinCosColumns[2 * degree], 1, UNKNOWN_TYPE_PATTERN);
}
setPatterns(newFieldsPatterns, sinColumn, 1, REAL_TYPE_PATTERN);
if (2 * degree + 1 < sinCosColumns.length) {
setPatterns(newFieldsPatterns, sinCosColumns[2 * degree + 1], 1, UNKNOWN_TYPE_PATTERN);
}
setPatterns(newFieldsPatterns, cosColumn, 1, REAL_TYPE_PATTERN);
return new PoissonSeriesParser(polynomialParser, newFieldsPatterns, optional,
gamma, firstDoodson, doodson, firstDelaunay,
firstPlanetary, newSinCosColumns, newSinCosFactors);
}
/** Parse a stream.
* @param stream stream containing the IERS table
* @param name name of the resource file (for error messages only)
* @return parsed Poisson series
*/
public PoissonSeries parse(final InputStream stream, final String name) {
if (stream == null) {
throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
}
// the degrees section header should read something like:
// j = 0 Nb of terms = 1306
// or something like:
// j = 0 Number of terms = 1037
final Pattern degreeSectionHeaderPattern =
Pattern.compile("^\\p{Space}*j\\p{Space}*=\\p{Space}*(\\p{Digit}+)" +
"[\\p{Alpha}\\p{Space}]+=\\p{Space}*(\\p{Digit}+)\\p{Space}*$");
// regular lines are simply a space separated list of numbers
final StringBuilder builder = new StringBuilder("^\\p{Space}*");
for (int i = 0; i < fieldsPatterns.length; ++i) {
builder.append("(");
builder.append(fieldsPatterns[i]);
builder.append(")");
builder.append((i < fieldsPatterns.length - 1) ? "\\p{Space}+" : "\\p{Space}*$");
}
final Pattern regularLinePattern = Pattern.compile(builder.toString());
// setup the reader
try (BufferedReader reader = new BufferedReader(new InputStreamReader(stream, StandardCharsets.UTF_8))) {
int lineNumber = 0;
int expectedIndex = -1;
int nTerms = -1;
int count = 0;
int degree = 0;
// prepare the container for the parsed data
PolynomialNutation polynomial;
if (polynomialParser == null) {
// we don't expect any polynomial, we directly set the zero polynomial
polynomial = new PolynomialNutation(new double[0]);
} else {
// the dedicated parser will fill in the polynomial later
polynomial = null;
}
final Map<Long, SeriesTerm> series = new HashMap<Long, SeriesTerm>();
for (String line = reader.readLine(); line != null; line = reader.readLine()) {
// replace unicode minus sign ('−') by regular hyphen ('-') for parsing
// such unicode characters occur in tables that are copy-pasted from PDF files
line = line.replace('\u2212', '-');
++lineNumber;
final Matcher regularMatcher = regularLinePattern.matcher(line);
if (regularMatcher.matches()) {
// we have found a regular data line
if (expectedIndex > 0) {
// we are in a file were terms are numbered, we check the index
if (Integer.parseInt(regularMatcher.group(1)) != expectedIndex) {
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, regularMatcher.group());
}
}
// get the Doodson multipliers as well as the Doodson number
final int cTau = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson));
final int cS = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 1));
final int cH = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 2));
final int cP = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 3));
final int cNprime = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 4));
final int cPs = (firstDoodson < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstDoodson + 5));
final int nDoodson = (doodson < 0) ? 0 : Integer.parseInt(PATTERN.matcher(regularMatcher.group(doodson)).replaceAll(""));
// get the tide multiplier
int cGamma = (gamma < 0) ? 0 : Integer.parseInt(regularMatcher.group(gamma));
// get the Delaunay multipliers
int cL = Integer.parseInt(regularMatcher.group(firstDelaunay));
int cLPrime = Integer.parseInt(regularMatcher.group(firstDelaunay + 1));
int cF = Integer.parseInt(regularMatcher.group(firstDelaunay + 2));
int cD = Integer.parseInt(regularMatcher.group(firstDelaunay + 3));
int cOmega = Integer.parseInt(regularMatcher.group(firstDelaunay + 4));
// get the planetary multipliers
final int cMe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary));
final int cVe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 1));
final int cE = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 2));
final int cMa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 3));
final int cJu = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 4));
final int cSa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 5));
final int cUr = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 6));
final int cNe = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 7));
final int cPa = (firstPlanetary < 0) ? 0 : Integer.parseInt(regularMatcher.group(firstPlanetary + 8));
if (nDoodson > 0) {
// set up the traditional parameters corresponding to the Doodson arguments
cGamma = cTau;
cL = -cL;
cLPrime = -cLPrime;
cF = -cF;
cD = -cD;
cOmega = -cOmega;
// check Doodson number, Doodson multipliers and Delaunay multipliers consistency
if (nDoodson != doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs) ||
nDoodson != delaunayToDoodsonNumber(cGamma, cL, cLPrime, cF, cD, cOmega)) {
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, regularMatcher.group());
}
}
final long key = NutationCodec.encode(cGamma, cL, cLPrime, cF, cD, cOmega,
cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
// retrieved the term, or build it if it's the first time it is encountered in the file
final SeriesTerm term;
if (series.containsKey(key)) {
// the term was already known, from another degree
term = series.get(key);
} else {
// the term is a new one
term = SeriesTerm.buildTerm(cGamma, cL, cLPrime, cF, cD, cOmega,
cMe, cVe, cE, cMa, cJu, cSa, cUr, cNe, cPa);
}
boolean nonZero = false;
for (int d = 0; d < sinCosColumns.length / 2; ++d) {
final double sinCoeff =
parseCoefficient(regularMatcher, sinCosColumns[2 * d], sinCosFactors[2 * d]);
final double cosCoeff =
parseCoefficient(regularMatcher, sinCosColumns[2 * d + 1], sinCosFactors[2 * d + 1]);
if (!Precision.equals(sinCoeff, 0.0, 0) || !Precision.equals(cosCoeff, 0.0, 0)) {
nonZero = true;
term.add(0, degree + d, sinCoeff, cosCoeff);
++count;
}
}
if (nonZero) {
series.put(key, term);
}
if (expectedIndex > 0) {
// we are in a file were terms are numbered
// we must update the expected value for next term
++expectedIndex;
}
} else {
final Matcher headerMatcher = degreeSectionHeaderPattern.matcher(line);
if (headerMatcher.matches()) {
// we have found a degree section header
final int nextDegree = Integer.parseInt(headerMatcher.group(1));
if (nextDegree != degree + 1 && (degree != 0 || nextDegree != 0)) {
throw new OrekitException(OrekitMessages.MISSING_SERIE_J_IN_FILE,
degree + 1, name, lineNumber);
}
if (nextDegree == 0) {
// in IERS files split in sections, all terms are numbered
// we can check the indices
expectedIndex = 1;
}
if (nextDegree > 0 && count != nTerms) {
// the previous degree does not have the expected number of terms
throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
}
// remember the number of terms the upcoming sublist should have
nTerms = Integer.parseInt(headerMatcher.group(2));
count = 0;
degree = nextDegree;
} else if (polynomial == null) {
// look for the polynomial part
final double[] coefficients = polynomialParser.parse(line);
if (coefficients != null) {
polynomial = new PolynomialNutation(coefficients);
}
}
}
}
if (polynomial == null || series.isEmpty()) {
throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
}
if (nTerms > 0 && count != nTerms) {
// the last degree does not have the expected number of terms
throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, name);
}
// build the series
return new PoissonSeries(polynomial, series);
} catch (IOException ioe) {
throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
}
}
/** Parse a scaled coefficient.
* @param matcher line matcher holding the coefficient
* @param group group number of the coefficient, or -1 if line does not contain coefficient
* @param scale scaling factor to apply
* @return scaled factor, or 0.0 if group is -1
*/
private double parseCoefficient(final Matcher matcher, final int group, final double scale) {
if (group < 0) {
return 0.0;
} else {
return scale * Double.parseDouble(matcher.group(group));
}
}
/** Compute Doodson number from Delaunay multipliers.
* @param cGamma coefficient for γ = GMST + π tide parameter
* @param cL coefficient for mean anomaly of the Moon
* @param cLPrime coefficient for mean anomaly of the Sun
* @param cF coefficient for L - Ω where L is the mean longitude of the Moon
* @param cD coefficient for mean elongation of the Moon from the Sun
* @param cOmega coefficient for mean longitude of the ascending node of the Moon
* @return computed Doodson number
*/
private int delaunayToDoodsonNumber(final int cGamma,
final int cL, final int cLPrime, final int cF,
final int cD, final int cOmega) {
// reconstruct Doodson multipliers from gamma and Delaunay multipliers
final int cTau = cGamma;
final int cS = cGamma + (cL + cF + cD);
final int cH = cLPrime - cD;
final int cP = -cL;
final int cNprime = cF - cOmega;
final int cPs = -cLPrime;
return doodsonToDoodsonNumber(cTau, cS, cH, cP, cNprime, cPs);
}
/** Compute Doodson number from Doodson multipliers.
* @param cTau coefficient for mean lunar time
* @param cS coefficient for mean longitude of the Moon
* @param cH coefficient for mean longitude of the Sun
* @param cP coefficient for longitude of Moon mean perigee
* @param cNprime negative of the longitude of the Moon's mean ascending node on the ecliptic
* @param cPs coefficient for longitude of Sun mean perigee
* @return computed Doodson number
*/
private int doodsonToDoodsonNumber(final int cTau,
final int cS, final int cH, final int cP,
final int cNprime, final int cPs) {
return ((((cTau * 10 + (cS + 5)) * 10 + (cH + 5)) * 10 + (cP + 5)) * 10 + (cNprime + 5)) * 10 + (cPs + 5);
}
}