1   /* Copyright 2002-2015 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.frames;
18  
19  import java.io.Serializable;
20  import java.util.ArrayList;
21  import java.util.Arrays;
22  import java.util.Collection;
23  import java.util.List;
24  
25  import org.apache.commons.math3.RealFieldElement;
26  import org.apache.commons.math3.geometry.euclidean.threed.FieldRotation;
27  import org.apache.commons.math3.geometry.euclidean.threed.FieldVector3D;
28  import org.apache.commons.math3.geometry.euclidean.threed.Line;
29  import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
30  import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
31  import org.orekit.errors.OrekitException;
32  import org.orekit.time.AbsoluteDate;
33  import org.orekit.time.TimeInterpolable;
34  import org.orekit.time.TimeShiftable;
35  import org.orekit.time.TimeStamped;
36  import org.orekit.utils.AngularCoordinates;
37  import org.orekit.utils.AngularDerivativesFilter;
38  import org.orekit.utils.CartesianDerivativesFilter;
39  import org.orekit.utils.FieldPVCoordinates;
40  import org.orekit.utils.PVCoordinates;
41  import org.orekit.utils.TimeStampedAngularCoordinates;
42  import org.orekit.utils.TimeStampedFieldPVCoordinates;
43  import org.orekit.utils.TimeStampedPVCoordinates;
44  
45  
46  /** Transformation class in three dimensional space.
47   *
48   * <p>This class represents the transformation engine between {@link Frame frames}.
49   * It is used both to define the relationship between each frame and its
50   * parent frame and to gather all individual transforms into one
51   * operation when converting between frames far away from each other.</p>
52   * <p> The convention used in OREKIT is vectorial transformation. It means
53   * that a transformation is defined as a transform to apply to the
54   * coordinates of a vector expressed in the old frame to obtain the
55   * same vector expressed in the new frame.<p>
56   *
57   * <p>Instances of this class are guaranteed to be immutable.</p>
58   *
59   *  <h5> Example </h5>
60   *
61   * <pre>
62   *
63   * 1 ) Example of translation from R<sub>A</sub> to R<sub>B</sub>:
64   * We want to transform the {@link PVCoordinates} PV<sub>A</sub> to PV<sub>B</sub>.
65   *
66   * With :  PV<sub>A</sub> = ({1, 0, 0}, {2, 0, 0}, {3, 0, 0});
67   * and  :  PV<sub>B</sub> = ({0, 0, 0}, {0, 0, 0}, {0, 0, 0});
68   *
69   * The transform to apply then is defined as follows :
70   *
71   * Vector3D translation  = new Vector3D(-1, 0, 0);
72   * Vector3D velocity     = new Vector3D(-2, 0, 0);
73   * Vector3D acceleration = new Vector3D(-3, 0, 0);
74   *
75   * Transform R1toR2 = new Transform(date, translation, velocity, acceleration);
76   *
77   * PV<sub>B</sub> = R1toR2.transformPVCoordinates(PV<sub>A</sub>);
78   *
79   *
80   * 2 ) Example of rotation from R<sub>A</sub> to R<sub>B</sub>:
81   * We want to transform the {@link PVCoordinates} PV<sub>A</sub> to PV<sub>B</sub>.
82   *
83   * With :  PV<sub>A</sub> = ({1, 0, 0}, {1, 0, 0});
84   * and  :  PV<sub>B</sub> = ({0, 1, 0}, {-2, 1, 0});
85   *
86   * The transform to apply then is defined as follows :
87   *
88   * Rotation rotation = new Rotation(Vector3D.PLUS_K, FastMath.PI / 2);
89   * Vector3D rotationRate = new Vector3D(0, 0, -2);
90   *
91   * Transform R1toR2 = new Transform(rotation, rotationRate);
92   *
93   * PV<sub>B</sub> = R1toR2.transformPVCoordinates(PV<sub>A</sub>);
94   *
95   * </pre>
96   *
97   * @author Luc Maisonobe
98   * @author Fabien Maussion
99   */
100 public class Transform
101     implements TimeStamped, TimeShiftable<Transform>, TimeInterpolable<Transform>, Serializable {
102 
103     /** Identity transform. */
104     public static final Transform IDENTITY = new IdentityTransform();
105 
106     /** Serializable UID. */
107     private static final long serialVersionUID = 210140410L;
108 
109     /** Date of the transform. */
110     private final AbsoluteDate date;
111 
112     /** Cartesian coordinates of the target frame with respect to the original frame. */
113     private final PVCoordinates cartesian;
114 
115     /** Angular coordinates of the target frame with respect to the original frame. */
116     private final AngularCoordinates angular;
117 
118     /** Build a transform from its primitive operations.
119      * @param date date of the transform
120      * @param cartesian Cartesian coordinates of the target frame with respect to the original frame
121      * @param angular angular coordinates of the target frame with respect to the original frame
122      */
123     private Transform(final AbsoluteDate date,
124                       final PVCoordinates cartesian, final AngularCoordinates angular) {
125         this.date      = date;
126         this.cartesian = cartesian;
127         this.angular   = angular;
128     }
129 
130     /** Build a translation transform.
131      * @param date date of the transform
132      * @param translation translation to apply (i.e. coordinates of
133      * the transformed origin, or coordinates of the origin of the
134      * old frame in the new frame)
135      */
136     public Transform(final AbsoluteDate date, final Vector3D translation) {
137         this(date,
138              new PVCoordinates(translation, Vector3D.ZERO, Vector3D.ZERO),
139              AngularCoordinates.IDENTITY);
140     }
141 
142     /** Build a rotation transform.
143      * @param date date of the transform
144      * @param rotation rotation to apply ( i.e. rotation to apply to the
145      * coordinates of a vector expressed in the old frame to obtain the
146      * same vector expressed in the new frame )
147      */
148     public Transform(final AbsoluteDate date, final Rotation rotation) {
149         this(date,
150              PVCoordinates.ZERO,
151              new AngularCoordinates(rotation, Vector3D.ZERO));
152     }
153 
154     /** Build a translation transform, with its first time derivative.
155      * @param date date of the transform
156      * @param translation translation to apply (i.e. coordinates of
157      * the transformed origin, or coordinates of the origin of the
158      * old frame in the new frame)
159      * @param velocity the velocity of the translation (i.e. origin
160      * of the old frame velocity in the new frame)
161      */
162     public Transform(final AbsoluteDate date, final Vector3D translation,
163                      final Vector3D velocity) {
164         this(date,
165              new PVCoordinates(translation, velocity, Vector3D.ZERO),
166              AngularCoordinates.IDENTITY);
167     }
168 
169     /** Build a translation transform, with its first and second time derivatives.
170      * @param date date of the transform
171      * @param translation translation to apply (i.e. coordinates of
172      * the transformed origin, or coordinates of the origin of the
173      * old frame in the new frame)
174      * @param velocity the velocity of the translation (i.e. origin
175      * of the old frame velocity in the new frame)
176      * @param acceleration the acceleration of the translation (i.e. origin
177      * of the old frame acceleration in the new frame)
178      */
179     public Transform(final AbsoluteDate date, final Vector3D translation,
180                      final Vector3D velocity, final Vector3D acceleration) {
181         this(date,
182              new PVCoordinates(translation, velocity, acceleration),
183              AngularCoordinates.IDENTITY);
184     }
185 
186     /** Build a translation transform, with its first time derivative.
187      * @param date date of the transform
188      * @param cartesian cartesian part of the transformation to apply (i.e. coordinates of
189      * the transformed origin, or coordinates of the origin of the
190      * old frame in the new frame, with their derivatives)
191      */
192     public Transform(final AbsoluteDate date, final PVCoordinates cartesian) {
193         this(date,
194              cartesian,
195              AngularCoordinates.IDENTITY);
196     }
197 
198     /** Build a rotation transform.
199      * @param date date of the transform
200      * @param rotation rotation to apply ( i.e. rotation to apply to the
201      * coordinates of a vector expressed in the old frame to obtain the
202      * same vector expressed in the new frame )
203      * @param rotationRate the axis of the instant rotation
204      * expressed in the new frame. (norm representing angular rate)
205      */
206     public Transform(final AbsoluteDate date, final Rotation rotation, final Vector3D rotationRate) {
207         this(date,
208              PVCoordinates.ZERO,
209              new AngularCoordinates(rotation, rotationRate, Vector3D.ZERO));
210     }
211 
212     /** Build a rotation transform.
213      * @param date date of the transform
214      * @param rotation rotation to apply ( i.e. rotation to apply to the
215      * coordinates of a vector expressed in the old frame to obtain the
216      * same vector expressed in the new frame )
217      * @param rotationRate the axis of the instant rotation
218      * @param rotationAcceleration the axis of the instant rotation
219      * expressed in the new frame. (norm representing angular rate)
220      */
221     public Transform(final AbsoluteDate date, final Rotation rotation, final Vector3D rotationRate,
222                      final Vector3D rotationAcceleration) {
223         this(date,
224              PVCoordinates.ZERO,
225              new AngularCoordinates(rotation, rotationRate, rotationAcceleration));
226     }
227 
228     /** Build a rotation transform.
229      * @param date date of the transform
230      * @param angular angular part of the transformation to apply (i.e. rotation to
231      * apply to the coordinates of a vector expressed in the old frame to obtain the
232      * same vector expressed in the new frame, with its rotation rate)
233      */
234     public Transform(final AbsoluteDate date, final AngularCoordinates angular) {
235         this(date, PVCoordinates.ZERO, angular);
236     }
237 
238     /** Build a transform by combining two existing ones.
239      * <p>
240      * Note that the dates of the two existing transformed are <em>ignored</em>,
241      * and the combined transform date is set to the date supplied in this constructor
242      * without any attempt to shift the raw transforms. This is a design choice allowing
243      * user full control of the combination.
244      * </p>
245      * @param date date of the transform
246      * @param first first transform applied
247      * @param second second transform applied
248      */
249     public Transform(final AbsoluteDate date, final Transform first, final Transform second) {
250         this(date,
251              new PVCoordinates(compositeTranslation(first, second),
252                                compositeVelocity(first, second),
253                                compositeAcceleration(first, second)),
254              new AngularCoordinates(compositeRotation(first, second),
255                                     compositeRotationRate(first, second),
256                                     compositeRotationAcceleration(first, second)));
257     }
258 
259     /** Compute a composite translation.
260      * @param first first applied transform
261      * @param second second applied transform
262      * @return translation part of the composite transform
263      */
264     private static Vector3D compositeTranslation(final Transform first, final Transform second) {
265 
266         final Vector3D p1 = first.cartesian.getPosition();
267         final Rotation r1 = first.angular.getRotation();
268         final Vector3D p2 = second.cartesian.getPosition();
269 
270         return p1.add(r1.applyInverseTo(p2));
271 
272     }
273 
274     /** Compute a composite velocity.
275      * @param first first applied transform
276      * @param second second applied transform
277      * @return velocity part of the composite transform
278      */
279     private static Vector3D compositeVelocity(final Transform first, final Transform second) {
280 
281         final Vector3D v1 = first.cartesian.getVelocity();
282         final Rotation r1 = first.angular.getRotation();
283         final Vector3D o1 = first.angular.getRotationRate();
284         final Vector3D p2 = second.cartesian.getPosition();
285         final Vector3D v2 = second.cartesian.getVelocity();
286 
287         final Vector3D crossP = Vector3D.crossProduct(o1, p2);
288 
289         return v1.add(r1.applyInverseTo(v2.add(crossP)));
290 
291     }
292 
293     /** Compute a composite acceleration.
294      * @param first first applied transform
295      * @param second second applied transform
296      * @return acceleration part of the composite transform
297      */
298     private static Vector3D compositeAcceleration(final Transform first, final Transform second) {
299 
300         final Vector3D a1    = first.cartesian.getAcceleration();
301         final Rotation r1    = first.angular.getRotation();
302         final Vector3D o1    = first.angular.getRotationRate();
303         final Vector3D oDot1 = first.angular.getRotationAcceleration();
304         final Vector3D p2    = second.cartesian.getPosition();
305         final Vector3D v2    = second.cartesian.getVelocity();
306         final Vector3D a2    = second.cartesian.getAcceleration();
307 
308         final Vector3D crossCrossP = Vector3D.crossProduct(o1,    Vector3D.crossProduct(o1, p2));
309         final Vector3D crossV      = Vector3D.crossProduct(o1,    v2);
310         final Vector3D crossDotP   = Vector3D.crossProduct(oDot1, p2);
311 
312         return a1.add(r1.applyInverseTo(new Vector3D(1, a2, 2, crossV, 1, crossCrossP, 1, crossDotP)));
313 
314     }
315 
316     /** Compute a composite rotation.
317      * @param first first applied transform
318      * @param second second applied transform
319      * @return rotation part of the composite transform
320      */
321     private static Rotation compositeRotation(final Transform first, final Transform second) {
322 
323         final Rotation r1 = first.angular.getRotation();
324         final Rotation r2 = second.angular.getRotation();
325 
326         return r2.applyTo(r1);
327 
328     }
329 
330     /** Compute a composite rotation rate.
331      * @param first first applied transform
332      * @param second second applied transform
333      * @return rotation rate part of the composite transform
334      */
335     private static Vector3D compositeRotationRate(final Transform first, final Transform second) {
336 
337         final Vector3D o1 = first.angular.getRotationRate();
338         final Rotation r2 = second.angular.getRotation();
339         final Vector3D o2 = second.angular.getRotationRate();
340 
341         return o2.add(r2.applyTo(o1));
342 
343     }
344 
345     /** Compute a composite rotation acceleration.
346      * @param first first applied transform
347      * @param second second applied transform
348      * @return rotation acceleration part of the composite transform
349      */
350     private static Vector3D compositeRotationAcceleration(final Transform first, final Transform second) {
351 
352         final Vector3D o1    = first.angular.getRotationRate();
353         final Vector3D oDot1 = first.angular.getRotationAcceleration();
354         final Rotation r2    = second.angular.getRotation();
355         final Vector3D o2    = second.angular.getRotationRate();
356         final Vector3D oDot2 = second.angular.getRotationAcceleration();
357 
358         return new Vector3D( 1, oDot2,
359                              1, r2.applyTo(oDot1),
360                             -1, Vector3D.crossProduct(o2, r2.applyTo(o1)));
361 
362     }
363 
364     /** {@inheritDoc} */
365     public AbsoluteDate getDate() {
366         return date;
367     }
368 
369     /** {@inheritDoc} */
370     public Transform shiftedBy(final double dt) {
371         return new Transform(date.shiftedBy(dt), cartesian.shiftedBy(dt), angular.shiftedBy(dt));
372     };
373 
374     /** {@inheritDoc}
375      * <p>
376      * Calling this method is equivalent to call {@link #interpolate(AbsoluteDate,
377      * CartesianDerivativesFilter, AngularDerivativesFilter, Collection)} with {@code cFilter}
378      * set to {@link CartesianDerivativesFilter#USE_PVA} and {@code aFilter} set to
379      * {@link AngularDerivativesFilter#USE_RRA}
380      * set to true.
381      * </p>
382      * @exception OrekitException if the number of point is too small for interpolating
383      */
384     public Transform interpolate(final AbsoluteDate interpolationDate, final Collection<Transform> sample)
385         throws OrekitException {
386         return interpolate(interpolationDate,
387                            CartesianDerivativesFilter.USE_PVA, AngularDerivativesFilter.USE_RRA,
388                            sample);
389     }
390 
391     /** Interpolate a transform from a sample set of existing transforms.
392      * <p>
393      * Note that even if first time derivatives (velocities and rotation rates)
394      * from sample can be ignored, the interpolated instance always includes
395      * interpolated derivatives. This feature can be used explicitly to
396      * compute these derivatives when it would be too complex to compute them
397      * from an analytical formula: just compute a few sample points from the
398      * explicit formula and set the derivatives to zero in these sample points,
399      * then use interpolation to add derivatives consistent with the positions
400      * and rotations.
401      * </p>
402      * <p>
403      * As this implementation of interpolation is polynomial, it should be used only
404      * with small samples (about 10-20 points) in order to avoid <a
405      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
406      * and numerical problems (including NaN appearing).
407      * </p>
408      * @param date interpolation date
409      * @param useVelocities if true, use sample transforms velocities,
410      * otherwise ignore them and use only positions
411      * @param useRotationRates if true, use sample points rotation rates,
412      * otherwise ignore them and use only rotations
413      * @param sample sample points on which interpolation should be done
414      * @return a new instance, interpolated at specified date
415      * @exception OrekitException if the number of point is too small for interpolating
416      * @deprecated as of 7.0, replaced with {@link #interpolate(AbsoluteDate, CartesianDerivativesFilter, AngularDerivativesFilter, Collection)}
417      */
418     @Deprecated
419     public static Transform interpolate(final AbsoluteDate date,
420                                         final boolean useVelocities, final boolean useRotationRates,
421                                         final Collection<Transform> sample)
422         throws OrekitException {
423         return interpolate(date,
424                            useVelocities    ? CartesianDerivativesFilter.USE_PV : CartesianDerivativesFilter.USE_P,
425                            useRotationRates ? AngularDerivativesFilter.USE_RR   : AngularDerivativesFilter.USE_R,
426                            sample);
427     }
428 
429     /** Interpolate a transform from a sample set of existing transforms.
430      * <p>
431      * Note that even if first time derivatives (velocities and rotation rates)
432      * from sample can be ignored, the interpolated instance always includes
433      * interpolated derivatives. This feature can be used explicitly to
434      * compute these derivatives when it would be too complex to compute them
435      * from an analytical formula: just compute a few sample points from the
436      * explicit formula and set the derivatives to zero in these sample points,
437      * then use interpolation to add derivatives consistent with the positions
438      * and rotations.
439      * </p>
440      * <p>
441      * As this implementation of interpolation is polynomial, it should be used only
442      * with small samples (about 10-20 points) in order to avoid <a
443      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
444      * and numerical problems (including NaN appearing).
445      * </p>
446      * @param date interpolation date
447      * @param cFilter filter for derivatives from the sample to use in interpolation
448      * @param aFilter filter for derivatives from the sample to use in interpolation
449      * @param sample sample points on which interpolation should be done
450      * @return a new instance, interpolated at specified date
451      * @exception OrekitException if the number of point is too small for interpolating
452      * @since 7.0
453      */
454     public static Transform interpolate(final AbsoluteDate date,
455                                         final CartesianDerivativesFilter cFilter,
456                                         final AngularDerivativesFilter aFilter,
457                                         final Collection<Transform> sample)
458         throws OrekitException {
459         final List<TimeStampedPVCoordinates>      datedPV = new ArrayList<TimeStampedPVCoordinates>(sample.size());
460         final List<TimeStampedAngularCoordinates> datedAC = new ArrayList<TimeStampedAngularCoordinates>(sample.size());
461         for (final Transform t : sample) {
462             datedPV.add(new TimeStampedPVCoordinates(t.getDate(), t.getTranslation(), t.getVelocity(), t.getAcceleration()));
463             datedAC.add(new TimeStampedAngularCoordinates(t.getDate(), t.getRotation(), t.getRotationRate(), t.getRotationAcceleration()));
464         }
465         final TimeStampedPVCoordinates      interpolatedPV = TimeStampedPVCoordinates.interpolate(date, cFilter, datedPV);
466         final TimeStampedAngularCoordinates interpolatedAC = TimeStampedAngularCoordinates.interpolate(date, aFilter, datedAC);
467         return new Transform(date, interpolatedPV, interpolatedAC);
468     }
469 
470     /** Get the inverse transform of the instance.
471      * @return inverse transform of the instance
472      */
473     public Transform getInverse() {
474 
475         final Rotation r    = angular.getRotation();
476         final Vector3D o    = angular.getRotationRate();
477         final Vector3D oDot = angular.getRotationAcceleration();
478         final Vector3D rp   = r.applyTo(cartesian.getPosition());
479         final Vector3D rv   = r.applyTo(cartesian.getVelocity());
480         final Vector3D ra   = r.applyTo(cartesian.getAcceleration());
481 
482         final Vector3D pInv        = rp.negate();
483         final Vector3D crossP      = Vector3D.crossProduct(o, rp);
484         final Vector3D vInv        = crossP.subtract(rv);
485         final Vector3D crossV      = Vector3D.crossProduct(o, rv);
486         final Vector3D crossDotP   = Vector3D.crossProduct(oDot, rp);
487         final Vector3D crossCrossP = Vector3D.crossProduct(o, crossP);
488         final Vector3D aInv        = new Vector3D(-1, ra,
489                                                    2, crossV,
490                                                    1, crossDotP,
491                                                   -1, crossCrossP);
492 
493         return new Transform(date, new PVCoordinates(pInv, vInv, aInv), angular.revert());
494 
495     }
496 
497     /** Get a frozen transform.
498      * <p>
499      * This method creates a copy of the instance but frozen in time,
500      * i.e. with velocity, acceleration and rotation rate forced to zero.
501      * </p>
502      * @return a new transform, without any time-dependent parts
503      */
504     public Transform freeze() {
505         return new Transform(date,
506                              new PVCoordinates(cartesian.getPosition(), Vector3D.ZERO, Vector3D.ZERO),
507                              new AngularCoordinates(angular.getRotation(), Vector3D.ZERO, Vector3D.ZERO));
508     }
509 
510     /** Transform a position vector (including translation effects).
511      * @param position vector to transform
512      * @return transformed position
513      */
514     public Vector3D transformPosition(final Vector3D position) {
515         return angular.getRotation().applyTo(cartesian.getPosition().add(position));
516     }
517 
518     /** Transform a position vector (including translation effects).
519      * @param position vector to transform
520      * @param <T> the type of the field elements
521      * @return transformed position
522      */
523     public <T extends RealFieldElement<T>> FieldVector3D<T> transformPosition(final FieldVector3D<T> position) {
524         return FieldRotation.applyTo(angular.getRotation(), position.add(cartesian.getPosition()));
525     }
526 
527     /** Transform a vector (ignoring translation effects).
528      * @param vector vector to transform
529      * @return transformed vector
530      */
531     public Vector3D transformVector(final Vector3D vector) {
532         return angular.getRotation().applyTo(vector);
533     }
534 
535     /** Transform a vector (ignoring translation effects).
536      * @param vector vector to transform
537      * @param <T> the type of the field elements
538      * @return transformed vector
539      */
540     public <T extends RealFieldElement<T>> FieldVector3D<T> transformVector(final FieldVector3D<T> vector) {
541         return FieldRotation.applyTo(angular.getRotation(), vector);
542     }
543 
544     /** Transform a line.
545      * @param line to transform
546      * @return transformed line
547      */
548     public Line transformLine(final Line line) {
549         final Vector3D transformedP0 = transformPosition(line.getOrigin());
550         final Vector3D transformedP1 = transformPosition(line.pointAt(1.0e6));
551         return new Line(transformedP0, transformedP1, 1.0e-10);
552     }
553 
554     /** Transform {@link PVCoordinates} including kinematic effects.
555      * @param pva the position-velocity-acceleration triplet to transform.
556      * @return transformed position-velocity-acceleration
557      */
558     public PVCoordinates transformPVCoordinates(final PVCoordinates pva) {
559         return angular.applyTo(new PVCoordinates(1, pva, 1, cartesian));
560     }
561 
562     /** Transform {@link TimeStampedPVCoordinates} including kinematic effects.
563      * <p>
564      * In order to allow the user more flexibility, this method does <em>not</em> check for
565      * consistency between the transform {@link #getDate() date} and the time-stamped
566      * position-velocity {@link TimeStampedPVCoordinates#getDate() date}. The returned
567      * value will always have the same {@link TimeStampedPVCoordinates#getDate() date} as
568      * the input argument, regardless of the instance {@link #getDate() date}.
569      * </p>
570      * @param pv time-stamped  position-velocity to transform.
571      * @return transformed time-stamped position-velocity
572      * @since 7.0
573      */
574     public TimeStampedPVCoordinates transformPVCoordinates(final TimeStampedPVCoordinates pv) {
575         return angular.applyTo(new TimeStampedPVCoordinates(pv.getDate(), 1, pv, 1, cartesian));
576     }
577 
578     /** Transform {@link FieldPVCoordinates} including kinematic effects.
579      * @param pv position-velocity to transform.
580      * @param <T> type of the field elements
581      * @return transformed position-velocity
582      */
583     public <T extends RealFieldElement<T>> FieldPVCoordinates<T> transformPVCoordinates(final FieldPVCoordinates<T> pv) {
584 
585         // apply translation
586         final FieldVector3D<T> intermediateP = pv.getPosition().add(cartesian.getPosition());
587         final FieldVector3D<T> intermediateV = pv.getVelocity().add(cartesian.getVelocity());
588         final FieldVector3D<T> intermediateA = pv.getAcceleration().add(cartesian.getAcceleration());
589 
590         // apply rotation
591         final FieldVector3D<T> transformedP = FieldRotation.applyTo(angular.getRotation(), intermediateP);
592         final FieldVector3D<T> crossP       = FieldVector3D.crossProduct(angular.getRotationRate(), transformedP);
593         final FieldVector3D<T> transformedV = FieldRotation.applyTo(angular.getRotation(), intermediateV).subtract(crossP);
594         final FieldVector3D<T> crossV       = FieldVector3D.crossProduct(angular.getRotationRate(), transformedV);
595         final FieldVector3D<T> crossCrossP  = FieldVector3D.crossProduct(angular.getRotationRate(), crossP);
596         final FieldVector3D<T> crossDotP    = FieldVector3D.crossProduct(angular.getRotationAcceleration(), transformedP);
597         final FieldVector3D<T> transformedA =
598                 new FieldVector3D<T>( 1, FieldRotation.applyTo(angular.getRotation(), intermediateA),
599                                      -2, crossV,
600                                      -1, crossCrossP,
601                                      -1, crossDotP);
602 
603         // build transformed object
604         return new FieldPVCoordinates<T>(transformedP, transformedV, transformedA);
605 
606     }
607 
608     /** Transform {@link TimeStampedFieldPVCoordinates} including kinematic effects.
609      * <p>
610      * In order to allow the user more flexibility, this method does <em>not</em> check for
611      * consistency between the transform {@link #getDate() date} and the time-stamped
612      * position-velocity {@link TimeStampedFieldPVCoordinates#getDate() date}. The returned
613      * value will always have the same {@link TimeStampedFieldPVCoordinates#getDate() date} as
614      * the input argument, regardless of the instance {@link #getDate() date}.
615      * </p>
616      * @param pv time-stamped position-velocity to transform.
617      * @param <T> type of the field elements
618      * @return transformed time-stamped position-velocity
619      * @since 7.0
620      */
621     public <T extends RealFieldElement<T>> TimeStampedFieldPVCoordinates<T> transformPVCoordinates(final TimeStampedFieldPVCoordinates<T> pv) {
622 
623         // apply translation
624         final FieldVector3D<T> intermediateP = pv.getPosition().add(cartesian.getPosition());
625         final FieldVector3D<T> intermediateV = pv.getVelocity().add(cartesian.getVelocity());
626         final FieldVector3D<T> intermediateA = pv.getAcceleration().add(cartesian.getAcceleration());
627 
628         // apply rotation
629         final FieldVector3D<T> transformedP = FieldRotation.applyTo(angular.getRotation(), intermediateP);
630         final FieldVector3D<T> crossP       = FieldVector3D.crossProduct(angular.getRotationRate(), transformedP);
631         final FieldVector3D<T> transformedV = FieldRotation.applyTo(angular.getRotation(), intermediateV).subtract(crossP);
632         final FieldVector3D<T> crossV       = FieldVector3D.crossProduct(angular.getRotationRate(), transformedV);
633         final FieldVector3D<T> crossCrossP  = FieldVector3D.crossProduct(angular.getRotationRate(), crossP);
634         final FieldVector3D<T> crossDotP    = FieldVector3D.crossProduct(angular.getRotationAcceleration(), transformedP);
635         final FieldVector3D<T> transformedA =
636                 new FieldVector3D<T>( 1, FieldRotation.applyTo(angular.getRotation(), intermediateA),
637                                      -2, crossV,
638                                      -1, crossCrossP,
639                                      -1, crossDotP);
640 
641         // build transformed object
642         return new TimeStampedFieldPVCoordinates<T>(pv.getDate(), transformedP, transformedV, transformedA);
643 
644     }
645 
646     /** Compute the Jacobian of the {@link #transformPVCoordinates(PVCoordinates)}
647      * method of the transform.
648      * <p>
649      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i
650      * of the transformed {@link PVCoordinates} with respect to Cartesian coordinate j
651      * of the input {@link PVCoordinates} in method {@link #transformPVCoordinates(PVCoordinates)}.
652      * </p>
653      * <p>
654      * This definition implies that if we define position-velocity coordinates
655      * <pre>
656      * PV₁ = transform.transformPVCoordinates(PV₀), then
657      * </pre>
658      * their differentials dPV₁ and dPV₀ will obey the following relation
659      * where J is the matrix computed by this method:<br/>
660      * <pre>
661      * dPV₁ = J &times; dPV₀
662      * </pre>
663      * </p>
664      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with
665      * the Jacobian, only the upper left 6x6 corner will be modified
666      * @deprecated as of 7.0, replaced with {@link #getJacobian(CartesianDerivativesFilter, double[][])}
667      */
668     @Deprecated
669     public void getJacobian(final double[][] jacobian) {
670         getJacobian(CartesianDerivativesFilter.USE_PV, jacobian);
671     }
672 
673     /** Compute the Jacobian of the {@link #transformPVCoordinates(PVCoordinates)}
674      * method of the transform.
675      * <p>
676      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i
677      * of the transformed {@link PVCoordinates} with respect to Cartesian coordinate j
678      * of the input {@link PVCoordinates} in method {@link #transformPVCoordinates(PVCoordinates)}.
679      * </p>
680      * <p>
681      * This definition implies that if we define position-velocity coordinates
682      * <pre>
683      * PV₁ = transform.transformPVCoordinates(PV₀), then
684      * </pre>
685      * their differentials dPV₁ and dPV₀ will obey the following relation
686      * where J is the matrix computed by this method:<br/>
687      * <pre>
688      * dPV₁ = J &times; dPV₀
689      * </pre>
690      * </p>
691      * @param selector selector specifying the size of the upper left corner that must be filled
692      * (either 3x3 for positions only, 6x6 for positions and velocities, 9x9 for positions,
693      * velocities and accelerations)
694      * @param jacobian placeholder matrix whose upper-left corner is to be filled with
695      * the Jacobian, the rest of the matrix remaining untouched
696      */
697     public void getJacobian(final CartesianDerivativesFilter selector, final double[][] jacobian) {
698 
699         // elementary matrix for rotation
700         final double[][] mData = angular.getRotation().getMatrix();
701 
702         // dP1/dP0
703         System.arraycopy(mData[0], 0, jacobian[0], 0, 3);
704         System.arraycopy(mData[1], 0, jacobian[1], 0, 3);
705         System.arraycopy(mData[2], 0, jacobian[2], 0, 3);
706 
707         if (selector.getMaxOrder() >= 1) {
708 
709             // dP1/dV0
710             Arrays.fill(jacobian[0], 3, 6, 0.0);
711             Arrays.fill(jacobian[1], 3, 6, 0.0);
712             Arrays.fill(jacobian[2], 3, 6, 0.0);
713 
714             // dV1/dP0
715             final Vector3D o = angular.getRotationRate();
716             final double ox = o.getX();
717             final double oy = o.getY();
718             final double oz = o.getZ();
719             for (int i = 0; i < 3; ++i) {
720                 jacobian[3][i] = -(oy * mData[2][i] - oz * mData[1][i]);
721                 jacobian[4][i] = -(oz * mData[0][i] - ox * mData[2][i]);
722                 jacobian[5][i] = -(ox * mData[1][i] - oy * mData[0][i]);
723             }
724 
725             // dV1/dV0
726             System.arraycopy(mData[0], 0, jacobian[3], 3, 3);
727             System.arraycopy(mData[1], 0, jacobian[4], 3, 3);
728             System.arraycopy(mData[2], 0, jacobian[5], 3, 3);
729 
730             if (selector.getMaxOrder() >= 2) {
731 
732                 // dP1/dA0
733                 Arrays.fill(jacobian[0], 6, 9, 0.0);
734                 Arrays.fill(jacobian[1], 6, 9, 0.0);
735                 Arrays.fill(jacobian[2], 6, 9, 0.0);
736 
737                 // dV1/dA0
738                 Arrays.fill(jacobian[3], 6, 9, 0.0);
739                 Arrays.fill(jacobian[4], 6, 9, 0.0);
740                 Arrays.fill(jacobian[5], 6, 9, 0.0);
741 
742                 // dA1/dP0
743                 final Vector3D oDot = angular.getRotationAcceleration();
744                 final double oDotx  = oDot.getX();
745                 final double oDoty  = oDot.getY();
746                 final double oDotz  = oDot.getZ();
747                 for (int i = 0; i < 3; ++i) {
748                     jacobian[6][i] = -(oDoty * mData[2][i] - oDotz * mData[1][i]) - (oy * jacobian[5][i] - oz * jacobian[4][i]);
749                     jacobian[7][i] = -(oDotz * mData[0][i] - oDotx * mData[2][i]) - (oz * jacobian[3][i] - ox * jacobian[5][i]);
750                     jacobian[8][i] = -(oDotx * mData[1][i] - oDoty * mData[0][i]) - (ox * jacobian[4][i] - oy * jacobian[3][i]);
751                 }
752 
753                 // dA1/dV0
754                 for (int i = 0; i < 3; ++i) {
755                     jacobian[6][i + 3] = -2 * (oy * mData[2][i] - oz * mData[1][i]);
756                     jacobian[7][i + 3] = -2 * (oz * mData[0][i] - ox * mData[2][i]);
757                     jacobian[8][i + 3] = -2 * (ox * mData[1][i] - oy * mData[0][i]);
758                 }
759 
760                 // dA1/dA0
761                 System.arraycopy(mData[0], 0, jacobian[6], 6, 3);
762                 System.arraycopy(mData[1], 0, jacobian[7], 6, 3);
763                 System.arraycopy(mData[2], 0, jacobian[8], 6, 3);
764 
765             }
766 
767         }
768 
769     }
770 
771     /** Get the underlying elementary cartesian part.
772      * <p>A transform can be uniquely represented as an elementary
773      * translation followed by an elementary rotation. This method
774      * returns this unique elementary translation with its derivative.</p>
775      * @return underlying elementary cartesian part
776      * @see #getTranslation()
777      * @see #getVelocity()
778      */
779     public PVCoordinates getCartesian() {
780         return cartesian;
781     }
782 
783     /** Get the underlying elementary translation.
784      * <p>A transform can be uniquely represented as an elementary
785      * translation followed by an elementary rotation. This method
786      * returns this unique elementary translation.</p>
787      * @return underlying elementary translation
788      * @see #getCartesian()
789      * @see #getVelocity()
790      * @see #getAcceleration()
791      */
792     public Vector3D getTranslation() {
793         return cartesian.getPosition();
794     }
795 
796     /** Get the first time derivative of the translation.
797      * @return first time derivative of the translation
798      * @see #getCartesian()
799      * @see #getTranslation()
800      * @see #getAcceleration()
801      */
802     public Vector3D getVelocity() {
803         return cartesian.getVelocity();
804     }
805 
806     /** Get the second time derivative of the translation.
807      * @return second time derivative of the translation
808      * @see #getCartesian()
809      * @see #getTranslation()
810      * @see #getVelocity()
811      */
812     public Vector3D getAcceleration() {
813         return cartesian.getAcceleration();
814     }
815 
816     /** Get the underlying elementary angular part.
817      * <p>A transform can be uniquely represented as an elementary
818      * translation followed by an elementary rotation. This method
819      * returns this unique elementary rotation with its derivative.</p>
820      * @return underlying elementary angular part
821      * @see #getRotation()
822      * @see #getRotationRate()
823      * @see #getRotationAcceleration()
824      */
825     public AngularCoordinates getAngular() {
826         return angular;
827     }
828 
829     /** Get the underlying elementary rotation.
830      * <p>A transform can be uniquely represented as an elementary
831      * translation followed by an elementary rotation. This method
832      * returns this unique elementary rotation.</p>
833      * @return underlying elementary rotation
834      * @see #getAngular()
835      * @see #getRotationRate()
836      * @see #getRotationAcceleration()
837      */
838     public Rotation getRotation() {
839         return angular.getRotation();
840     }
841 
842     /** Get the first time derivative of the rotation.
843      * <p>The norm represents the angular rate.</p>
844      * @return First time derivative of the rotation
845      * @see #getAngular()
846      * @see #getRotation()
847      * @see #getRotationAcceleration()
848      */
849     public Vector3D getRotationRate() {
850         return angular.getRotationRate();
851     }
852 
853     /** Get the second time derivative of the rotation.
854      * @return Second time derivative of the rotation
855      * @see #getAngular()
856      * @see #getRotation()
857      * @see #getRotationRate()
858      */
859     public Vector3D getRotationAcceleration() {
860         return angular.getRotationAcceleration();
861     }
862 
863     /** Specialized class for identity transform. */
864     private static class IdentityTransform extends Transform {
865 
866         /** Serializable UID. */
867         private static final long serialVersionUID = -9042082036141830517L;
868 
869         /** Simple constructor. */
870         public IdentityTransform() {
871             super(AbsoluteDate.J2000_EPOCH, PVCoordinates.ZERO, AngularCoordinates.IDENTITY);
872         }
873 
874         /** {@inheritDoc} */
875         @Override
876         public Transform shiftedBy(final double dt) {
877             return this;
878         }
879 
880         /** {@inheritDoc} */
881         @Override
882         public Transform getInverse() {
883             return this;
884         };
885 
886         /** {@inheritDoc} */
887         @Override
888         public Vector3D transformPosition(final Vector3D position) {
889             return position;
890         }
891 
892         /** {@inheritDoc} */
893         @Override
894         public Vector3D transformVector(final Vector3D vector) {
895             return vector;
896         }
897 
898         /** {@inheritDoc} */
899         @Override
900         public Line transformLine(final Line line) {
901             return line;
902         }
903 
904         /** {@inheritDoc} */
905         @Override
906         public PVCoordinates transformPVCoordinates(final PVCoordinates pv) {
907             return pv;
908         }
909 
910         /** {@inheritDoc} */
911         @Override
912         public void getJacobian(final CartesianDerivativesFilter selector, final double[][] jacobian) {
913             final int n = 3 * (selector.getMaxOrder() + 1);
914             for (int i = 0; i < n; ++i) {
915                 Arrays.fill(jacobian[i], 0, n, 0.0);
916                 jacobian[i][i] = 1.0;
917             }
918         }
919 
920     }
921 
922 }