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