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 GroundObserver {
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     @Override
148     public BodyShape getParentShape() {
149         return getBaseFrame().getParentShape();
150     }
151 
152     /** Get a driver allowing to change station position along East axis.
153      * @return driver for station position offset along East axis
154      */
155     public ParameterDriver getEastOffsetDriver() {
156         return eastOffsetDriver;
157     }
158 
159     /** Get a driver allowing to change station position along North axis.
160      * @return driver for station position offset along North axis
161      */
162     public ParameterDriver getNorthOffsetDriver() {
163         return northOffsetDriver;
164     }
165 
166     /** Get a driver allowing to change station position along Zenith axis.
167      * @return driver for station position offset along Zenith axis
168      */
169     public ParameterDriver getZenithOffsetDriver() {
170         return zenithOffsetDriver;
171     }
172 
173     /** Get the base frame associated with the station.
174      * <p>
175      * The base frame corresponds to a null position offset, null
176      * polar motion, null meridian shift
177      * </p>
178      * @return base frame associated with the station
179      */
180     public TopocentricFrame getBaseFrame() {
181         return baseFrame;
182     }
183 
184     /** Get the station displacement.
185      * @param date current date
186      * @param position raw position of the station in Earth frame
187      * before displacement is applied
188      * @return station displacement
189      * @since 9.1
190      */
191     protected Vector3D computeDisplacement(final AbsoluteDate date, final Vector3D position) {
192         return Vector3D.ZERO;
193     }
194 
195     /** Get the geodetic point at the center of the offset frame.
196      * @param date current date (may be null if displacements are ignored)
197      * @return geodetic point at the center of the offset frame
198      * @since 9.1
199      */
200     public GeodeticPoint getOffsetGeodeticPoint(final AbsoluteDate date) {
201 
202         // take station offset into account
203         final double    x          = eastOffsetDriver.getValue();
204         final double    y          = northOffsetDriver.getValue();
205         final double    z          = zenithOffsetDriver.getValue();
206         final BodyShape baseShape  = baseFrame.getParentShape();
207         final StaticTransform baseToBody = baseFrame.getStaticTransformTo(baseShape.getBodyFrame(), date);
208         Vector3D        origin     = baseToBody.transformPosition(new Vector3D(x, y, z));
209 
210         if (date != null) {
211             origin = origin.add(computeDisplacement(date, origin));
212         }
213 
214         return baseShape.transform(origin, baseShape.getBodyFrame(), date);
215 
216     }
217 
218     /** Get the geodetic point at the center of the offset frame.
219      * @param <T> type of the field elements
220      * @param date current date(<em>must</em> be non-null, which is a more stringent condition
221      *      *                    than in {@link #getOffsetGeodeticPoint(AbsoluteDate)}
222      * @return geodetic point at the center of the offset frame
223      * @since 12.1
224      */
225     public <T extends CalculusFieldElement<T>> FieldGeodeticPoint<T> getOffsetGeodeticPoint(final FieldAbsoluteDate<T> date) {
226 
227         // take station offset into account
228         final double    x          = eastOffsetDriver.getValue();
229         final double    y          = northOffsetDriver.getValue();
230         final double    z          = zenithOffsetDriver.getValue();
231         final BodyShape baseShape  = baseFrame.getParentShape();
232         final FieldStaticTransform<T> baseToBody = baseFrame.getStaticTransformTo(baseShape.getBodyFrame(), date);
233         FieldVector3D<T> origin    = baseToBody.transformPosition(new Vector3D(x, y, z));
234         origin = origin.add(computeDisplacement(date.toAbsoluteDate(), origin.toVector3D()));
235 
236         return baseShape.transform(origin, baseShape.getBodyFrame(), date);
237 
238     }
239 
240     /** {@inheritDoc} */
241     @Override
242     public PVCoordinatesProvider getPVCoordinatesProvider() {
243         final GeodeticPoint offsetPoint = getOffsetGeodeticPoint(AbsoluteDate.ARBITRARY_EPOCH);
244         return new TopocentricFrame(baseFrame.getParentShape(), offsetPoint, "offset");
245     }
246 
247     /** {@inheritDoc} */
248     @Override
249     public FieldPVCoordinatesProvider<Gradient> getFieldPVCoordinatesProvider(final int freeParameters,
250                                                                               final Map<String, Integer> parameterIndices) {
251         return new FieldPVCoordinatesProvider<>() {
252             @Override
253             public TimeStampedFieldPVCoordinates<Gradient> getPVCoordinates(final FieldAbsoluteDate<Gradient> date,
254                                                                             final Frame frame) {
255                 // take station offsets into account
256                 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
257 
258                 // body-fixed body-centered to target (with linear approximation for performance)
259                 final Transform bodyToInertNonField = baseFrame.getParent().getTransformTo(frame, date.toAbsoluteDate());
260                 final FieldTransform<Gradient> bodyToInert = new FieldTransform<>(date.getField(),
261                         bodyToInertNonField).shiftedBy(date.durationFrom(date.toAbsoluteDate()));
262 
263                 final TimeStampedFieldPVCoordinates<Gradient> zeroPV = new TimeStampedFieldPVCoordinates<>(date,
264                         new FieldPVCoordinates<>(origin, FieldVector3D.getZero(date.getField())));
265                 return bodyToInert.transformPVCoordinates(zeroPV);
266             }
267 
268             @Override
269             public FieldVector3D<Gradient> getPosition(final FieldAbsoluteDate<Gradient> date, final Frame frame) {
270                 // take station offsets into account
271                 final FieldVector3D<Gradient> origin = getOrigin(date, parameterIndices);
272 
273                 // body-fixed body-centered to target (with linear approximation for performance)
274                 final KinematicTransform bodyToInertNonField = baseFrame.getParent().getKinematicTransformTo(frame,
275                         date.toAbsoluteDate());
276                 final FieldStaticTransform<Gradient> bodyToInert = shiftKinematicTransform(bodyToInertNonField,
277                         date.durationFrom(date.toAbsoluteDate()));
278 
279                 // combine by hand for performance reasons
280                 return bodyToInert.getRotation().applyTo(bodyToInert.getTranslation().add(origin));
281             }
282         };
283     }
284 
285     /**
286      * Retrieve station's position in body shape frame.
287      * @param date date
288      * @param indices mapping from parameters' name to derivatives' index.
289      * @return origin position
290      */
291     protected FieldVector3D<Gradient> getOrigin(final FieldAbsoluteDate<Gradient> date,
292                                                 final Map<String, Integer> indices) {
293         // compute position in topocentric frame
294         final int freeParameters = date.getField().getZero().getFreeParameters();
295         final AbsoluteDate absoluteDate = date.toAbsoluteDate();
296         final Gradient x          = eastOffsetDriver.getValue(freeParameters, indices, absoluteDate);
297         final Gradient                       y          = northOffsetDriver.getValue(freeParameters, indices, absoluteDate);
298         final Gradient                       z          = zenithOffsetDriver.getValue(freeParameters, indices, absoluteDate);
299         final FieldVector3D<Gradient> position = new FieldVector3D<>(x, y, z);
300         // approximate linearly (for performance) static transform from topocentric to body shape frame
301         final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
302         final KinematicTransform kinematicTopoToBody = baseFrame.getKinematicTransformTo(bodyFrame, absoluteDate);
303         final FieldStaticTransform<Gradient> staticTopoToBody = shiftKinematicTransform(kinematicTopoToBody,
304                 date.durationFrom(absoluteDate));
305         // apply transform and displacement
306         final FieldVector3D<Gradient>        originBeforeDisplacement     = staticTopoToBody.transformPosition(position);
307         return originBeforeDisplacement.add(computeDisplacement(absoluteDate, originBeforeDisplacement.toVector3D()));
308     }
309 
310     /**
311      * Shift a kinematic transform by a Gradient time into a FieldStaticTransform.
312      * @param kinematicTransform kinematic transform to shift
313      * @param dt time to shift by
314      * @return Field static transform shifted by dt
315      * @since 14.0
316      */
317     protected FieldStaticTransform<Gradient> shiftKinematicTransform(final KinematicTransform kinematicTransform,
318                                                                      final Gradient dt) {
319         // shift translation
320         final Field<Gradient> field = dt.getField();
321         final AbsoluteDate date = kinematicTransform.getDate();
322         final FieldVector3D<Gradient> fieldVelocity = new FieldVector3D<>(field, kinematicTransform.getVelocity());
323         final FieldVector3D<Gradient> shiftedTranslation = fieldVelocity.scalarMultiply(dt).add(kinematicTransform.getTranslation());
324         // shift rotation
325         final FieldAngularCoordinates<Gradient> fieldAngularCoordinates = new FieldAngularCoordinates<>(field,
326                 new AngularCoordinates(kinematicTransform.getRotation(), kinematicTransform.getRotationRate()));
327         final FieldVector3D<Gradient> rotationRate = fieldAngularCoordinates.getRotationRate();
328         final Gradient rate = rotationRate.getNorm();
329         final FieldRotation<Gradient> shiftedRotation = (rate.getReal() == 0.0) ?
330                 fieldAngularCoordinates.getRotation() :
331                 new FieldRotation<>(rotationRate, rate.multiply(dt), RotationConvention.FRAME_TRANSFORM)
332                         .compose(fieldAngularCoordinates.getRotation(), RotationConvention.VECTOR_OPERATOR);
333         return FieldStaticTransform.of(new FieldAbsoluteDate<>(field, date).shiftedBy(dt), shiftedTranslation,
334                 shiftedRotation);
335     }
336 
337     /** {@inheritDoc} */
338     @Override
339     public Transform getOffsetToInertial(final Frame inertial, final AbsoluteDate date,
340                                          final boolean clockOffsetAlreadyApplied) {
341 
342         // take clock offset into account
343         final AbsoluteDate offsetCompensatedDate = clockOffsetAlreadyApplied ?
344                 date :
345                 new AbsoluteDate(date, -getOffsetValue(date));
346 
347         final TopocentricFrame topocentricFrame = (TopocentricFrame) getPVCoordinatesProvider();
348         return topocentricFrame.getTransformTo(inertial, offsetCompensatedDate);
349     }
350 
351     /** {@inheritDoc} */
352     @Override
353     public FieldTransform<Gradient> getOffsetToInertial(final Frame inertial,
354                                                         final FieldAbsoluteDate<Gradient> offsetCompensatedDate,
355                                                         final int freeParameters,
356                                                         final Map<String, Integer> indices) {
357         // take station offsets into account
358         final FieldVector3D<Gradient> origin = getOrigin(offsetCompensatedDate, indices);
359         final FieldGeodeticPoint<Gradient> originGP = baseFrame.getParentShape().transform(origin, baseFrame.getParent(),
360                 offsetCompensatedDate);
361         final FieldStaticTransform<Gradient> staticOffsetToBody = TopocentricTransformProvider.getTransform(baseFrame.getParentShape(),
362                 offsetCompensatedDate, originGP).getStaticInverse();
363         final FieldTransform<Gradient> offsetToBody = new FieldTransform<>(offsetCompensatedDate,
364                 staticOffsetToBody.getTranslation(), staticOffsetToBody.getRotation());
365 
366         // Body-fixed, body-centered frame to target one
367         final FieldTransform<Gradient> bodyToInert = baseFrame.getParent().getTransformTo(inertial, offsetCompensatedDate);
368 
369         // combine all transforms together
370         return new FieldTransform<>(offsetCompensatedDate, offsetToBody, bodyToInert);
371     }
372 
373 }