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.hipparchus.analysis.differentiation.DerivativeStructure;
33  import org.hipparchus.exception.DummyLocalizable;
34  import org.hipparchus.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 }