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.gnss;
18  
19  import java.util.Arrays;
20  import java.util.Collections;
21  import java.util.Map;
22  
23  import org.hipparchus.analysis.differentiation.Gradient;
24  import org.hipparchus.analysis.differentiation.GradientField;
25  import org.orekit.estimation.measurements.AbstractMeasurement;
26  import org.orekit.estimation.measurements.AbstractParticipant;
27  import org.orekit.estimation.measurements.EstimatedMeasurement;
28  import org.orekit.estimation.measurements.EstimatedMeasurementBase;
29  import org.orekit.estimation.measurements.MeasurementQuality;
30  import org.orekit.estimation.measurements.ObservableSatellite;
31  import org.orekit.estimation.measurements.Observer;
32  import org.orekit.estimation.measurements.SignalBasedMeasurement;
33  import org.orekit.frames.FieldTransform;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.Transform;
36  import org.orekit.propagation.SpacecraftState;
37  import org.orekit.signal.AdjustableEmitterSignalTimer;
38  import org.orekit.signal.FieldAdjustableEmitterSignalTimer;
39  import org.orekit.signal.FieldSignalReceptionCondition;
40  import org.orekit.signal.SignalReceptionCondition;
41  import org.orekit.signal.SignalTravelTimeModel;
42  import org.orekit.time.AbsoluteDate;
43  import org.orekit.time.FieldAbsoluteDate;
44  import org.orekit.utils.Constants;
45  import org.orekit.utils.FieldPVCoordinates;
46  import org.orekit.utils.FieldPVCoordinatesProvider;
47  import org.orekit.utils.PVCoordinates;
48  import org.orekit.utils.PVCoordinatesProvider;
49  import org.orekit.utils.drivers.ParameterDriver;
50  import org.orekit.utils.TimeStampedFieldPVCoordinates;
51  import org.orekit.utils.TimeStampedPVCoordinates;
52  
53  /** Class modeling a phase measurement from a ground station.
54   * <p>
55   * The measurement is considered to be a signal emitted from
56   * a spacecraft and received on a ground station.
57   * Its value is the number of cycles between emission and
58   * reception. The motion of both the station and the
59   * spacecraft during the signal flight time are taken into
60   * account. The date of the measurement corresponds to the
61   * reception on ground of the emitted signal.
62   * </p>
63   * @author Thierry Ceolin
64   * @author Luc Maisonobe
65   * @author Maxime Journot
66   * @since 9.2
67   */
68  public class Phase extends SignalBasedMeasurement<Phase> {
69  
70      /** Type of the measurement. */
71      public static final String MEASUREMENT_TYPE = "Phase";
72  
73      /** Driver for ambiguity. */
74      private final AmbiguityDriver ambiguityDriver;
75  
76      /** Wavelength of the phase observed value [m]. */
77      private final double wavelength;
78  
79      /** Observer that receives signal from satellite. */
80      private final Observer observer;
81  
82      /** Simple constructor.
83       * @param observer observer that performs the measurement
84       * @param date date of the measurement
85       * @param phase observed value (cycles)
86       * @param wavelength phase observed value wavelength (m)
87       * @param sigma theoretical standard deviation
88       * @param baseWeight base weight
89       * @param satellite satellite related to this measurement
90       * @param cache from which ambiguity drive should come
91       * @since 12.1
92       */
93      public Phase(final Observer observer, final AbsoluteDate date,
94                   final double phase, final double wavelength, final double sigma,
95                   final double baseWeight, final ObservableSatellite satellite,
96                   final AmbiguityCache cache) {
97          this(observer, date, phase, wavelength, new MeasurementQuality(sigma, baseWeight), new SignalTravelTimeModel(),
98                  satellite, cache);
99      }
100 
101     /** Simple constructor.
102      * @param observer observer that performs the measurement
103      * @param date date of the measurement
104      * @param phase observed value (cycles)
105      * @param wavelength phase observed value wavelength (m)
106      * @param measurementQuality measurement quality data as used in orbit determination
107      * @param signalTravelTimeModel signal model
108      * @param satellite satellite related to this measurement
109      * @param cache from which ambiguity drive should come
110      * @since 14.0
111      */
112     public Phase(final Observer observer, final AbsoluteDate date,
113                  final double phase, final double wavelength, final MeasurementQuality measurementQuality,
114                  final SignalTravelTimeModel signalTravelTimeModel, final ObservableSatellite satellite,
115                  final AmbiguityCache cache) {
116         super(date, false, phase, measurementQuality, signalTravelTimeModel,
117                 Collections.singletonList(satellite));
118         ambiguityDriver = cache.getAmbiguity(satellite.getName(), observer.getName(), wavelength);
119         addParametersDrivers(observer.getParametersDrivers());
120         addParameterDriver(ambiguityDriver);
121         this.observer = observer;
122         this.wavelength = wavelength;
123     }
124 
125     /** Get receiving observer.
126      * @return observer
127      */
128     public final Observer getObserver() {
129         return observer;
130     }
131 
132     /** Get the wavelength.
133      * @return wavelength (m)
134      */
135     public double getWavelength() {
136         return wavelength;
137     }
138 
139     /** Get the driver for phase ambiguity.
140      * @return the driver for phase ambiguity
141      * @since 10.3
142      */
143     public AmbiguityDriver getAmbiguityDriver() {
144         return ambiguityDriver;
145     }
146 
147     /** {@inheritDoc} */
148     @Override
149     protected EstimatedMeasurementBase<Phase> theoreticalEvaluationWithoutDerivatives(final int iteration,
150                                                                                       final int evaluation,
151                                                                                       final SpacecraftState[] states,
152                                                                                       final boolean fillParticipants) {
153         // Coordinates of the measured spacecraft
154         final SpacecraftState state = states[0];
155         final Frame frame = state.getFrame();
156         final TimeStampedPVCoordinates pva   = state.getPVCoordinates();
157 
158         // transform between remote observer frame and inertial frame
159         final AbsoluteDate measurementDate = getDate();
160         final Transform offsetToInertialDownlink = getObserver().getOffsetToInertial(frame, measurementDate, false);
161         final AbsoluteDate downlinkDate             = offsetToInertialDownlink.getDate();
162 
163         // Observer position in inertial frame at end of the downlink leg
164         final TimeStampedPVCoordinates origin = new TimeStampedPVCoordinates(downlinkDate, PVCoordinates.ZERO);
165         final TimeStampedPVCoordinates satelliteDownlink = offsetToInertialDownlink.transformPVCoordinates(origin);
166 
167         // Coordinates provider for emitting object (observed spacecraft)
168         final PVCoordinatesProvider pvCoordinatesProvider = AbstractParticipant.extractPVCoordinatesProvider(states[0], pva);
169 
170         // Downlink delay / determine time of emission of signal by ObservableSatellite
171         final AdjustableEmitterSignalTimer signalTimeOfFlight = getSignalTravelTimeModel()
172                 .getAdjustableEmitterComputer(pvCoordinatesProvider);
173         final SignalReceptionCondition receptionCondition = new SignalReceptionCondition(downlinkDate,
174                 satelliteDownlink.getPosition(), frame);
175         final double tauD = signalTimeOfFlight.computeDelay(receptionCondition, pva.getDate());
176 
177         // Transit state & Transit state (re)computed with gradients
178         final double          delta             = downlinkDate.durationFrom(state.getDate());
179         final double          deltaMTauD        = delta - tauD;
180         final SpacecraftState transitState      = states[0].shiftedBy(deltaMTauD);
181 
182         // prepare the evaluation
183         final EstimatedMeasurementBase<Phase> estimated = new EstimatedMeasurementBase<>(this, iteration, evaluation,
184                                                        new SpacecraftState[] { transitState }, fillParticipants ? new TimeStampedPVCoordinates[] {
185                                                            transitState.getPVCoordinates(), satelliteDownlink } : new TimeStampedPVCoordinates[0]);
186 
187         // Clock offsets
188         final ObservableSatellite satellite = getSatellites().getFirst();
189 
190         final double dts = satellite.getOffsetValue(state.getDate());
191         final double dtg = getObserver().getOffsetValue(getDate());
192 
193         // Phase value
194         final double cOverLambda = Constants.SPEED_OF_LIGHT / wavelength;
195         final double ambiguity   = ambiguityDriver.getValue();
196         final double phase       = (tauD + dtg - dts) * cOverLambda + ambiguity;
197 
198         estimated.setEstimatedValue(phase);
199 
200         return estimated;
201 
202     }
203 
204     /** {@inheritDoc} */
205     @Override
206     protected EstimatedMeasurement<Phase> theoreticalEvaluation(final int iteration,
207                                                                 final int evaluation,
208                                                                 final SpacecraftState[] states) {
209         // Create the parameter indices map
210         final SpacecraftState state = states[0];
211         final Frame                frame        = state.getFrame();
212         final Map<String, Integer> paramIndices = getParameterIndices(states);
213         final int                  nbParams     = 6 * states.length + paramIndices.size();
214 
215         // Coordinates of the spacecraft expressed as a gradient
216         final TimeStampedFieldPVCoordinates<Gradient> pva = AbstractMeasurement.getCoordinates(states[0], 0, nbParams);
217 
218         // transform between Observer object and inertial frame, expressed as a gradient
219         // The components of the Observer's position in offset frame are the 3 last derivative parameters
220         final FieldTransform<Gradient> offsetToInertialDownlink = getObserver().
221                 getOffsetToInertial(frame, getDate(), nbParams, paramIndices);
222         final FieldAbsoluteDate<Gradient> downlinkDate = offsetToInertialDownlink.getFieldDate();
223 
224         // Observer position in inertial frame at end of the downlink leg
225         final GradientField field = GradientField.getField(nbParams);
226         final TimeStampedFieldPVCoordinates<Gradient> satelliteDownlink =
227                 offsetToInertialDownlink.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(downlinkDate,
228                         FieldPVCoordinates.getZero(field)));
229 
230         // Form coordinates provider
231         final FieldPVCoordinatesProvider<Gradient> fieldPVCoordinatesProvider = AbstractParticipant.extractFieldPVCoordinatesProvider(states[0], pva);
232 
233         // Downlink delay
234         final FieldAdjustableEmitterSignalTimer<Gradient> fieldComputer = getSignalTravelTimeModel()
235                 .getFieldAdjustableEmitterComputer(field, fieldPVCoordinatesProvider);
236         final FieldSignalReceptionCondition<Gradient> receptionCondition = new FieldSignalReceptionCondition<>(downlinkDate,
237                 satelliteDownlink.getPosition(), frame);
238         final Gradient tauD = fieldComputer.computeDelay(receptionCondition, pva.getDate());
239 
240         // Transit state & Transit state (re)computed with gradients
241         final Gradient        delta        = downlinkDate.durationFrom(states[0].getDate());
242         final Gradient        deltaMTauD   = tauD.negate().add(delta);
243         final SpacecraftState transitState = states[0].shiftedBy(deltaMTauD.getValue());
244         final FieldAbsoluteDate<Gradient> fieldDate = new FieldAbsoluteDate<>(field, states[0].getDate()).shiftedBy(deltaMTauD);
245         final TimeStampedFieldPVCoordinates<Gradient> transitPV = fieldPVCoordinatesProvider.getPVCoordinates(fieldDate, frame);
246 
247         // prepare the evaluation
248         final EstimatedMeasurement<Phase> estimated =
249                         new EstimatedMeasurement<>(this, iteration, evaluation,
250                                 new SpacecraftState[] { transitState},
251                                 new TimeStampedPVCoordinates[] {
252                                         transitPV.toTimeStampedPVCoordinates(),
253                                         satelliteDownlink.toTimeStampedPVCoordinates()});
254 
255         // Clock offsets
256         final ObservableSatellite satellite = getSatellites().getFirst();
257 
258         final Gradient dts = satellite.getFieldOffsetValue(nbParams, paramIndices, state.getDate());
259         final Gradient dtg = getObserver().getFieldOffsetValue(nbParams, paramIndices, getDate());
260 
261         // Phase value
262         final double   cOverLambda = Constants.SPEED_OF_LIGHT / wavelength;
263         final Gradient ambiguity   = ambiguityDriver.getValue(nbParams, paramIndices);
264         final Gradient phase       = tauD.add(dtg).subtract(dts).multiply(cOverLambda).add(ambiguity);
265 
266         estimated.setEstimatedValue(phase.getValue());
267 
268         // Phase first order derivatives with respect to state
269         final double[] derivatives = phase.getGradient();
270         estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 0, 6));
271 
272         // Set first order derivatives with respect to parameters
273         for (final ParameterDriver driver : getParametersDrivers()) {
274             final Integer index = paramIndices.get(driver.getName());
275             if (index != null) {
276                 estimated.setParameterDerivatives(driver, derivatives[index]);
277             }
278         }
279 
280         return estimated;
281 
282     }
283 
284 }