1 /* Copyright 2002-2021 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
19 import org.hipparchus.Field;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Rotation;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.orekit.frames.FieldTransform;
26 import org.orekit.frames.Frame;
27 import org.orekit.frames.Transform;
28 import org.orekit.time.AbsoluteDate;
29 import org.orekit.time.FieldAbsoluteDate;
30 import org.orekit.utils.FieldPVCoordinates;
31 import org.orekit.utils.FieldPVCoordinatesProvider;
32 import org.orekit.utils.PVCoordinates;
33 import org.orekit.utils.PVCoordinatesProvider;
34
35
36 /**
37 * This class handles a celestial body pointed attitude provider.
38 * <p>The celestial body pointed law is defined by two main elements:
39 * <ul>
40 * <li>a celestial body towards which some satellite axis is exactly aimed</li>
41 * <li>a phasing reference defining the rotation around the pointing axis</li>
42 * </ul>
43 *
44 * <p>
45 * The celestial body implicitly defines two of the three degrees of freedom
46 * and the phasing reference defines the remaining degree of freedom. This definition
47 * can be represented as first aligning exactly the satellite pointing axis to
48 * the current direction of the celestial body, and then to find the rotation
49 * around this axis such that the satellite phasing axis is in the half-plane
50 * defined by a cut line on the pointing axis and containing the celestial
51 * phasing reference.
52 * </p>
53 * <p>
54 * In order for this definition to work, the user must ensure that the phasing
55 * reference is <strong>never</strong> aligned with the pointing reference.
56 * Since the pointed body moves as the date changes, this should be ensured
57 * regardless of the date. A simple way to do this for Sun, Moon or any planet
58 * pointing is to choose a phasing reference far from the ecliptic plane. Using
59 * <code>Vector3D.PLUS_K</code>, the equatorial pole, is perfect in these cases.
60 * </p>
61 * <p>Instances of this class are guaranteed to be immutable.</p>
62 * @author Luc Maisonobe
63 */
64 public class CelestialBodyPointed implements AttitudeProvider {
65
66 /** Frame in which {@link #phasingCel} is defined. */
67 private final Frame celestialFrame;
68
69 /** Celestial body to point at. */
70 private final PVCoordinatesProvider pointedBody;
71
72 /** Phasing reference, in celestial frame. */
73 private final Vector3D phasingCel;
74
75 /** Satellite axis aiming at the pointed body, in satellite frame. */
76 private final Vector3D pointingSat;
77
78 /** Phasing reference, in satellite frame. */
79 private final Vector3D phasingSat;
80
81 /** Creates new instance.
82 * @param celestialFrame frame in which <code>phasingCel</code> is defined
83 * @param pointedBody celestial body to point at
84 * @param phasingCel phasing reference, in celestial frame
85 * @param pointingSat satellite vector defining the pointing direction
86 * @param phasingSat phasing reference, in satellite frame
87 */
88 public CelestialBodyPointed(final Frame celestialFrame,
89 final PVCoordinatesProvider pointedBody,
90 final Vector3D phasingCel,
91 final Vector3D pointingSat,
92 final Vector3D phasingSat) {
93 this.celestialFrame = celestialFrame;
94 this.pointedBody = pointedBody;
95 this.phasingCel = phasingCel;
96 this.pointingSat = pointingSat;
97 this.phasingSat = phasingSat;
98 }
99
100 /** {@inheritDoc} */
101 public Attitude getAttitude(final PVCoordinatesProvider pvProv,
102 final AbsoluteDate date, final Frame frame) {
103
104 final PVCoordinates satPV = pvProv.getPVCoordinates(date, celestialFrame);
105
106 // compute celestial references at the specified date
107 final PVCoordinates bodyPV = pointedBody.getPVCoordinates(date, celestialFrame);
108 final PVCoordinates pointing = new PVCoordinates(satPV, bodyPV);
109 final Vector3D pointingP = pointing.getPosition();
110 final double r2 = Vector3D.dotProduct(pointingP, pointingP);
111
112 // evaluate instant rotation axis due to sat and body motion only (no phasing yet)
113 final Vector3D rotAxisCel =
114 new Vector3D(1 / r2, Vector3D.crossProduct(pointingP, pointing.getVelocity()));
115
116 // fix instant rotation to take phasing constraint into account
117 // (adding a rotation around pointing axis ensuring the motion of the phasing axis
118 // is constrained in the pointing-phasing plane)
119 final Vector3D v1 = Vector3D.crossProduct(rotAxisCel, phasingCel);
120 final Vector3D v2 = Vector3D.crossProduct(pointingP, phasingCel);
121 final double compensation = -Vector3D.dotProduct(v1, v2) / v2.getNormSq();
122 final Vector3D phasedRotAxisCel = new Vector3D(1.0, rotAxisCel, compensation, pointingP);
123
124 // compute transform from celestial frame to satellite frame
125 final Rotation celToSatRotation =
126 new Rotation(pointingP, phasingCel, pointingSat, phasingSat);
127
128 // build transform combining rotation and instant rotation axis
129 Transform transform = new Transform(date, celToSatRotation, celToSatRotation.applyTo(phasedRotAxisCel));
130 if (frame != celestialFrame) {
131 // prepend transform from specified frame to celestial frame
132 transform = new Transform(date, frame.getTransformTo(celestialFrame, date), transform);
133 }
134
135 // build the attitude
136 return new Attitude(date, frame, transform.getRotation(), transform.getRotationRate(), transform.getRotationAcceleration());
137
138 }
139
140 /** {@inheritDoc} */
141 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(final FieldPVCoordinatesProvider<T> pvProv,
142 final FieldAbsoluteDate<T> date,
143 final Frame frame) {
144
145 final Field<T> field = date.getField();
146 final FieldPVCoordinates<T> satPV = pvProv.getPVCoordinates(date, celestialFrame);
147
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
156 // evaluate instant rotation axis due to sat and body motion only (no phasing yet)
157 final FieldVector3D<T> rotAxisCel =
158 new FieldVector3D<>(r2.reciprocal(), FieldVector3D.crossProduct(pointingP, pointing.getVelocity()));
159
160 // fix instant rotation to take phasing constraint into account
161 // (adding a rotation around pointing axis ensuring the motion of the phasing axis
162 // is constrained in the pointing-phasing plane)
163 final FieldVector3D<T> v1 = FieldVector3D.crossProduct(rotAxisCel, phasingCel);
164 final FieldVector3D<T> v2 = FieldVector3D.crossProduct(pointingP, phasingCel);
165 final T compensation = FieldVector3D.dotProduct(v1, v2).negate().divide(v2.getNormSq());
166 final FieldVector3D<T> phasedRotAxisCel = new FieldVector3D<>(field.getOne(), rotAxisCel, compensation, pointingP);
167
168 // compute transform from celestial frame to satellite frame
169 final FieldRotation<T> celToSatRotation =
170 new FieldRotation<>(pointingP, new FieldVector3D<>(field, phasingCel),
171 new FieldVector3D<>(field, pointingSat), new FieldVector3D<>(field, phasingSat));
172
173 // build transform combining rotation and instant rotation axis
174 FieldTransform<T> transform = new FieldTransform<>(date, celToSatRotation, celToSatRotation.applyTo(phasedRotAxisCel));
175 if (frame != celestialFrame) {
176 // prepend transform from specified frame to celestial frame
177 transform = new FieldTransform<>(date, frame.getTransformTo(celestialFrame, date), transform);
178 }
179
180 // build the attitude
181 return new FieldAttitude<>(date, frame,
182 transform.getRotation(), transform.getRotationRate(), transform.getRotationAcceleration());
183
184 }
185
186 }