1   /* Copyright 2022-2026 Romain Serra
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 org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
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.bodies.BodyShape;
26  import org.orekit.bodies.FieldGeodeticPoint;
27  import org.orekit.bodies.GeodeticPoint;
28  import org.orekit.data.BodiesElements;
29  import org.orekit.data.FundamentalNutationArguments;
30  import org.orekit.frames.FieldStaticTransform;
31  import org.orekit.frames.FieldTransform;
32  import org.orekit.frames.Frame;
33  import org.orekit.frames.StaticTransform;
34  import org.orekit.frames.TopocentricFrame;
35  import org.orekit.frames.Transform;
36  import org.orekit.frames.TransformProvider;
37  import org.orekit.models.earth.displacement.StationDisplacement;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.time.FieldAbsoluteDate;
40  import org.orekit.utils.ParameterDriver;
41  
42  /** Class modeling a Earth-based station frame transform.
43   * <p>
44   * This class considers a position offset parameter w.r.t. a base {@link TopocentricFrame
45   * topocentric frame}.
46   * </p>
47   * <p>
48   * This class also adds parameters for an additional polar motion
49   * and an additional prime meridian orientation. Since these parameters will
50   * have the same name for all ground stations, they will be managed consistently
51   * and allow to estimate Earth orientation precisely (this is needed for precise
52   * orbit determination). The polar motion and prime meridian orientation will
53   * be applied <em>after</em> regular Earth orientation parameters, so the value
54   * of the estimated parameters will be correction to EOP, they will not be the
55   * complete EOP values by themselves. Basically, this means that for Earth, the
56   * following transforms are applied in order, between inertial frame and ground
57   * station frame (for non-Earth based ground stations, different precession nutation
58   * models and associated planet orientation parameters would be applied, if available):
59   * </p>
60   * @author Romain Serra
61   * @see EarthBasedStation
62   * @since 14.0
63   */
64  class EarthBasedStationTransformProvider implements TransformProvider {
65  
66      /**
67       * Target frame.
68       */
69      private final Frame frame;
70  
71      /** Provider for Earth frame whose EOP parameters can be estimated. */
72      private final EstimatedEarthFrameProvider estimatedEarthFrameProvider;
73  
74      /** Earth frame whose EOP parameters can be estimated. */
75      private final Frame estimatedEarthFrame;
76  
77      /** Base frame associated with the station. */
78      private final TopocentricFrame baseFrame;
79  
80      /** Fundamental nutation arguments. */
81      private final FundamentalNutationArguments arguments;
82  
83      /** Displacement models. */
84      private final StationDisplacement[] displacements;
85  
86      /** Driver for position offset along the East axis. */
87      private final ParameterDriver eastOffsetDriver;
88  
89      /** Driver for position offset along the North axis. */
90      private final ParameterDriver northOffsetDriver;
91  
92      /** Driver for position offset along the zenith axis. */
93      private final ParameterDriver zenithOffsetDriver;
94  
95       /**
96       * Constructor.
97       * @param frame target frame
98       * @param baseFrame     base frame associated with the station, without *any* parametric model (no station offset,
99       *                      no polar motion, no meridian shift)
100      * @param eastOffsetDriver driver for position offset along the East axis
101      * @param northOffsetDriver driver for position offset along the North axis
102      * @param zenithOffsetDriver driver for position offset along the zenith axis
103      * @param estimatedEarthFrameProvider provider for Earth frame whose EOP parameters can be estimated
104      * @param fundamentalNutationArguments fundamental nutation arguments
105      * @param displacements ground station displacement model (tides, ocean loading, atmospheric loading, thermal
106      *                      effects...)
107      */
108     EarthBasedStationTransformProvider(final Frame frame, final TopocentricFrame baseFrame,
109                                        final ParameterDriver eastOffsetDriver,
110                                        final ParameterDriver northOffsetDriver,
111                                        final ParameterDriver zenithOffsetDriver,
112                                        final EstimatedEarthFrameProvider estimatedEarthFrameProvider,
113                                        final FundamentalNutationArguments fundamentalNutationArguments,
114                                        final StationDisplacement... displacements) {
115         this.frame = frame;
116         this.baseFrame = baseFrame;
117 
118         this.estimatedEarthFrameProvider = estimatedEarthFrameProvider;
119         this.estimatedEarthFrame = new Frame(baseFrame.getParent(), estimatedEarthFrameProvider,
120                                              baseFrame.getParent() + "-estimated");
121         this.arguments = fundamentalNutationArguments;
122 
123         this.displacements = displacements.clone();
124         this.eastOffsetDriver   = eastOffsetDriver;
125         this.northOffsetDriver  = northOffsetDriver;
126         this.zenithOffsetDriver = zenithOffsetDriver;
127     }
128 
129     /** {@inheritDoc} */
130     @Override
131     public Transform getTransform(final AbsoluteDate date) {
132         // take Earth offsets into account
133         final Transform intermediateToBody = estimatedEarthFrameProvider.getTransform(date).getInverse();
134 
135         // take station offsets into account
136         final BodyShape baseShape = baseFrame.getParentShape();
137         final Vector3D  origin    = getOrigin(date);
138 
139         final GeodeticPoint originGP = baseShape.transform(origin, baseShape.getBodyFrame(), date);
140         final Transform offsetToIntermediate =
141                         new Transform(date,
142                                       new Transform(date,
143                                                     new Rotation(Vector3D.PLUS_I, Vector3D.PLUS_K,
144                                                                  originGP.getEast(), originGP.getZenith()),
145                                                     Vector3D.ZERO),
146                                       new Transform(date, origin));
147         if (baseFrame.getParent() == frame) {
148             return new Transform(date, offsetToIntermediate, intermediateToBody);
149         }
150 
151         // combine all transforms together
152         final Transform bodyToInert = baseFrame.getParent().getTransformTo(frame, date);
153 
154         return new Transform(date, offsetToIntermediate, new Transform(date, intermediateToBody, bodyToInert));
155     }
156 
157     /** {@inheritDoc} */
158     @Override
159     public StaticTransform getStaticTransform(final AbsoluteDate date) {
160         // take Earth offsets into account
161         final StaticTransform intermediateToBody = estimatedEarthFrameProvider.getStaticTransform(date).getInverse();
162 
163         // take station offsets into account
164         final BodyShape baseShape = baseFrame.getParentShape();
165         final Vector3D  origin    = getOrigin(date);
166 
167         final GeodeticPoint originGP = baseShape.transform(origin, baseShape.getBodyFrame(), date);
168         final StaticTransform offsetToIntermediate = StaticTransform.compose(date,
169                         StaticTransform.of(date,
170                                 new Rotation(Vector3D.PLUS_I, Vector3D.PLUS_K,
171                                         originGP.getEast(), originGP.getZenith())),
172                         StaticTransform.of(date, origin));
173         if (baseFrame.getParent() == frame) {
174             return StaticTransform.compose(date, offsetToIntermediate, intermediateToBody);
175         }
176 
177         // combine all transforms together
178         final StaticTransform bodyToInert = baseFrame.getParent().getStaticTransformTo(frame, date);
179 
180         return StaticTransform.compose(date, offsetToIntermediate, StaticTransform.compose(date, intermediateToBody, bodyToInert));
181     }
182 
183     /**
184      * Retrieve station's position in body shape frame.
185      * @param date date
186      * @return origin position
187      */
188     private Vector3D getOrigin(final AbsoluteDate date) {
189         final double    x          = eastOffsetDriver.getValue(date);
190         final double    y          = northOffsetDriver.getValue(date);
191         final double    z          = zenithOffsetDriver.getValue(date);
192         final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
193         final StaticTransform staticTopoToBody = baseFrame.getStaticTransformTo(bodyFrame, date);
194         final Vector3D        originBeforeDisplacement     = staticTopoToBody.transformPosition(new Vector3D(x, y, z));
195         return originBeforeDisplacement.add(computeDisplacement(date, originBeforeDisplacement));
196     }
197 
198     /** Get the station displacement.
199      * @param date current date
200      * @param position raw position of the station in Earth frame
201      * before displacement is applied
202      * @return station displacement
203      */
204     private Vector3D computeDisplacement(final AbsoluteDate date, final Vector3D position) {
205         Vector3D displacement = Vector3D.ZERO;
206         if (arguments != null) {
207             final BodiesElements elements = arguments.evaluateAll(date);
208             for (final StationDisplacement sd : displacements) {
209                 // we consider all displacements apply to the same initial position,
210                 // i.e. they apply simultaneously, not according to some order
211                 displacement = displacement.add(sd.displacement(elements, estimatedEarthFrame, position));
212             }
213         }
214         return displacement;
215     }
216 
217     /** {@inheritDoc} */
218     @Override
219     public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
220 
221         // take Earth offsets into account
222         final FieldTransform<T> intermediateToBody = estimatedEarthFrameProvider.getTransform(date).getInverse();
223 
224         // take station offsets into account
225         final FieldVector3D<T> origin = getOrigin(date);
226         return getTransform(date, origin, intermediateToBody);
227 
228     }
229 
230     <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date,
231                                                                        final FieldVector3D<T> origin,
232                                                                        final FieldTransform<T> intermediateToBody) {
233         final Field<T>         field = date.getField();
234         final FieldVector3D<T> zero  = FieldVector3D.getZero(field);
235         final FieldVector3D<T> plusI = FieldVector3D.getPlusI(field);
236         final FieldVector3D<T> plusK = FieldVector3D.getPlusK(field);
237         final FieldGeodeticPoint<T> originGP = baseFrame.getParentShape().transform(origin,
238                 baseFrame.getParentShape().getBodyFrame(), date);
239         final FieldTransform<T> offsetToIntermediate =
240                 new FieldTransform<>(date,
241                         new FieldTransform<>(date,
242                                 new FieldRotation<>(plusI, plusK,
243                                         originGP.getEast(), originGP.getZenith()),
244                                 zero),
245                         new FieldTransform<>(date, origin));
246 
247         // combine all transforms together
248         if (baseFrame.getParent() == frame) {
249             return new FieldTransform<>(date, offsetToIntermediate, intermediateToBody);
250         }
251         final FieldTransform<T> bodyToInert = baseFrame.getParent().getTransformTo(frame, date);
252 
253         return new FieldTransform<>(date, offsetToIntermediate,
254                 new FieldTransform<>(date, intermediateToBody, bodyToInert));
255 
256     }
257 
258     /** {@inheritDoc} */
259     @Override
260     public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
261 
262         // take Earth offsets into account
263         final FieldStaticTransform<T> intermediateToBody = estimatedEarthFrameProvider.getStaticTransform(date).getInverse();
264 
265         final FieldVector3D<T> origin = new FieldVector3D<>(date.getField(), getOrigin(date.toAbsoluteDate()));
266         // take station offsets into account
267         final Field<T>         field = date.getField();
268         final FieldVector3D<T> plusI = FieldVector3D.getPlusI(field);
269         final FieldVector3D<T> plusK = FieldVector3D.getPlusK(field);
270         final FieldGeodeticPoint<T> originGP = baseFrame.getParentShape().transform(origin,
271                 baseFrame.getParentShape().getBodyFrame(), date);
272         final FieldStaticTransform<T> offsetToIntermediate = FieldStaticTransform.compose(date,
273                 FieldStaticTransform.of(date,
274                         new FieldRotation<>(plusI, plusK, originGP.getEast(), originGP.getZenith())),
275                 FieldStaticTransform.of(date, origin));
276 
277         // combine all transforms together
278         if (baseFrame.getParent() == frame) {
279             return FieldStaticTransform.compose(date, offsetToIntermediate, intermediateToBody);
280         }
281         final FieldStaticTransform<T> bodyToInert = baseFrame.getParent().getStaticTransformTo(frame, date);
282 
283         return FieldStaticTransform.compose(date, offsetToIntermediate,
284                 FieldStaticTransform.compose(date, intermediateToBody, bodyToInert));
285     }
286 
287     /**
288      * Retrieve station's position in body shape frame.
289      * @param <T> field type
290      * @param date date
291      * @return origin position
292      */
293     private <T extends CalculusFieldElement<T>> FieldVector3D<T> getOrigin(final FieldAbsoluteDate<T> date) {
294         final AbsoluteDate absoluteDate = date.toAbsoluteDate();
295         final Field<T> field = date.getField();
296         final T x          = field.getZero().newInstance(eastOffsetDriver.getValue(absoluteDate));
297         final T                       y          = field.getZero().newInstance(northOffsetDriver.getValue(absoluteDate));
298         final T                       z          = field.getZero().newInstance(zenithOffsetDriver.getValue(absoluteDate));
299         final Frame bodyFrame = baseFrame.getParentShape().getBodyFrame();
300         final FieldStaticTransform<T> staticTopoToBody = baseFrame.getStaticTransformTo(bodyFrame, date);
301         final FieldVector3D<T>        originBeforeDisplacement     = staticTopoToBody.transformPosition(new FieldVector3D<>(x, y, z));
302         return originBeforeDisplacement.add(computeDisplacement(absoluteDate, originBeforeDisplacement.toVector3D()));
303     }
304 
305 }