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