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 × 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 × 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 }