1 /* Copyright 2002-2016 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.io.Serializable;
24 import java.util.ArrayList;
25 import java.util.Arrays;
26 import java.util.HashMap;
27 import java.util.List;
28 import java.util.Map;
29 import java.util.regex.Matcher;
30 import java.util.regex.Pattern;
31
32 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
33 import org.apache.commons.math3.exception.util.DummyLocalizable;
34 import org.apache.commons.math3.util.FastMath;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.errors.OrekitInternalError;
37 import org.orekit.errors.OrekitMessages;
38 import org.orekit.time.AbsoluteDate;
39 import org.orekit.time.TimeFunction;
40 import org.orekit.time.TimeScale;
41 import org.orekit.utils.IERSConventions;
42
43 /**
44 * Class computing the fundamental arguments for nutation and tides.
45 * <p>
46 * The fundamental arguments are split in two sets:
47 * </p>
48 * <ul>
49 * <li>the Delaunay arguments for Moon and Sun effects</li>
50 * <li>the planetary arguments for other planets</li>
51 * </ul>
52 *
53 * @author Luc Maisonobe
54 * @see SeriesTerm
55 * @see PoissonSeries
56 * @see BodiesElements
57 */
58 public class FundamentalNutationArguments implements Serializable {
59
60 /** Serializable UID. */
61 private static final long serialVersionUID = 20131209L;
62
63 /** IERS conventions to use. */
64 private final IERSConventions conventions;
65
66 /** Time scale for GMST computation. */
67 private final TimeScale timeScale;
68
69 /** Function computing Greenwich Mean Sidereal Time. */
70 private final transient TimeFunction<DerivativeStructure> gmstFunction;
71
72 // luni-solar Delaunay arguments
73
74 /** Coefficients for mean anomaly of the Moon. */
75 private final double[] lCoefficients;
76
77 /** Coefficients for mean anomaly of the Sun. */
78 private final double[] lPrimeCoefficients;
79
80 /** Coefficients for L - Ω where L is the mean longitude of the Moon. */
81 private final double[] fCoefficients;
82
83 /** Coefficients for mean elongation of the Moon from the Sun. */
84 private final double[] dCoefficients;
85
86 /** Coefficients for mean longitude of the ascending node of the Moon. */
87 private final double[] omegaCoefficients;
88
89 // planetary nutation arguments
90
91 /** Coefficients for mean Mercury longitude. */
92 private final double[] lMeCoefficients;
93
94 /** Coefficients for mean Venus longitude. */
95 private final double[] lVeCoefficients;
96
97 /** Coefficients for mean Earth longitude. */
98 private final double[] lECoefficients;
99
100 /** Coefficients for mean Mars longitude. */
101 private final double[] lMaCoefficients;
102
103 /** Coefficients for mean Jupiter longitude. */
104 private final double[] lJCoefficients;
105
106 /** Coefficients for mean Saturn longitude. */
107 private final double[] lSaCoefficients;
108
109 /** Coefficients for mean Uranus longitude. */
110 private final double[] lUCoefficients;
111
112 /** Coefficients for mean Neptune longitude. */
113 private final double[] lNeCoefficients;
114
115 /** Coefficients for general accumulated precession. */
116 private final double[] paCoefficients;
117
118 /** Build a model of fundamental arguments from an IERS table file.
119 * @param conventions IERS conventions to use
120 * @param timeScale time scale for GMST computation
121 * (may be null if tide parameter γ = GMST + π is not needed)
122 * @param stream stream containing the IERS table
123 * @param name name of the resource file (for error messages only)
124 * @exception OrekitException if stream is null or the table cannot be parsed
125 */
126 public FundamentalNutationArguments(final IERSConventions conventions,
127 final TimeScale timeScale,
128 final InputStream stream, final String name)
129 throws OrekitException {
130 this(conventions, timeScale, parseCoefficients(stream, name));
131 }
132
133 /** Build a model of fundamental arguments from an IERS table file.
134 * @param conventions IERS conventions to use
135 * @param timeScale time scale for GMST computation
136 * (may be null if tide parameter γ = GMST + π is not needed)
137 * @param coefficients list of coefficients arrays (all 14 arrays must be provided,
138 * the 5 Delaunay first and the 9 planetary afterwards)
139 * @exception OrekitException if GMST function cannot be retrieved
140 * @since 6.1
141 */
142 public FundamentalNutationArguments(final IERSConventions conventions, final TimeScale timeScale,
143 final List<double[]> coefficients)
144 throws OrekitException {
145 this.conventions = conventions;
146 this.timeScale = timeScale;
147 this.gmstFunction = (timeScale == null) ? null : conventions.getGMSTFunction(timeScale);
148 this.lCoefficients = coefficients.get( 0);
149 this.lPrimeCoefficients = coefficients.get( 1);
150 this.fCoefficients = coefficients.get( 2);
151 this.dCoefficients = coefficients.get( 3);
152 this.omegaCoefficients = coefficients.get( 4);
153 this.lMeCoefficients = coefficients.get( 5);
154 this.lVeCoefficients = coefficients.get( 6);
155 this.lECoefficients = coefficients.get( 7);
156 this.lMaCoefficients = coefficients.get( 8);
157 this.lJCoefficients = coefficients.get( 9);
158 this.lSaCoefficients = coefficients.get(10);
159 this.lUCoefficients = coefficients.get(11);
160 this.lNeCoefficients = coefficients.get(12);
161 this.paCoefficients = coefficients.get(13);
162 }
163
164 /** Parse coefficients.
165 * @param stream stream containing the IERS table
166 * @param name name of the resource file (for error messages only)
167 * @return list of coefficients arrays
168 * @exception OrekitException if stream is null or the table cannot be parsed
169 */
170 private static List<double[]> parseCoefficients(final InputStream stream, final String name)
171 throws OrekitException {
172
173 if (stream == null) {
174 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, name);
175 }
176
177 try {
178
179 final DefinitionParser definitionParser = new DefinitionParser();
180
181 // setup the reader
182 final BufferedReader reader = new BufferedReader(new InputStreamReader(stream, "UTF-8"));
183 int lineNumber = 0;
184
185 // look for the reference date and the 14 polynomials
186 final int n = FundamentalName.values().length;
187 final Map<FundamentalName, double[]> polynomials = new HashMap<FundamentalName, double[]>(n);
188 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
189 lineNumber++;
190 if (definitionParser.parseDefinition(line, lineNumber, name)) {
191 polynomials.put(definitionParser.getParsedName(),
192 definitionParser.getParsedPolynomial());
193 }
194 }
195
196 final List<double[]> coefficients = new ArrayList<double[]>(n);
197 coefficients.add(getCoefficients(FundamentalName.L, polynomials, name));
198 coefficients.add(getCoefficients(FundamentalName.L_PRIME, polynomials, name));
199 coefficients.add(getCoefficients(FundamentalName.F, polynomials, name));
200 coefficients.add(getCoefficients(FundamentalName.D, polynomials, name));
201 coefficients.add(getCoefficients(FundamentalName.OMEGA, polynomials, name));
202 if (polynomials.containsKey(FundamentalName.L_ME)) {
203 // IERS conventions 2003 and later provide planetary nutation arguments
204 coefficients.add(getCoefficients(FundamentalName.L_ME, polynomials, name));
205 coefficients.add(getCoefficients(FundamentalName.L_VE, polynomials, name));
206 coefficients.add(getCoefficients(FundamentalName.L_E, polynomials, name));
207 coefficients.add(getCoefficients(FundamentalName.L_MA, polynomials, name));
208 coefficients.add(getCoefficients(FundamentalName.L_J, polynomials, name));
209 coefficients.add(getCoefficients(FundamentalName.L_SA, polynomials, name));
210 coefficients.add(getCoefficients(FundamentalName.L_U, polynomials, name));
211 coefficients.add(getCoefficients(FundamentalName.L_NE, polynomials, name));
212 coefficients.add(getCoefficients(FundamentalName.PA, polynomials, name));
213 } else {
214 // IERS conventions 1996 and earlier don't provide planetary nutation arguments
215 final double[] zero = new double[] {
216 0.0
217 };
218 while (coefficients.size() < n) {
219 coefficients.add(zero);
220 }
221 }
222
223 return coefficients;
224
225 } catch (IOException ioe) {
226 throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
227 }
228
229 }
230
231 /** Get the coefficients for a fundamental argument.
232 * @param argument fundamental argument
233 * @param polynomials map of the polynomials
234 * @param fileName name of the file from which the coefficients have been read
235 * @return polynomials coefficients (ordered from high degrees to low degrees)
236 * @exception OrekitException if the argument is not found
237 */
238 private static double[] getCoefficients(final FundamentalName argument,
239 final Map<FundamentalName, double[]> polynomials,
240 final String fileName)
241 throws OrekitException {
242 if (!polynomials.containsKey(argument)) {
243 throw new OrekitException(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, fileName);
244 }
245 return polynomials.get(argument);
246 }
247
248 /** Evaluate a polynomial.
249 * @param tc offset in Julian centuries
250 * @param coefficients polynomial coefficients (ordered from low degrees to high degrees)
251 * @return value of the polynomial
252 */
253 private double value(final double tc, final double[] coefficients) {
254 double value = 0;
255 for (int i = coefficients.length - 1; i >= 0; --i) {
256 value = coefficients[i] + tc * value;
257 }
258 return value;
259 }
260
261 /** Evaluate all fundamental arguments for the current date (Delaunay plus planetary).
262 * @param date current date
263 * @return all fundamental arguments for the current date (Delaunay plus planetary)
264 */
265 public BodiesElements evaluateAll(final AbsoluteDate date) {
266
267 final double tc = conventions.evaluateTC(date);
268 final double gamma = gmstFunction == null ?
269 Double.NaN : gmstFunction.value(date).getValue() + FastMath.PI;
270
271 return new BodiesElements(date, tc, gamma,
272 value(tc, lCoefficients), // mean anomaly of the Moon
273 value(tc, lPrimeCoefficients), // mean anomaly of the Sun
274 value(tc, fCoefficients), // L - Ω where L is the mean longitude of the Moon
275 value(tc, dCoefficients), // mean elongation of the Moon from the Sun
276 value(tc, omegaCoefficients), // mean longitude of the ascending node of the Moon
277 value(tc, lMeCoefficients), // mean Mercury longitude
278 value(tc, lVeCoefficients), // mean Venus longitude
279 value(tc, lECoefficients), // mean Earth longitude
280 value(tc, lMaCoefficients), // mean Mars longitude
281 value(tc, lJCoefficients), // mean Jupiter longitude
282 value(tc, lSaCoefficients), // mean Saturn longitude
283 value(tc, lUCoefficients), // mean Uranus longitude
284 value(tc, lNeCoefficients), // mean Neptune longitude
285 value(tc, paCoefficients)); // general accumulated precession in longitude
286
287 }
288
289 /** Evaluate a polynomial.
290 * @param tc offset in Julian centuries
291 * @param coefficients polynomial coefficients (ordered from low degrees to high degrees)
292 * @return value of the polynomial
293 */
294 private DerivativeStructure value(final DerivativeStructure tc, final double[] coefficients) {
295 DerivativeStructure value = tc.getField().getZero();
296 for (int i = coefficients.length - 1; i >= 0; --i) {
297 value = value.multiply(tc).add(coefficients[i]);
298 }
299 return value;
300 }
301
302 /** Evaluate all fundamental arguments for the current date (Delaunay plus planetary),
303 * including the first time derivative.
304 * @param date current date
305 * @return all fundamental arguments for the current date (Delaunay plus planetary),
306 * including the first time derivative
307 */
308 public FieldBodiesElements<DerivativeStructure> evaluateDerivative(final AbsoluteDate date) {
309
310 final DerivativeStructure tc = conventions.dsEvaluateTC(date);
311
312 return new FieldBodiesElements<DerivativeStructure>(date, tc,
313 gmstFunction.value(date).add(FastMath.PI),
314 value(tc, lCoefficients), // mean anomaly of the Moon
315 value(tc, lPrimeCoefficients), // mean anomaly of the Sun
316 value(tc, fCoefficients), // L - Ω where L is the mean longitude of the Moon
317 value(tc, dCoefficients), // mean elongation of the Moon from the Sun
318 value(tc, omegaCoefficients), // mean longitude of the ascending node of the Moon
319 value(tc, lMeCoefficients), // mean Mercury longitude
320 value(tc, lVeCoefficients), // mean Venus longitude
321 value(tc, lECoefficients), // mean Earth longitude
322 value(tc, lMaCoefficients), // mean Mars longitude
323 value(tc, lJCoefficients), // mean Jupiter longitude
324 value(tc, lSaCoefficients), // mean Saturn longitude
325 value(tc, lUCoefficients), // mean Uranus longitude
326 value(tc, lNeCoefficients), // mean Neptune longitude
327 value(tc, paCoefficients)); // general accumulated precession in longitude
328
329 }
330
331 /** Replace the instance with a data transfer object for serialization.
332 * <p>
333 * This intermediate class serializes only the frame key.
334 * </p>
335 * @return data transfer object that will be serialized
336 */
337 private Object writeReplace() {
338 return new DataTransferObject(conventions, timeScale,
339 Arrays.asList(lCoefficients, lPrimeCoefficients, fCoefficients,
340 dCoefficients, omegaCoefficients,
341 lMeCoefficients, lVeCoefficients, lECoefficients,
342 lMaCoefficients, lJCoefficients, lSaCoefficients,
343 lUCoefficients, lNeCoefficients, paCoefficients));
344 }
345
346 /** Internal class used only for serialization. */
347 private static class DataTransferObject implements Serializable {
348
349 /** Serializable UID. */
350 private static final long serialVersionUID = 20131209L;
351
352 /** IERS conventions to use. */
353 private final IERSConventions conventions;
354
355 /** Time scale for GMST computation. */
356 private final TimeScale timeScale;
357
358 /** All coefficients. */
359 private final List<double[]> coefficients;
360
361 /** Simple constructor.
362 * @param conventions IERS conventions to use
363 * @param timeScale time scale for GMST computation
364 * @param coefficients all coefficients
365 */
366 DataTransferObject(final IERSConventions conventions, final TimeScale timeScale,
367 final List<double[]> coefficients) {
368 this.conventions = conventions;
369 this.timeScale = timeScale;
370 this.coefficients = coefficients;
371 }
372
373 /** Replace the deserialized data transfer object with a {@link TIRFProvider}.
374 * @return replacement {@link TIRFProvider}
375 */
376 private Object readResolve() {
377 try {
378 // retrieve a managed frame
379 return new FundamentalNutationArguments(conventions, timeScale, coefficients);
380 } catch (OrekitException oe) {
381 throw new OrekitInternalError(oe);
382 }
383 }
384
385 }
386
387 /** Enumerate for the fundamental names. */
388 private enum FundamentalName {
389
390 /** Constant for Mean anomaly of the Moon. */
391 L() {
392 /** {@inheritDoc} */
393 public String getArgumentName() {
394 return "l";
395 }
396 },
397
398 /** Constant for Mean anomaly of the Sun. */
399 L_PRIME() {
400 /** {@inheritDoc} */
401 public String getArgumentName() {
402 return "l'";
403 }
404 },
405
406 /** Constant for L - Ω where L is the mean longitude of the Moon. */
407 F() {
408 /** {@inheritDoc} */
409 public String getArgumentName() {
410 return "F";
411 }
412 },
413
414 /** Constant for mean elongation of the Moon from the Sun. */
415 D() {
416 /** {@inheritDoc} */
417 public String getArgumentName() {
418 return "D";
419 }
420 },
421
422 /** Constant for longitude of the ascending node of the Moon. */
423 OMEGA() {
424 /** {@inheritDoc} */
425 public String getArgumentName() {
426 return "\u03a9";
427 }
428 },
429
430 /** Constant for mean Mercury longitude. */
431 L_ME() {
432 /** {@inheritDoc} */
433 public String getArgumentName() {
434 return "LMe";
435 }
436 },
437
438 /** Constant for mean Venus longitude. */
439 L_VE() {
440 /** {@inheritDoc} */
441 public String getArgumentName() {
442 return "LVe";
443 }
444 },
445
446 /** Constant for mean Earth longitude. */
447 L_E() {
448 /** {@inheritDoc} */
449 public String getArgumentName() {
450 return "LE";
451 }
452 },
453
454 /** Constant for mean Mars longitude. */
455 L_MA() {
456 /** {@inheritDoc} */
457 public String getArgumentName() {
458 return "LMa";
459 }
460 },
461
462 /** Constant for mean Jupiter longitude. */
463 L_J() {
464 /** {@inheritDoc} */
465 public String getArgumentName() {
466 return "LJ";
467 }
468 },
469
470 /** Constant for mean Saturn longitude. */
471 L_SA() {
472 /** {@inheritDoc} */
473 public String getArgumentName() {
474 return "LSa";
475 }
476 },
477
478 /** Constant for mean Uranus longitude. */
479 L_U() {
480 /** {@inheritDoc} */
481 public String getArgumentName() {
482 return "LU";
483 }
484 },
485
486 /** Constant for mean Neptune longitude. */
487 L_NE() {
488 /** {@inheritDoc} */
489 public String getArgumentName() {
490 return "LNe";
491 }
492 },
493
494 /** Constant for general accumulated precession in longitude. */
495 PA() {
496 /** {@inheritDoc} */
497 public String getArgumentName() {
498 return "pA";
499 }
500 };
501
502 /** Get the fundamental name.
503 * @return fundamental name
504 */
505 public abstract String getArgumentName();
506
507 }
508
509 /** Local parser for argument definition lines. */
510 private static class DefinitionParser {
511
512 /** Regular expression pattern for definitions. */
513 private final Pattern pattern;
514
515 /** Parser for polynomials. */
516 private PolynomialParser polynomialParser;
517
518 /** Last parsed fundamental name. */
519 private FundamentalName parsedName;
520
521 /** Last parsed polynomial. */
522 private double[] parsedPolynomial;
523
524 /** Simple constructor. */
525 DefinitionParser() {
526
527 // the luni-solar Delaunay arguments polynomial parts should read something like:
528 // F5 ≡ Ω = 125.04455501° − 6962890.5431″t + 7.4722″t² + 0.007702″t³ − 0.00005939″t⁴
529 // whereas the planetary arguments polynomial parts should read something like:
530 // F14 ≡ pA = 0.02438175 × t + 0.00000538691 × t²
531 final String unicodeIdenticalTo = "\u2261";
532
533 // pattern for the global line
534 final StringBuilder builder = new StringBuilder();
535 for (final FundamentalName fn : FundamentalName.values()) {
536 if (builder.length() > 0) {
537 builder.append('|');
538 }
539 builder.append(fn.getArgumentName());
540 }
541 final String fundamentalName = "\\p{Space}*((?:" + builder.toString() + ")+)";
542 pattern = Pattern.compile("\\p{Space}*F\\p{Digit}+\\p{Space}*" + unicodeIdenticalTo +
543 fundamentalName + "\\p{Space}*=\\p{Space}*(.*)");
544
545 polynomialParser = new PolynomialParser('t', PolynomialParser.Unit.NO_UNITS);
546
547 }
548
549 /** Parse a definition line.
550 * @param line line to parse
551 * @param lineNumber line number
552 * @param fileName name of the file
553 * @return true if a definition has been parsed
554 */
555 public boolean parseDefinition(final String line, final int lineNumber, final String fileName) {
556
557 parsedName = null;
558 parsedPolynomial = null;
559
560 final Matcher matcher = pattern.matcher(line);
561 if (matcher.matches()) {
562 for (FundamentalName fn : FundamentalName.values()) {
563 if (fn.getArgumentName().equals(matcher.group(1))) {
564 parsedName = fn;
565 }
566 }
567
568 // parse the polynomial
569 parsedPolynomial = polynomialParser.parse(matcher.group(2));
570
571 return true;
572
573 } else {
574 return false;
575 }
576
577 }
578
579 /** Get the last parsed fundamental name.
580 * @return last parsed fundamental name
581 */
582 public FundamentalName getParsedName() {
583 return parsedName;
584 }
585
586 /** Get the last parsed polynomial.
587 * @return last parsed polynomial
588 */
589 public double[] getParsedPolynomial() {
590 return parsedPolynomial.clone();
591 }
592
593 }
594
595 }