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.propagation;
18  
19  import java.io.Serializable;
20  import java.util.ArrayList;
21  import java.util.Collection;
22  import java.util.Collections;
23  import java.util.HashMap;
24  import java.util.List;
25  import java.util.Map;
26  
27  import org.apache.commons.math3.analysis.interpolation.HermiteInterpolator;
28  import org.apache.commons.math3.exception.DimensionMismatchException;
29  import org.orekit.attitudes.Attitude;
30  import org.orekit.attitudes.LofOffset;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.errors.OrekitMessages;
33  import org.orekit.frames.Frame;
34  import org.orekit.frames.LOFType;
35  import org.orekit.frames.Transform;
36  import org.orekit.orbits.Orbit;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.TimeInterpolable;
39  import org.orekit.time.TimeShiftable;
40  import org.orekit.time.TimeStamped;
41  import org.orekit.utils.PVCoordinates;
42  
43  
44  /** This class is the representation of a complete state holding orbit, attitude
45   * and mass information at a given date.
46   *
47   * <p>It contains an {@link Orbit orbital state} at a current
48   * {@link AbsoluteDate} both handled by an {@link Orbit}, plus the current
49   * mass and attitude. Orbit and state are guaranteed to be consistent in terms
50   * of date and reference frame. The spacecraft state may also contain additional
51   * states, which are simply named double arrays which can hold any user-defined
52   * data.
53   * </p>
54   * <p>
55   * The state can be slightly shifted to close dates. This shift is based on
56   * a simple keplerian model for orbit, a linear extrapolation for attitude
57   * taking the spin rate into account and no mass change. It is <em>not</em>
58   * intended as a replacement for proper orbit and attitude propagation but
59   * should be sufficient for either small time shifts or coarse accuracy.
60   * </p>
61   * <p>
62   * The instance <code>SpacecraftState</code> is guaranteed to be immutable.
63   * </p>
64   * @see org.orekit.propagation.numerical.NumericalPropagator
65   * @author Fabien Maussion
66   * @author V&eacute;ronique Pommier-Maurussane
67   * @author Luc Maisonobe
68   */
69  public class SpacecraftState
70      implements TimeStamped, TimeShiftable<SpacecraftState>, TimeInterpolable<SpacecraftState>, Serializable {
71  
72      /** Serializable UID. */
73      private static final long serialVersionUID = 20130407L;
74  
75      /** Default mass. */
76      private static final double DEFAULT_MASS = 1000.0;
77  
78      /** Orbital state. */
79      private final Orbit orbit;
80  
81      /** Attitude. */
82      private final Attitude attitude;
83  
84      /** Current mass (kg). */
85      private final double mass;
86  
87      /** Additional states. */
88      private final Map<String, double[]> additional;
89  
90      /** Build a spacecraft state from orbit only.
91       * <p>Attitude and mass are set to unspecified non-null arbitrary values.</p>
92       * @param orbit the orbit
93       * @exception OrekitException if default attitude cannot be computed
94       */
95      public SpacecraftState(final Orbit orbit)
96          throws OrekitException {
97          this(orbit,
98               new LofOffset(orbit.getFrame(), LOFType.VVLH).getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
99               DEFAULT_MASS, null);
100     }
101 
102     /** Build a spacecraft state from orbit and attitude provider.
103      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
104      * @param orbit the orbit
105      * @param attitude attitude
106      * @exception IllegalArgumentException if orbit and attitude dates
107      * or frames are not equal
108      */
109     public SpacecraftState(final Orbit orbit, final Attitude attitude)
110         throws IllegalArgumentException {
111         this(orbit, attitude, DEFAULT_MASS, null);
112     }
113 
114     /** Create a new instance from orbit and mass.
115      * <p>Attitude law is set to an unspecified default attitude.</p>
116      * @param orbit the orbit
117      * @param mass the mass (kg)
118      * @exception OrekitException if default attitude cannot be computed
119      */
120     public SpacecraftState(final Orbit orbit, final double mass)
121         throws OrekitException {
122         this(orbit,
123              new LofOffset(orbit.getFrame(), LOFType.VVLH).getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
124              mass, null);
125     }
126 
127     /** Build a spacecraft state from orbit, attitude provider and mass.
128      * @param orbit the orbit
129      * @param attitude attitude
130      * @param mass the mass (kg)
131      * @exception IllegalArgumentException if orbit and attitude dates
132      * or frames are not equal
133      */
134     public SpacecraftState(final Orbit orbit, final Attitude attitude, final double mass)
135         throws IllegalArgumentException {
136         this(orbit, attitude, mass, null);
137     }
138 
139     /** Build a spacecraft state from orbit only.
140      * <p>Attitude and mass are set to unspecified non-null arbitrary values.</p>
141      * @param orbit the orbit
142      * @param additional additional states
143      * @exception OrekitException if default attitude cannot be computed
144      */
145     public SpacecraftState(final Orbit orbit, final Map<String, double[]> additional)
146         throws OrekitException {
147         this(orbit,
148              new LofOffset(orbit.getFrame(), LOFType.VVLH).getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
149              DEFAULT_MASS, additional);
150     }
151 
152     /** Build a spacecraft state from orbit and attitude provider.
153      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
154      * @param orbit the orbit
155      * @param attitude attitude
156      * @param additional additional states
157      * @exception IllegalArgumentException if orbit and attitude dates
158      * or frames are not equal
159      */
160     public SpacecraftState(final Orbit orbit, final Attitude attitude, final Map<String, double[]> additional)
161         throws IllegalArgumentException {
162         this(orbit, attitude, DEFAULT_MASS, additional);
163     }
164 
165     /** Create a new instance from orbit and mass.
166      * <p>Attitude law is set to an unspecified default attitude.</p>
167      * @param orbit the orbit
168      * @param mass the mass (kg)
169      * @param additional additional states
170      * @exception OrekitException if default attitude cannot be computed
171      */
172     public SpacecraftState(final Orbit orbit, final double mass, final Map<String, double[]> additional)
173         throws OrekitException {
174         this(orbit,
175              new LofOffset(orbit.getFrame(), LOFType.VVLH).getAttitude(orbit, orbit.getDate(), orbit.getFrame()),
176              mass, additional);
177     }
178 
179     /** Build a spacecraft state from orbit, attitude provider and mass.
180      * @param orbit the orbit
181      * @param attitude attitude
182      * @param mass the mass (kg)
183      * @param additional additional states (may be null if no additional states are available)
184      * @exception IllegalArgumentException if orbit and attitude dates
185      * or frames are not equal
186      */
187     public SpacecraftState(final Orbit orbit, final Attitude attitude,
188                            final double mass, final Map<String, double[]> additional)
189         throws IllegalArgumentException {
190         checkConsistency(orbit, attitude);
191         this.orbit      = orbit;
192         this.attitude   = attitude;
193         this.mass       = mass;
194         if (additional == null) {
195             this.additional = Collections.emptyMap();
196         } else {
197             this.additional = new HashMap<String, double[]>(additional.size());
198             for (final Map.Entry<String, double[]> entry : additional.entrySet()) {
199                 this.additional.put(entry.getKey(), entry.getValue().clone());
200             }
201         }
202     }
203 
204     /** Add an additional state.
205      * <p>
206      * {@link SpacecraftState SpacecraftState} instances are immutable,
207      * so this method does <em>not</em> change the instance, but rather
208      * creates a new instance, which has the same orbit, attitude, mass
209      * and additional states as the original instance, except it also
210      * has the specified state. If the original instance already had an
211      * additional state with the same name, it will be overridden. If it
212      * did not have any additional state with that name, the new instance
213      * will have one more additional state than the original instance.
214      * </p>
215      * @param name name of the additional state
216      * @param value value of the additional state
217      * @return a new instance, with the additional state added
218      * @see #hasAdditionalState(String)
219      * @see #getAdditionalState(String)
220      * @see #getAdditionalStates()
221      */
222     public SpacecraftState addAdditionalState(final String name, final double ... value) {
223         final Map<String, double[]> newMap = new HashMap<String, double[]>(additional.size() + 1);
224         newMap.putAll(additional);
225         newMap.put(name, value.clone());
226         return new SpacecraftState(orbit, attitude, mass, newMap);
227     }
228 
229     /** Check orbit and attitude dates are equal.
230      * @param orbit the orbit
231      * @param attitude attitude
232      * @exception IllegalArgumentException if orbit and attitude dates
233      * are not equal
234      */
235     private static void checkConsistency(final Orbit orbit, final Attitude attitude)
236         throws IllegalArgumentException {
237         if (!orbit.getDate().equals(attitude.getDate())) {
238             throw OrekitException.createIllegalArgumentException(
239                   OrekitMessages.ORBIT_AND_ATTITUDE_DATES_MISMATCH,
240                   orbit.getDate(), attitude.getDate());
241         }
242         if (orbit.getFrame() != attitude.getReferenceFrame()) {
243             throw OrekitException.createIllegalArgumentException(
244                   OrekitMessages.FRAMES_MISMATCH,
245                   orbit.getFrame().getName(), attitude.getReferenceFrame().getName());
246         }
247     }
248 
249     /** Get a time-shifted state.
250      * <p>
251      * The state can be slightly shifted to close dates. This shift is based on
252      * a simple keplerian model for orbit, a linear extrapolation for attitude
253      * taking the spin rate into account and neither mass nor additional states
254      * changes. It is <em>not</em> intended as a replacement for proper orbit
255      * and attitude propagation but should be sufficient for small time shifts
256      * or coarse accuracy.
257      * </p>
258      * <p>
259      * As a rough order of magnitude, the following table shows the interpolation
260      * errors obtained between this simple shift method and an {@link
261      * org.orekit.propagation.analytical.EcksteinHechlerPropagator Eckstein-Heschler
262      * propagator} for an 800km altitude nearly circular polar Earth orbit with
263      * {@link org.orekit.attitudes.BodyCenterPointing body center pointing}. Beware
264      * that these results may be different for other orbits.
265      * </p>
266      * <table border="1" cellpadding="5">
267      * <tr bgcolor="#ccccff"><th>interpolation time (s)</th>
268      * <th>position error (m)</th><th>velocity error (m/s)</th>
269      * <th>attitude error (&deg;)</th></tr>
270      * <tr><td bgcolor="#eeeeff"> 60</td><td>  20</td><td>1</td><td>0.001</td></tr>
271      * <tr><td bgcolor="#eeeeff">120</td><td> 100</td><td>2</td><td>0.002</td></tr>
272      * <tr><td bgcolor="#eeeeff">300</td><td> 600</td><td>4</td><td>0.005</td></tr>
273      * <tr><td bgcolor="#eeeeff">600</td><td>2000</td><td>6</td><td>0.008</td></tr>
274      * <tr><td bgcolor="#eeeeff">900</td><td>4000</td><td>6</td><td>0.010</td></tr>
275      * </table>
276      * @param dt time shift in seconds
277      * @return a new state, shifted with respect to the instance (which is immutable)
278      * except for the mass which stay unchanged
279      */
280     public SpacecraftState shiftedBy(final double dt) {
281         return new SpacecraftState(orbit.shiftedBy(dt), attitude.shiftedBy(dt),
282                                    mass, additional);
283     }
284 
285     /** {@inheritDoc}
286      * <p>
287      * The additional states that are interpolated are the ones already present
288      * in the instance. The sample instances must therefore have at least the same
289      * additional states has the instance. They may have more additional states,
290      * but the extra ones will be ignored.
291      * </p>
292      * <p>
293      * As this implementation of interpolation is polynomial, it should be used only
294      * with small samples (about 10-20 points) in order to avoid <a
295      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
296      * and numerical problems (including NaN appearing).
297      * </p>
298      */
299     public SpacecraftState interpolate(final AbsoluteDate date,
300                                        final Collection<SpacecraftState> sample)
301         throws OrekitException {
302 
303         // prepare interpolators
304         final List<Orbit> orbits = new ArrayList<Orbit>(sample.size());
305         final List<Attitude> attitudes = new ArrayList<Attitude>(sample.size());
306         final HermiteInterpolator massInterpolator = new HermiteInterpolator();
307         final Map<String, HermiteInterpolator> additionalInterpolators =
308                 new HashMap<String, HermiteInterpolator>(additional.size());
309         for (final String name : additional.keySet()) {
310             additionalInterpolators.put(name, new HermiteInterpolator());
311         }
312 
313         // extract sample data
314         for (final SpacecraftState state : sample) {
315             final double deltaT = state.getDate().durationFrom(date);
316             orbits.add(state.getOrbit());
317             attitudes.add(state.getAttitude());
318             massInterpolator.addSamplePoint(deltaT,
319                                             new double[] {
320                                                 state.getMass()
321                                             });
322             for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
323                 entry.getValue().addSamplePoint(deltaT, state.getAdditionalState(entry.getKey()));
324             }
325         }
326 
327         // perform interpolations
328         final Orbit interpolatedOrbit       = orbit.interpolate(date, orbits);
329         final Attitude interpolatedAttitude = attitude.interpolate(date, attitudes);
330         final double interpolatedMass       = massInterpolator.value(0)[0];
331         final Map<String, double[]> interpolatedAdditional;
332         if (additional.isEmpty()) {
333             interpolatedAdditional = null;
334         } else {
335             interpolatedAdditional = new HashMap<String, double[]>(additional.size());
336             for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
337                 interpolatedAdditional.put(entry.getKey(), entry.getValue().value(0));
338             }
339         }
340 
341         // create the complete interpolated state
342         return new SpacecraftState(interpolatedOrbit, interpolatedAttitude,
343                                    interpolatedMass, interpolatedAdditional);
344 
345     }
346 
347     /** Gets the current orbit.
348      * @return the orbit
349      */
350     public Orbit getOrbit() {
351         return orbit;
352     }
353 
354     /** Get the date.
355      * @return date
356      */
357     public AbsoluteDate getDate() {
358         return orbit.getDate();
359     }
360 
361     /** Get the inertial frame.
362      * @return the frame
363      */
364     public Frame getFrame() {
365         return orbit.getFrame();
366     }
367 
368     /** Check if an additional state is available.
369      * @param name name of the additional state
370      * @return true if the additional state is available
371      * @see #addAdditionalState(String, double[])
372      * @see #getAdditionalState(String)
373      * @see #getAdditionalStates()
374      */
375     public boolean hasAdditionalState(final String name) {
376         return additional.containsKey(name);
377     }
378 
379     /** Check if two instances have the same set of additional states available.
380      * <p>
381      * Only the names and dimensions of the additional states are compared,
382      * not their values.
383      * </p>
384      * @param state state to compare to instance
385      * @exception OrekitException if either instance or state supports an additional
386      * state not supported by the other one
387      * @exception DimensionMismatchException if an additional state does not have
388      * the same dimension in both states
389      */
390     public void ensureCompatibleAdditionalStates(final SpacecraftState state)
391         throws OrekitException, DimensionMismatchException {
392 
393         // check instance additional states is a subset of the other one
394         for (final Map.Entry<String, double[]> entry : additional.entrySet()) {
395             final double[] other = state.additional.get(entry.getKey());
396             if (other == null) {
397                 throw new OrekitException(OrekitMessages.UNKNOWN_ADDITIONAL_STATE,
398                                           entry.getKey());
399             }
400             if (other.length != entry.getValue().length) {
401                 throw new DimensionMismatchException(other.length, entry.getValue().length);
402             }
403         }
404 
405         if (state.additional.size() > additional.size()) {
406             // the other state has more additional states
407             for (final String name : state.additional.keySet()) {
408                 if (!additional.containsKey(name)) {
409                     throw new OrekitException(OrekitMessages.UNKNOWN_ADDITIONAL_STATE,
410                                               name);
411                 }
412             }
413         }
414 
415     }
416 
417     /** Get an additional state.
418      * @param name name of the additional state
419      * @return value of the additional state
420      * @exception OrekitException if no additional state with that name exists
421      * @see #addAdditionalState(String, double[])
422      * @see #hasAdditionalState(String)
423      * @see #getAdditionalStates()
424      */
425     public double[] getAdditionalState(final String name) throws OrekitException {
426         if (!additional.containsKey(name)) {
427             throw new OrekitException(OrekitMessages.UNKNOWN_ADDITIONAL_STATE, name);
428         }
429         return additional.get(name).clone();
430     }
431 
432     /** Get an unmodifiable map of additional states.
433      * @return unmodifiable map of additional states
434      * @see #addAdditionalState(String, double[])
435      * @see #hasAdditionalState(String)
436      * @see #getAdditionalState(String)
437      */
438     public Map<String, double[]> getAdditionalStates() {
439         return Collections.unmodifiableMap(additional);
440     }
441 
442     /** Compute the transform from orbite/attitude reference frame to spacecraft frame.
443      * <p>The spacecraft frame origin is at the point defined by the orbit,
444      * and its orientation is defined by the attitude.</p>
445      * @return transform from specified frame to current spacecraft frame
446      */
447     public Transform toTransform() {
448         final AbsoluteDate date = orbit.getDate();
449         return new Transform(date,
450                              new Transform(date, orbit.getPVCoordinates().negate()),
451                              new Transform(date, attitude.getOrientation()));
452     }
453 
454     /** Get the central attraction coefficient.
455      * @return mu central attraction coefficient (m^3/s^2)
456      */
457     public double getMu() {
458         return orbit.getMu();
459     }
460 
461     /** Get the keplerian period.
462      * <p>The keplerian period is computed directly from semi major axis
463      * and central acceleration constant.</p>
464      * @return keplerian period in seconds
465      */
466     public double getKeplerianPeriod() {
467         return orbit.getKeplerianPeriod();
468     }
469 
470     /** Get the keplerian mean motion.
471      * <p>The keplerian mean motion is computed directly from semi major axis
472      * and central acceleration constant.</p>
473      * @return keplerian mean motion in radians per second
474      */
475     public double getKeplerianMeanMotion() {
476         return orbit.getKeplerianMeanMotion();
477     }
478 
479     /** Get the semi-major axis.
480      * @return semi-major axis (m)
481      */
482     public double getA() {
483         return orbit.getA();
484     }
485 
486     /** Get the first component of the eccentricity vector (as per equinoctial parameters).
487      * @return e cos(&omega; + &Omega;), first component of eccentricity vector
488      * @see #getE()
489      */
490     public double getEquinoctialEx() {
491         return orbit.getEquinoctialEx();
492     }
493 
494     /** Get the second component of the eccentricity vector (as per equinoctial parameters).
495      * @return e sin(&omega; + &Omega;), second component of the eccentricity vector
496      * @see #getE()
497      */
498     public double getEquinoctialEy() {
499         return orbit.getEquinoctialEy();
500     }
501 
502     /** Get the first component of the inclination vector (as per equinoctial parameters).
503      * @return tan(i/2) cos(&Omega;), first component of the inclination vector
504      * @see #getI()
505      */
506     public double getHx() {
507         return orbit.getHx();
508     }
509 
510     /** Get the second component of the inclination vector (as per equinoctial parameters).
511      * @return tan(i/2) sin(&Omega;), second component of the inclination vector
512      * @see #getI()
513      */
514     public double getHy() {
515         return orbit.getHy();
516     }
517 
518     /** Get the true latitude argument (as per equinoctial parameters).
519      * @return v + &omega; + &Omega; true latitude argument (rad)
520      * @see #getLE()
521      * @see #getLM()
522      */
523     public double getLv() {
524         return orbit.getLv();
525     }
526 
527     /** Get the eccentric latitude argument (as per equinoctial parameters).
528      * @return E + &omega; + &Omega; eccentric latitude argument (rad)
529      * @see #getLv()
530      * @see #getLM()
531      */
532     public double getLE() {
533         return orbit.getLE();
534     }
535 
536     /** Get the mean latitude argument (as per equinoctial parameters).
537      * @return M + &omega; + &Omega; mean latitude argument (rad)
538      * @see #getLv()
539      * @see #getLE()
540      */
541     public double getLM() {
542         return orbit.getLM();
543     }
544 
545     // Additional orbital elements
546 
547     /** Get the eccentricity.
548      * @return eccentricity
549      * @see #getEquinoctialEx()
550      * @see #getEquinoctialEy()
551      */
552     public double getE() {
553         return orbit.getE();
554     }
555 
556     /** Get the inclination.
557      * @return inclination (rad)
558      * @see #getHx()
559      * @see #getHy()
560      */
561     public double getI() {
562         return orbit.getI();
563     }
564 
565     /** Get the {@link PVCoordinates} in orbit definition frame.
566      * Compute the position and velocity of the satellite. This method caches its
567      * results, and recompute them only when the method is called with a new value
568      * for mu. The result is provided as a reference to the internally cached
569      * {@link PVCoordinates}, so the caller is responsible to copy it in a separate
570      * {@link PVCoordinates} if it needs to keep the value for a while.
571      * @return pvCoordinates in orbit definition frame
572      */
573     public PVCoordinates getPVCoordinates() {
574         return orbit.getPVCoordinates();
575     }
576 
577     /** Get the {@link PVCoordinates} in given output frame.
578      * Compute the position and velocity of the satellite. This method caches its
579      * results, and recompute them only when the method is called with a new value
580      * for mu. The result is provided as a reference to the internally cached
581      * {@link PVCoordinates}, so the caller is responsible to copy it in a separate
582      * {@link PVCoordinates} if it needs to keep the value for a while.
583      * @param outputFrame frame in which coordinates should be defined
584      * @return pvCoordinates in orbit definition frame
585      * @exception OrekitException if the transformation between frames cannot be computed
586      */
587     public PVCoordinates getPVCoordinates(final Frame outputFrame)
588         throws OrekitException {
589         return orbit.getPVCoordinates(outputFrame);
590     }
591 
592     /** Get the attitude.
593      * @return the attitude.
594      */
595     public Attitude getAttitude() {
596         return attitude;
597     }
598 
599     /** Gets the current mass.
600      * @return the mass (kg)
601      */
602     public double getMass() {
603         return mass;
604     }
605 
606 }