CelestialBodyPointed.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (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.attitudes;

  18. import org.hipparchus.Field;
  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.geometry.euclidean.threed.Rotation;
  23. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  24. import org.orekit.frames.FieldStaticTransform;
  25. import org.orekit.frames.FieldTransform;
  26. import org.orekit.frames.Frame;
  27. import org.orekit.frames.StaticTransform;
  28. import org.orekit.frames.Transform;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.time.FieldAbsoluteDate;
  31. import org.orekit.utils.FieldPVCoordinates;
  32. import org.orekit.utils.FieldPVCoordinatesProvider;
  33. import org.orekit.utils.PVCoordinates;
  34. import org.orekit.utils.PVCoordinatesProvider;


  35. /**
  36.  * This class handles a celestial body pointed attitude provider.
  37.  * <p>The celestial body pointed law is defined by two main elements:
  38.  * <ul>
  39.  *   <li>a celestial body towards which some satellite axis is exactly aimed</li>
  40.  *   <li>a phasing reference defining the rotation around the pointing axis</li>
  41.  * </ul>
  42.  *
  43.  * <p>
  44.  * The celestial body implicitly defines two of the three degrees of freedom
  45.  * and the phasing reference defines the remaining degree of freedom. This definition
  46.  * can be represented as first aligning exactly the satellite pointing axis to
  47.  * the current direction of the celestial body, and then to find the rotation
  48.  * around this axis such that the satellite phasing axis is in the half-plane
  49.  * defined by a cut line on the pointing axis and containing the celestial
  50.  * phasing reference.
  51.  * </p>
  52.  * <p>
  53.  * In order for this definition to work, the user must ensure that the phasing
  54.  * reference is <strong>never</strong> aligned with the pointing reference.
  55.  * Since the pointed body moves as the date changes, this should be ensured
  56.  * regardless of the date. A simple way to do this for Sun, Moon or any planet
  57.  * pointing is to choose a phasing reference far from the ecliptic plane. Using
  58.  * <code>Vector3D.PLUS_K</code>, the equatorial pole, is perfect in these cases.
  59.  * </p>
  60.  * <p>Instances of this class are guaranteed to be immutable.</p>
  61.  * @author Luc Maisonobe
  62.  */
  63. public class CelestialBodyPointed implements AttitudeProvider {

  64.     /** Frame in which {@link #phasingCel} is defined. */
  65.     private final Frame celestialFrame;

  66.     /** Celestial body to point at. */
  67.     private final PVCoordinatesProvider pointedBody;

  68.     /** Phasing reference, in celestial frame. */
  69.     private final Vector3D phasingCel;

  70.     /** Satellite axis aiming at the pointed body, in satellite frame. */
  71.     private final Vector3D pointingSat;

  72.     /** Phasing reference, in satellite frame. */
  73.     private final Vector3D phasingSat;

  74.     /** Creates new instance.
  75.      * @param celestialFrame frame in which <code>phasingCel</code> is defined
  76.      * @param pointedBody celestial body to point at
  77.      * @param phasingCel phasing reference, in celestial frame
  78.      * @param pointingSat satellite vector defining the pointing direction
  79.      * @param phasingSat phasing reference, in satellite frame
  80.      */
  81.     public CelestialBodyPointed(final Frame celestialFrame,
  82.                                 final PVCoordinatesProvider pointedBody,
  83.                                 final Vector3D phasingCel,
  84.                                 final Vector3D pointingSat,
  85.                                 final Vector3D phasingSat) {
  86.         this.celestialFrame = celestialFrame;
  87.         this.pointedBody    = pointedBody;
  88.         this.phasingCel     = phasingCel;
  89.         this.pointingSat    = pointingSat;
  90.         this.phasingSat     = phasingSat;
  91.     }

  92.     /** {@inheritDoc} */
  93.     @Override
  94.     public Attitude getAttitude(final PVCoordinatesProvider pvProv,
  95.                                 final AbsoluteDate date, final Frame frame) {

  96.         final PVCoordinates satPV = pvProv.getPVCoordinates(date, celestialFrame);

  97.         // compute celestial references at the specified date
  98.         final PVCoordinates bodyPV    = pointedBody.getPVCoordinates(date, celestialFrame);
  99.         final PVCoordinates pointing  = new PVCoordinates(satPV, bodyPV);
  100.         final Vector3D      pointingP = pointing.getPosition();
  101.         final double        r2        = Vector3D.dotProduct(pointingP, pointingP);

  102.         // evaluate instant rotation axis due to sat and body motion only (no phasing yet)
  103.         final Vector3D rotAxisCel = new Vector3D(1 / r2, Vector3D.crossProduct(pointingP, pointing.getVelocity()));

  104.         // fix instant rotation to take phasing constraint into account
  105.         // (adding a rotation around pointing axis ensuring the motion of the phasing axis
  106.         //  is constrained in the pointing-phasing plane)
  107.         final Vector3D v1    = Vector3D.crossProduct(rotAxisCel, phasingCel);
  108.         final Vector3D v2    = Vector3D.crossProduct(pointingP,  phasingCel);
  109.         final double   compensation = -Vector3D.dotProduct(v1, v2) / v2.getNormSq();
  110.         final Vector3D phasedRotAxisCel = new Vector3D(1.0, rotAxisCel, compensation, pointingP);

  111.         // compute transform from celestial frame to satellite frame
  112.         final Rotation celToSatRotation =
  113.             new Rotation(pointingP, phasingCel, pointingSat, phasingSat);

  114.         // build transform combining rotation and instant rotation axis
  115.         Transform transform = new Transform(date, celToSatRotation, celToSatRotation.applyTo(phasedRotAxisCel));
  116.         if (frame != celestialFrame) {
  117.             // prepend transform from specified frame to celestial frame
  118.             transform = new Transform(date, frame.getTransformTo(celestialFrame, date), transform);
  119.         }

  120.         // build the attitude
  121.         return new Attitude(date, frame, transform.getRotation(), transform.getRotationRate(), transform.getRotationAcceleration());

  122.     }

  123.     /** {@inheritDoc} */
  124.     @Override
  125.     public Rotation getAttitudeRotation(final PVCoordinatesProvider pvProv,
  126.                                         final AbsoluteDate date,
  127.                                         final Frame frame) {
  128.         final Vector3D satPosition = pvProv.getPosition(date, celestialFrame);

  129.         // compute celestial references at the specified date
  130.         final Vector3D bodyPosition    = pointedBody.getPosition(date, celestialFrame);
  131.         final Vector3D      pointingP  = bodyPosition.subtract(satPosition);

  132.         // compute static transform from celestial frame to satellite frame
  133.         final Rotation celToSatRotation = new Rotation(pointingP, phasingCel, pointingSat, phasingSat);
  134.         StaticTransform staticTransform = StaticTransform.of(date, celToSatRotation);

  135.         if (frame != celestialFrame) {
  136.             // prepend static transform from specified frame to celestial frame
  137.             staticTransform = StaticTransform.compose(date, frame.getStaticTransformTo(celestialFrame, date), staticTransform);
  138.         }
  139.         return staticTransform.getRotation();
  140.     }

  141.     /** {@inheritDoc} */
  142.     @Override
  143.     public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(final FieldPVCoordinatesProvider<T> pvProv,
  144.                                                                             final FieldAbsoluteDate<T> date,
  145.                                                                             final Frame frame) {

  146.         final Field<T> field = date.getField();
  147.         final FieldPVCoordinates<T> satPV = pvProv.getPVCoordinates(date, celestialFrame);

  148.         // compute celestial references at the specified date
  149.         final FieldPVCoordinates<T> bodyPV    = new FieldPVCoordinates<>(field,
  150.                                                                          pointedBody.getPVCoordinates(date.toAbsoluteDate(),
  151.                                                                                                       celestialFrame));
  152.         final FieldPVCoordinates<T> pointing  = new FieldPVCoordinates<>(satPV, bodyPV);
  153.         final FieldVector3D<T>      pointingP = pointing.getPosition();
  154.         final T                     r2        = FieldVector3D.dotProduct(pointingP, pointingP);

  155.         // evaluate instant rotation axis due to sat and body motion only (no phasing yet)
  156.         final FieldVector3D<T> rotAxisCel =
  157.             new FieldVector3D<>(r2.reciprocal(), FieldVector3D.crossProduct(pointingP, pointing.getVelocity()));

  158.         // fix instant rotation to take phasing constraint into account
  159.         // (adding a rotation around pointing axis ensuring the motion of the phasing axis
  160.         //  is constrained in the pointing-phasing plane)
  161.         final FieldVector3D<T> v1           = FieldVector3D.crossProduct(rotAxisCel, phasingCel);
  162.         final FieldVector3D<T> v2           = FieldVector3D.crossProduct(pointingP,  phasingCel);
  163.         final T                compensation = FieldVector3D.dotProduct(v1, v2).negate().divide(v2.getNormSq());
  164.         final FieldVector3D<T> phasedRotAxisCel = new FieldVector3D<>(field.getOne(), rotAxisCel, compensation, pointingP);

  165.         // compute transform from celestial frame to satellite frame
  166.         final FieldRotation<T> celToSatRotation =
  167.             new FieldRotation<>(pointingP, new FieldVector3D<>(field, phasingCel),
  168.                             new FieldVector3D<>(field, pointingSat), new FieldVector3D<>(field, phasingSat));

  169.         // build transform combining rotation and instant rotation axis
  170.         FieldTransform<T> transform = new FieldTransform<>(date, celToSatRotation, celToSatRotation.applyTo(phasedRotAxisCel));
  171.         if (frame != celestialFrame) {
  172.             // prepend transform from specified frame to celestial frame
  173.             transform = new FieldTransform<>(date, frame.getTransformTo(celestialFrame, date), transform);
  174.         }

  175.         // build the attitude
  176.         return new FieldAttitude<>(date, frame,
  177.                         transform.getRotation(), transform.getRotationRate(), transform.getRotationAcceleration());

  178.     }

  179.     /** {@inheritDoc} */
  180.     @Override
  181.     public <T extends CalculusFieldElement<T>> FieldRotation<T> getAttitudeRotation(final FieldPVCoordinatesProvider<T> pvProv,
  182.                                                                                     final FieldAbsoluteDate<T> date,
  183.                                                                                     final Frame frame) {
  184.         final Field<T> field = date.getField();
  185.         final FieldVector3D<T> satPosition = pvProv.getPosition(date, celestialFrame);

  186.         // compute celestial references at the specified date
  187.         final FieldVector3D<T> bodyPosition    = new FieldVector3D<>(field,
  188.                 pointedBody.getPosition(date.toAbsoluteDate(), celestialFrame));
  189.         final FieldVector3D<T>      pointingP = bodyPosition.subtract(satPosition);

  190.         // compute rotation from celestial frame to satellite frame
  191.         final FieldRotation<T> celToSatRotation =
  192.                 new FieldRotation<>(pointingP, new FieldVector3D<>(field, phasingCel),
  193.                         new FieldVector3D<>(field, pointingSat), new FieldVector3D<>(field, phasingSat));

  194.         // build static transform combining rotation and instant rotation axis
  195.         FieldStaticTransform<T> staticTransform = FieldStaticTransform.of(date, celToSatRotation);
  196.         if (frame != celestialFrame) {
  197.             // prepend static transform from specified frame to celestial frame
  198.             staticTransform = FieldStaticTransform.compose(date, frame.getStaticTransformTo(celestialFrame, date), staticTransform);
  199.         }
  200.         return staticTransform.getRotation();
  201.     }
  202. }