1   /* Copyright 2002-2026 Airbus Defence and Space
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    * Airbus Defence and Space licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.analytical.intelsat;
18  
19  import java.util.Collections;
20  import java.util.List;
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.FieldSinCos;
27  import org.orekit.annotation.DefaultDataContext;
28  import org.orekit.attitudes.AttitudeProvider;
29  import org.orekit.attitudes.FieldAttitude;
30  import org.orekit.attitudes.FrameAlignedProvider;
31  import org.orekit.data.DataContext;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.errors.OrekitMessages;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.orbits.FieldCartesianOrbit;
37  import org.orekit.orbits.FieldOrbit;
38  import org.orekit.propagation.FieldSpacecraftState;
39  import org.orekit.propagation.Propagator;
40  import org.orekit.propagation.analytical.FieldAbstractAnalyticalPropagator;
41  import org.orekit.time.FieldAbsoluteDate;
42  import org.orekit.utils.Constants;
43  import org.orekit.utils.FieldPVCoordinates;
44  import org.orekit.utils.IERSConventions;
45  import org.orekit.utils.ParameterDriver;
46  import org.orekit.utils.units.Unit;
47  
48  /**
49   * This class provides elements to propagate Intelsat's 11 elements.
50   * <p>
51   * Intelsat's 11 elements propagation is defined in ITU-R S.1525 standard.
52   * </p>
53   *
54   * @param <T> type of the field elements
55   * @author Bryan Cazabonne
56   * @since 12.1
57   */
58  public class FieldIntelsatElevenElementsPropagator<T extends CalculusFieldElement<T>> extends FieldAbstractAnalyticalPropagator<T> {
59  
60      /**
61       * Intelsat's 11 elements.
62       */
63      private final FieldIntelsatElevenElements<T> elements;
64  
65      /**
66       * Inertial frame for the output orbit.
67       */
68      private final Frame inertialFrame;
69  
70      /**
71       * ECEF frame related to the Intelsat's 11 elements.
72       */
73      private final Frame ecefFrame;
74  
75      /**
76       * Spacecraft mass in kilograms.
77       */
78      private final T mass;
79  
80      /**
81       * Compute spacecraft's east longitude.
82       */
83      private FieldUnivariateDerivative2<T> eastLongitudeDegrees;
84  
85      /**
86       * Compute spacecraft's geocentric latitude.
87       */
88      private FieldUnivariateDerivative2<T> geocentricLatitudeDegrees;
89  
90      /**
91       * Compute spacecraft's orbit radius.
92       */
93      private FieldUnivariateDerivative2<T> orbitRadius;
94  
95      /**
96       * Default constructor.
97       * <p>
98       * This constructor uses the {@link DataContext#getDefault() default data context}.
99       * </p>
100      * <p> The attitude provider is set by default to be aligned with the inertial frame.<br>
101      * The mass is set by default to the
102      * {@link org.orekit.propagation.Propagator#DEFAULT_MASS DEFAULT_MASS}.<br>
103      * The inertial frame is set by default to the
104      * {@link org.orekit.frames.Predefined#TOD_CONVENTIONS_2010_SIMPLE_EOP TOD frame} in the default data
105      * context.<br>
106      * The ECEF frame is set by default to the
107      * {@link org.orekit.frames.Predefined#ITRF_CIO_CONV_2010_SIMPLE_EOP
108      * CIO/2010-based ITRF simple EOP} in the default data context.
109      * </p>
110      *
111      * @param elements Intelsat's 11 elements
112      */
113     @DefaultDataContext
114     public FieldIntelsatElevenElementsPropagator(final FieldIntelsatElevenElements<T> elements) {
115         this(elements, FramesFactory.getTOD(IERSConventions.IERS_2010, true), FramesFactory.getITRF(IERSConventions.IERS_2010, true));
116     }
117 
118     /**
119      * Constructor.
120      *
121      * <p> The attitude provider is set by default to be aligned with the inertial frame.<br>
122      * The mass is set by default to the
123      * {@link org.orekit.propagation.Propagator#DEFAULT_MASS DEFAULT_MASS}.<br>
124      * </p>
125      *
126      * @param elements      Intelsat's 11 elements
127      * @param inertialFrame inertial frame for the output orbit
128      * @param ecefFrame     ECEF frame related to the Intelsat's 11 elements
129      */
130     public FieldIntelsatElevenElementsPropagator(final FieldIntelsatElevenElements<T> elements, final Frame inertialFrame, final Frame ecefFrame) {
131         this(elements, inertialFrame, ecefFrame, FrameAlignedProvider.of(inertialFrame), elements.getEpoch().getField().getZero().add(Propagator.DEFAULT_MASS));
132     }
133 
134     /**
135      * Constructor.
136      *
137      * @param elements         Intelsat's 11 elements
138      * @param inertialFrame    inertial frame for the output orbit
139      * @param ecefFrame        ECEF frame related to the Intelsat's 11 elements
140      * @param attitudeProvider attitude provider
141      * @param mass             spacecraft mass
142      */
143     public FieldIntelsatElevenElementsPropagator(final FieldIntelsatElevenElements<T> elements, final Frame inertialFrame, final Frame ecefFrame,
144                                                  final AttitudeProvider attitudeProvider, final T mass) {
145         super(elements.getEpoch().getField(), attitudeProvider);
146         this.elements = elements;
147         this.inertialFrame = inertialFrame;
148         this.ecefFrame = ecefFrame;
149         this.mass = mass;
150         setStartDate(elements.getEpoch());
151         final FieldOrbit<T> orbitAtElementsDate = propagateOrbit(elements.getEpoch(), getParameters(elements.getEpoch().getField()));
152         final FieldAttitude<T> attitude = attitudeProvider.getAttitude(orbitAtElementsDate, elements.getEpoch(), inertialFrame);
153         super.resetInitialState(new FieldSpacecraftState<>(orbitAtElementsDate, attitude).withMass(mass));
154     }
155 
156     /**
157      * Converts the Intelsat's 11 elements into Position/Velocity coordinates in ECEF.
158      *
159      * @param date computation epoch
160      * @return Position/Velocity coordinates in ECEF
161      */
162     public FieldPVCoordinates<T> propagateInEcef(final FieldAbsoluteDate<T> date) {
163         final Field<T> field = date.getField();
164         final FieldUnivariateDerivative2<T> tDays = new FieldUnivariateDerivative2<>(date.durationFrom(elements.getEpoch()), field.getOne(), field.getZero()).divide(
165                 Constants.JULIAN_DAY);
166         final T wDegreesPerDay = elements.getLm1().add(IntelsatElevenElements.DRIFT_RATE_SHIFT_DEG_PER_DAY);
167         final FieldUnivariateDerivative2<T> wt = FastMath.toRadians(tDays.multiply(wDegreesPerDay));
168         final FieldSinCos<FieldUnivariateDerivative2<T>> scWt = FastMath.sinCos(wt);
169         final FieldSinCos<FieldUnivariateDerivative2<T>> sc2Wt = FastMath.sinCos(wt.multiply(2.0));
170         final FieldUnivariateDerivative2<T> satelliteEastLongitudeDegrees = computeSatelliteEastLongitudeDegrees(tDays, scWt, sc2Wt);
171         final FieldUnivariateDerivative2<T> satelliteGeocentricLatitudeDegrees = computeSatelliteGeocentricLatitudeDegrees(tDays, scWt);
172         final FieldUnivariateDerivative2<T> satelliteRadius = computeSatelliteRadiusKilometers(wDegreesPerDay, scWt).multiply(Unit.KILOMETRE.getScale());
173         this.eastLongitudeDegrees = satelliteEastLongitudeDegrees;
174         this.geocentricLatitudeDegrees = satelliteGeocentricLatitudeDegrees;
175         this.orbitRadius = satelliteRadius;
176         final FieldSinCos<FieldUnivariateDerivative2<T>> scLongitude = FastMath.sinCos(FastMath.toRadians(satelliteEastLongitudeDegrees));
177         final FieldSinCos<FieldUnivariateDerivative2<T>> scLatitude = FastMath.sinCos(FastMath.toRadians(satelliteGeocentricLatitudeDegrees));
178         final FieldVector3D<FieldUnivariateDerivative2<T>> positionWithDerivatives = new FieldVector3D<>(satelliteRadius.multiply(scLatitude.cos()).multiply(scLongitude.cos()),
179                                                                                                          satelliteRadius.multiply(scLatitude.cos()).multiply(scLongitude.sin()),
180                                                                                                          satelliteRadius.multiply(scLatitude.sin()));
181         return new FieldPVCoordinates<>(new FieldVector3D<>(positionWithDerivatives.getX().getValue(), //
182                                                             positionWithDerivatives.getY().getValue(), //
183                                                             positionWithDerivatives.getZ().getValue()), //
184                                         new FieldVector3D<>(positionWithDerivatives.getX().getFirstDerivative(), //
185                                                             positionWithDerivatives.getY().getFirstDerivative(), //
186                                                             positionWithDerivatives.getZ().getFirstDerivative()), //
187                                         new FieldVector3D<>(positionWithDerivatives.getX().getSecondDerivative(), //
188                                                             positionWithDerivatives.getY().getSecondDerivative(), //
189                                                             positionWithDerivatives.getZ().getSecondDerivative()));
190     }
191 
192     /**
193      * {@inheritDoc}.
194      */
195     @Override
196     public void resetInitialState(final FieldSpacecraftState<T> state) {
197         throw new OrekitException(OrekitMessages.NON_RESETABLE_STATE);
198     }
199 
200     /**
201      * {@inheritDoc}.
202      */
203     @Override
204     protected T getMass(final FieldAbsoluteDate<T> date) {
205         return mass;
206     }
207 
208     /**
209      * {@inheritDoc}.
210      */
211     @Override
212     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
213         throw new OrekitException(OrekitMessages.NON_RESETABLE_STATE);
214     }
215 
216     /**
217      * {@inheritDoc}.
218      */
219     @Override
220     public FieldOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date,
221                                         final T[] parameters) {
222         return new FieldCartesianOrbit<>(ecefFrame.getTransformTo(inertialFrame, date).transformPVCoordinates(propagateInEcef(date)), inertialFrame, date,
223                                          date.getField().getZero().add(Constants.WGS84_EARTH_MU));
224     }
225 
226     /**
227      * Computes the satellite's east longitude.
228      *
229      * @param tDays delta time in days
230      * @param scW   sin/cos of the W angle
231      * @param sc2W  sin/cos of the 2xW angle
232      * @return the satellite's east longitude in degrees
233      */
234     private FieldUnivariateDerivative2<T> computeSatelliteEastLongitudeDegrees(final FieldUnivariateDerivative2<T> tDays, final FieldSinCos<FieldUnivariateDerivative2<T>> scW,
235                                                                                final FieldSinCos<FieldUnivariateDerivative2<T>> sc2W) {
236         final FieldUnivariateDerivative2<T> longitude = tDays.multiply(tDays).multiply(elements.getLm2()) //
237                                                              .add(tDays.multiply(elements.getLm1())) //
238                                                              .add(elements.getLm0());
239         final FieldUnivariateDerivative2<T> cosineLongitudeTerm = scW.cos().multiply(tDays.multiply(elements.getLonC1()).add(elements.getLonC()));
240         final FieldUnivariateDerivative2<T> sineLongitudeTerm = scW.sin().multiply(tDays.multiply(elements.getLonS1()).add(elements.getLonS()));
241         final FieldUnivariateDerivative2<T> latitudeTerm = sc2W.sin()
242                                                                .multiply(elements.getLatC()
243                                                                                  .multiply(elements.getLatC())
244                                                                                  .subtract(elements.getLatS().multiply(elements.getLatS()))
245                                                                                  .multiply(0.5)) //
246                                                                .subtract(sc2W.cos().multiply(elements.getLatC().multiply(elements.getLatS()))) //
247                                                                .multiply(IntelsatElevenElements.K);
248         return longitude.add(cosineLongitudeTerm).add(sineLongitudeTerm).add(latitudeTerm);
249     }
250 
251     /**
252      * Computes the satellite's geocentric latitude.
253      *
254      * @param tDays delta time in days
255      * @param scW   sin/cos of the W angle
256      * @return he satellite geocentric latitude in degrees
257      */
258     private FieldUnivariateDerivative2<T> computeSatelliteGeocentricLatitudeDegrees(final FieldUnivariateDerivative2<T> tDays,
259                                                                                     final FieldSinCos<FieldUnivariateDerivative2<T>> scW) {
260         final FieldUnivariateDerivative2<T> cosineTerm = scW.cos().multiply(tDays.multiply(elements.getLatC1()).add(elements.getLatC()));
261         final FieldUnivariateDerivative2<T> sineTerm = scW.sin().multiply(tDays.multiply(elements.getLatS1()).add(elements.getLatS()));
262         return cosineTerm.add(sineTerm);
263     }
264 
265     /**
266      * Computes the satellite's orbit radius.
267      *
268      * @param wDegreesPerDay W angle in degrees/day
269      * @param scW            sin/cos of the W angle
270      * @return the satellite's orbit radius in kilometers
271      */
272     private FieldUnivariateDerivative2<T> computeSatelliteRadiusKilometers(final T wDegreesPerDay, final FieldSinCos<FieldUnivariateDerivative2<T>> scW) {
273         final T coefficient = elements.getLm1()
274                                       .multiply(2.0)
275                                       .divide(wDegreesPerDay.subtract(elements.getLm1()).multiply(3.0))
276                                       .negate()
277                                       .add(1.0)
278                                       .multiply(IntelsatElevenElements.SYNCHRONOUS_RADIUS_KM);
279         return scW.sin()
280                   .multiply(elements.getLonC().multiply(IntelsatElevenElements.K))
281                   .add(1.0)
282                   .subtract(scW.cos().multiply(elements.getLonS().multiply(IntelsatElevenElements.K)))
283                   .multiply(coefficient);
284     }
285 
286     /**
287      * Get the computed satellite's east longitude.
288      *
289      * @return the satellite's east longitude in degrees
290      */
291     public FieldUnivariateDerivative2<T> getEastLongitudeDegrees() {
292         return eastLongitudeDegrees;
293     }
294 
295     /**
296      * Get the computed satellite's geocentric latitude.
297      *
298      * @return the satellite's geocentric latitude in degrees
299      */
300     public FieldUnivariateDerivative2<T> getGeocentricLatitudeDegrees() {
301         return geocentricLatitudeDegrees;
302     }
303 
304     /**
305      * Get the computed satellite's orbit.
306      *
307      * @return satellite's orbit radius in meters
308      */
309     public FieldUnivariateDerivative2<T> getOrbitRadius() {
310         return orbitRadius;
311     }
312 
313     /**
314      * {@inheritDoc}.
315      */
316     @Override
317     public Frame getFrame() {
318         return inertialFrame;
319     }
320 
321     /**
322      * {@inheritDoc}.
323      */
324     @Override
325     public List<ParameterDriver> getParametersDrivers() {
326         return Collections.emptyList();
327     }
328 
329     /**
330      * Get the Intelsat's 11 elements used by the propagator.
331      *
332      * @return the Intelsat's 11 elements used by the propagator
333      */
334     public FieldIntelsatElevenElements<T> getIntelsatElevenElements() {
335         return elements;
336     }
337 }