1   /* Copyright 2002-2026 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.estimation.measurements;
18  
19  import java.util.Map;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.differentiation.Gradient;
24  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
25  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
27  import org.hipparchus.geometry.euclidean.threed.Vector3D;
28  import org.hipparchus.util.FastMath;
29  import org.orekit.bodies.BodyShape;
30  import org.orekit.bodies.FieldGeodeticPoint;
31  import org.orekit.bodies.GeodeticPoint;
32  import org.orekit.frames.FieldStaticTransform;
33  import org.orekit.frames.FieldTransform;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.KinematicTransform;
36  import org.orekit.frames.StaticTransform;
37  import org.orekit.frames.TopocentricFrame;
38  import org.orekit.frames.TopocentricTransformProvider;
39  import org.orekit.frames.Transform;
40  import org.orekit.models.earth.displacement.StationDisplacement;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.time.clocks.QuadraticClockModel;
44  import org.orekit.utils.AngularCoordinates;
45  import org.orekit.utils.FieldAngularCoordinates;
46  import org.orekit.utils.FieldPVCoordinates;
47  import org.orekit.utils.FieldPVCoordinatesProvider;
48  import org.orekit.utils.PVCoordinatesProvider;
49  import org.orekit.utils.ParameterDriver;
50  import org.orekit.utils.TimeStampedFieldPVCoordinates;
51  
52  /** Class modeling a ground station that can perform some measurements.
53   * <p>
54   * This class adds a position offset parameter to a base {@link TopocentricFrame
55   * topocentric frame}.
56   * </p>
57   * <p>
58   * Since 9.3, this class also adds a station clock offset parameter, which manages
59   * the value that must be subtracted from the observed measurement date to get the real
60   * physical date at which the measurement was performed (i.e. the offset is negative
61   * if the ground station clock is slow and positive if it is fast).
62   * </p>
63   * <ol>
64   *   <li>station clock offset, controlled by {@link #getClockBiasDriver()}</li>
65   *   <li>station position offset, controlled by {@link #getEastOffsetDriver()},
66   *   {@link #getNorthOffsetDriver()} and {@link #getZenithOffsetDriver()}</li>
67   * </ol>
68   * @author Luc Maisonobe
69   * @since 8.0
70   */
71  public class GroundStation extends AbstractParticipant implements Observer {
72  
73      /** Position offsets scaling factor.
74       * <p>
75       * We use a power of 2 (in fact really 1.0 here) to avoid numeric noise introduction
76       * in the multiplications/divisions sequences.
77       * </p>
78       */
79      private static final double POSITION_OFFSET_SCALE = FastMath.scalb(1.0, 0);
80  
81      /** Base frame associated with the station. */
82      private final TopocentricFrame baseFrame;
83  
84      /** Driver for position offset along the East axis. */
85      private final ParameterDriver eastOffsetDriver;
86  
87      /** Driver for position offset along the North axis. */
88      private final ParameterDriver northOffsetDriver;
89  
90      /** Driver for position offset along the zenith axis. */
91      private final ParameterDriver zenithOffsetDriver;
92  
93      /**
94       * Build a ground station ignoring {@link StationDisplacement station displacements}.
95       * <p> The initial values for the station offset model
96       * ({@link #getClockBiasDriver()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
97       * {@link #getZenithOffsetDriver()}) are set to 0. This implies that as long as these values are not changed, the
98       * offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as some of these models are changed,
99       * the offset frame moves away from the {@link #getBaseFrame() base frame}.
100      * </p>
101      *
102      * @param baseFrame base frame associated with the station, without *any* parametric model (no station offset)
103      * @see #GroundStation(TopocentricFrame, QuadraticClockModel)
104      * @since 13.0
105      */
106     public GroundStation(final TopocentricFrame baseFrame) {
107         this(baseFrame, createEmptyQuadraticClock(baseFrame.getName()));
108     }
109 
110      /**
111      * Simple constructor.
112      * <p>
113      * The initial values for the station offset model
114      * ({@link #getClockBiasDriver()}, {@link #getEastOffsetDriver()}, {@link #getNorthOffsetDriver()},
115      * {@link #getZenithOffsetDriver()}, {@link #getClockBiasDriver()}) are set to 0. This implies that as long as
116      * these values are not changed, the offset frame is the same as the {@link #getBaseFrame() base frame}. As soon as
117      * some of these models are changed, the offset frame moves away from the {@link #getBaseFrame() base frame}.
118      * </p>
119      *
120      * @param baseFrame     base frame associated with the station, without *any* parametric model (no station offset)
121      * @param clock         new quadratic clock model with user-supplied displacements
122      * @since 12.1
123      */
124     public GroundStation(final TopocentricFrame baseFrame, final QuadraticClockModel clock) {
125         super(baseFrame.getName(), clock);
126         this.baseFrame = baseFrame;
127 
128         this.eastOffsetDriver = new ParameterDriver(baseFrame.getName() + OFFSET_SUFFIX + "-East",
129                                                     0.0, POSITION_OFFSET_SCALE,
130                                                     Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
131 
132         this.northOffsetDriver = new ParameterDriver(baseFrame.getName() + OFFSET_SUFFIX + "-North",
133                                                      0.0, POSITION_OFFSET_SCALE,
134                                                      Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
135 
136         this.zenithOffsetDriver = new ParameterDriver(baseFrame.getName() + OFFSET_SUFFIX + "-Zenith",
137                                                       0.0, POSITION_OFFSET_SCALE,
138                                                       Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
139 
140         // Add the ground station parameters to the master list.
141         addParameterDriver(this.eastOffsetDriver);
142         addParameterDriver(this.northOffsetDriver);
143         addParameterDriver(this.zenithOffsetDriver);
144 
145     }
146 
147     /** {@inheritDoc} */
148     @Override
149     public final boolean isSpaceBased() {
150         return false;
151     }
152 
153     /** Get a driver allowing to change station position along East axis.
154      * @return driver for station position offset along East axis
155      */
156     public ParameterDriver getEastOffsetDriver() {
157         return eastOffsetDriver;
158     }
159 
160     /** Get a driver allowing to change station position along North axis.
161      * @return driver for station position offset along North axis
162      */
163     public ParameterDriver getNorthOffsetDriver() {
164         return northOffsetDriver;
165     }
166 
167     /** Get a driver allowing to change station position along Zenith axis.
168      * @return driver for station position offset along Zenith axis
169      */
170     public ParameterDriver getZenithOffsetDriver() {
171         return zenithOffsetDriver;
172     }
173 
174     /** Get the base frame associated with the station.
175      * <p>
176      * The base frame corresponds to a null position offset, null
177      * polar motion, null meridian shift
178      * </p>
179      * @return base frame associated with the station
180      */
181     public TopocentricFrame getBaseFrame() {
182         return baseFrame;
183     }
184 
185     /** Get the station displacement.
186      * @param date current date
187      * @param position raw position of the station in Earth frame
188      * before displacement is applied
189      * @return station displacement
190      * @since 9.1
191      */
192     protected Vector3D computeDisplacement(final AbsoluteDate date, final Vector3D position) {
193         return Vector3D.ZERO;
194     }
195 
196     /** Get the geodetic point at the center of the offset frame.
197      * @param date current date (may be null if displacements are ignored)
198      * @return geodetic point at the center of the offset frame
199      * @since 9.1
200      */
201     public GeodeticPoint getOffsetGeodeticPoint(final AbsoluteDate date) {
202 
203         // take station offset into account
204         final double    x          = eastOffsetDriver.getValue();
205         final double    y          = northOffsetDriver.getValue();
206         final double    z          = zenithOffsetDriver.getValue();
207         final BodyShape baseShape  = baseFrame.getParentShape();
208         final StaticTransform baseToBody = baseFrame.getStaticTransformTo(baseShape.getBodyFrame(), date);
209         Vector3D        origin     = baseToBody.transformPosition(new Vector3D(x, y, z));
210 
211         if (date != null) {
212             origin = origin.add(computeDisplacement(date, origin));
213         }
214 
215         return baseShape.transform(origin, baseShape.getBodyFrame(), date);
216 
217     }
218 
219     /** Get the geodetic point at the center of the offset frame.
220      * @param <T> type of the field elements
221      * @param date current date(<em>must</em> be non-null, which is a more stringent condition
222      *      *                    than in {@link #getOffsetGeodeticPoint(AbsoluteDate)}
223      * @return geodetic point at the center of the offset frame
224      * @since 12.1
225      */
226     public <T extends CalculusFieldElement<T>> FieldGeodeticPoint<T> getOffsetGeodeticPoint(final FieldAbsoluteDate<T> date) {
227 
228         // take station offset into account
229         final double    x          = eastOffsetDriver.getValue();
230         final double    y          = northOffsetDriver.getValue();
231         final double    z          = zenithOffsetDriver.getValue();
232         final BodyShape baseShape  = baseFrame.getParentShape();
233         final FieldStaticTransform<T> baseToBody = baseFrame.getStaticTransformTo(baseShape.getBodyFrame(), date);
234         FieldVector3D<T> origin    = baseToBody.transformPosition(new Vector3D(x, y, z));
235         origin = origin.add(computeDisplacement(date.toAbsoluteDate(), origin.toVector3D()));
236 
237         return baseShape.transform(origin, baseShape.getBodyFrame(), date);
238 
239     }
240 
241     /** {@inheritDoc} */
242     @Override
243     public PVCoordinatesProvider getPVCoordinatesProvider() {
244         final GeodeticPoint offsetPoint = getOffsetGeodeticPoint(AbsoluteDate.ARBITRARY_EPOCH);
245         return new TopocentricFrame(baseFrame.getParentShape(), offsetPoint, "offset");
246     }
247 
248     /** {@inheritDoc} */
249     @Override
250     public FieldPVCoordinatesProvider<Gradient> getFieldPVCoordinatesProvider(final int freeParameters,
251                                                                               final Map<String, Integer> parameterIndices) {
252         return new FieldPVCoordinatesProvider<>() {
253             @Override
254             public TimeStampedFieldPVCoordinates<Gradient> getPVCoordinates(final FieldAbsoluteDate<Gradient> date,
255                                                                             final Frame frame) {
256                 // take station offsets into account
257                 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
258 
259                 // body-fixed body-centered to target (with linear approximation for performance)
260                 final Transform bodyToInertNonField = baseFrame.getParent().getTransformTo(frame, date.toAbsoluteDate());
261                 final FieldTransform<Gradient> bodyToInert = new FieldTransform<>(date.getField(),
262                         bodyToInertNonField).shiftedBy(date.durationFrom(date.toAbsoluteDate()));
263 
264                 final TimeStampedFieldPVCoordinates<Gradient> zeroPV = new TimeStampedFieldPVCoordinates<>(date,
265                         new FieldPVCoordinates<>(origin, FieldVector3D.getZero(date.getField())));
266                 return bodyToInert.transformPVCoordinates(zeroPV);
267             }
268 
269             @Override
270             public FieldVector3D<Gradient> getPosition(final FieldAbsoluteDate<Gradient> date, final Frame frame) {
271                 // take station offsets into account
272                 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
273 
274                 // body-fixed body-centered to target (with linear approximation for performance)
275                 final KinematicTransform bodyToInertNonField = baseFrame.getParent().getKinematicTransformTo(frame,
276                         date.toAbsoluteDate());
277                 final FieldStaticTransform<Gradient> bodyToInert = shiftKinematicTransform(bodyToInertNonField,
278                         date.durationFrom(date.toAbsoluteDate()));
279 
280                 // combine by hand for performance reasons
281                 return bodyToInert.getRotation().applyTo(bodyToInert.getTranslation().add(origin));
282             }
283         };
284     }
285 
286     /**
287      * Retrieve station's position in body shape frame.
288      * @param date date
289      * @param indices mapping from parameters' name to derivatives' index.
290      * @return origin position
291      */
292     protected FieldVector3D<Gradient> getOrigin(final FieldAbsoluteDate<Gradient> date,
293                                                 final Map<String, Integer> indices) {
294         // compute position in topocentric frame
295         final int freeParameters = date.getField().getZero().getFreeParameters();
296         final AbsoluteDate absoluteDate = date.toAbsoluteDate();
297         final Gradient x          = eastOffsetDriver.getValue(freeParameters, indices, absoluteDate);
298         final Gradient                       y          = northOffsetDriver.getValue(freeParameters, indices, absoluteDate);
299         final Gradient                       z          = zenithOffsetDriver.getValue(freeParameters, indices, absoluteDate);
300         final FieldVector3D<Gradient> position = new FieldVector3D<>(x, y, z);
301         // approximate linearly (for performance) static transform from topocentric to body shape frame
302         final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
303         final KinematicTransform kinematicTopoToBody = baseFrame.getKinematicTransformTo(bodyFrame, absoluteDate);
304         final FieldStaticTransform<Gradient> staticTopoToBody = shiftKinematicTransform(kinematicTopoToBody,
305                 date.durationFrom(absoluteDate));
306         // apply transform and displacement
307         final FieldVector3D<Gradient>        originBeforeDisplacement     = staticTopoToBody.transformPosition(position);
308         return originBeforeDisplacement.add(computeDisplacement(absoluteDate, originBeforeDisplacement.toVector3D()));
309     }
310 
311     /**
312      * Shift a kinematic transform by a Gradient time into a FieldStaticTransform.
313      * @param kinematicTransform kinematic transform to shift
314      * @param dt time to shift by
315      * @return Field static transform shifted by dt
316      * @since 14.0
317      */
318     protected FieldStaticTransform<Gradient> shiftKinematicTransform(final KinematicTransform kinematicTransform,
319                                                                      final Gradient dt) {
320         // shift translation
321         final Field<Gradient> field = dt.getField();
322         final AbsoluteDate date = kinematicTransform.getDate();
323         final FieldVector3D<Gradient> fieldVelocity = new FieldVector3D<>(field, kinematicTransform.getVelocity());
324         final FieldVector3D<Gradient> shiftedTranslation = fieldVelocity.scalarMultiply(dt).add(kinematicTransform.getTranslation());
325         // shift rotation
326         final FieldAngularCoordinates<Gradient> fieldAngularCoordinates = new FieldAngularCoordinates<>(field,
327                 new AngularCoordinates(kinematicTransform.getRotation(), kinematicTransform.getRotationRate()));
328         final FieldVector3D<Gradient> rotationRate = fieldAngularCoordinates.getRotationRate();
329         final Gradient rate = rotationRate.getNorm();
330         final FieldRotation<Gradient> shiftedRotation = (rate.getReal() == 0.0) ?
331                 fieldAngularCoordinates.getRotation() :
332                 new FieldRotation<>(rotationRate, rate.multiply(dt), RotationConvention.FRAME_TRANSFORM)
333                         .compose(fieldAngularCoordinates.getRotation(), RotationConvention.VECTOR_OPERATOR);
334         return FieldStaticTransform.of(new FieldAbsoluteDate<>(field, date).shiftedBy(dt), shiftedTranslation,
335                 shiftedRotation);
336     }
337 
338     /** {@inheritDoc} */
339     @Override
340     public Transform getOffsetToInertial(final Frame inertial, final AbsoluteDate date,
341                                          final boolean clockOffsetAlreadyApplied) {
342 
343         // take clock offset into account
344         final AbsoluteDate offsetCompensatedDate = clockOffsetAlreadyApplied ?
345                 date :
346                 new AbsoluteDate(date, -getOffsetValue(date));
347 
348         final TopocentricFrame topocentricFrame = (TopocentricFrame) getPVCoordinatesProvider();
349         return topocentricFrame.getTransformTo(inertial, offsetCompensatedDate);
350     }
351 
352     /** {@inheritDoc} */
353     @Override
354     public FieldTransform<Gradient> getOffsetToInertial(final Frame inertial,
355                                                         final FieldAbsoluteDate<Gradient> offsetCompensatedDate,
356                                                         final int freeParameters,
357                                                         final Map<String, Integer> indices) {
358         // take station offsets into account
359         final FieldVector3D<Gradient> origin = getOrigin(offsetCompensatedDate, indices);
360         final FieldGeodeticPoint<Gradient> originGP = baseFrame.getParentShape().transform(origin, baseFrame.getParent(),
361                 offsetCompensatedDate);
362         final FieldStaticTransform<Gradient> staticOffsetToBody = TopocentricTransformProvider.getTransform(baseFrame.getParentShape(),
363                 offsetCompensatedDate, originGP).getStaticInverse();
364         final FieldTransform<Gradient> offsetToBody = new FieldTransform<>(offsetCompensatedDate,
365                 staticOffsetToBody.getTranslation(), staticOffsetToBody.getRotation());
366 
367         // Body-fixed, body-centered frame to target one
368         final FieldTransform<Gradient> bodyToInert = baseFrame.getParent().getTransformTo(inertial, offsetCompensatedDate);
369 
370         // combine all transforms together
371         return new FieldTransform<>(offsetCompensatedDate, offsetToBody, bodyToInert);
372     }
373 
374 }