1   /* Copyright 2002-2021 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.HashMap;
22  import java.util.Map;
23  
24  import org.hipparchus.analysis.differentiation.Gradient;
25  import org.hipparchus.analysis.differentiation.GradientField;
26  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27  import org.orekit.estimation.measurements.AbstractMeasurement;
28  import org.orekit.estimation.measurements.EstimatedMeasurement;
29  import org.orekit.estimation.measurements.GroundStation;
30  import org.orekit.estimation.measurements.ObservableSatellite;
31  import org.orekit.frames.FieldTransform;
32  import org.orekit.propagation.SpacecraftState;
33  import org.orekit.time.AbsoluteDate;
34  import org.orekit.time.FieldAbsoluteDate;
35  import org.orekit.utils.Constants;
36  import org.orekit.utils.ParameterDriver;
37  import org.orekit.utils.TimeStampedFieldPVCoordinates;
38  import org.orekit.utils.TimeStampedPVCoordinates;
39  
40  /** Class modeling a phase measurement from a ground station.
41   * <p>
42   * The measurement is considered to be a signal emitted from
43   * a spacecraft and received on a ground station.
44   * Its value is the number of cycles between emission and
45   * reception. The motion of both the station and the
46   * spacecraft during the signal flight time are taken into
47   * account. The date of the measurement corresponds to the
48   * reception on ground of the emitted signal.
49   * </p>
50   * @author Thierry Ceolin
51   * @author Luc Maisonobe
52   * @author Maxime Journot
53   * @since 9.2
54   */
55  public class Phase extends AbstractMeasurement<Phase> {
56  
57      /** Name for ambiguity driver. */
58      public static final String AMBIGUITY_NAME = "ambiguity";
59  
60      /** Driver for ambiguity. */
61      private final ParameterDriver ambiguityDriver;
62  
63      /** Ground station from which measurement is performed. */
64      private final GroundStation station;
65  
66      /** Wavelength of the phase observed value [m]. */
67      private final double wavelength;
68  
69      /** Simple constructor.
70       * @param station ground station from which measurement is performed
71       * @param date date of the measurement
72       * @param phase observed value (cycles)
73       * @param wavelength phase observed value wavelength (m)
74       * @param sigma theoretical standard deviation
75       * @param baseWeight base weight
76       * @param satellite satellite related to this measurement
77       * @since 9.3
78       */
79      public Phase(final GroundStation station, final AbsoluteDate date,
80                   final double phase, final double wavelength, final double sigma,
81                   final double baseWeight, final ObservableSatellite satellite) {
82          super(date, phase, sigma, baseWeight, Collections.singletonList(satellite));
83          ambiguityDriver = new ParameterDriver(AMBIGUITY_NAME,
84                                                 0.0, 1.0,
85                                                 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
86          addParameterDriver(ambiguityDriver);
87          addParameterDriver(satellite.getClockOffsetDriver());
88          addParameterDriver(station.getClockOffsetDriver());
89          addParameterDriver(station.getEastOffsetDriver());
90          addParameterDriver(station.getNorthOffsetDriver());
91          addParameterDriver(station.getZenithOffsetDriver());
92          addParameterDriver(station.getPrimeMeridianOffsetDriver());
93          addParameterDriver(station.getPrimeMeridianDriftDriver());
94          addParameterDriver(station.getPolarOffsetXDriver());
95          addParameterDriver(station.getPolarDriftXDriver());
96          addParameterDriver(station.getPolarOffsetYDriver());
97          addParameterDriver(station.getPolarDriftYDriver());
98          this.station    = station;
99          this.wavelength = wavelength;
100     }
101 
102     /** Get the ground station from which measurement is performed.
103      * @return ground station from which measurement is performed
104      */
105     public GroundStation getStation() {
106         return station;
107     }
108 
109     /** Get the wavelength.
110      * @return wavelength (m)
111      */
112     public double getWavelength() {
113         return wavelength;
114     }
115 
116     /** Get the driver for phase ambiguity.
117      * @return the driver for phase ambiguity
118      * @since 10.3
119      */
120     public ParameterDriver getAmbiguityDriver() {
121         return ambiguityDriver;
122     }
123 
124     /** {@inheritDoc} */
125     @Override
126     protected EstimatedMeasurement<Phase> theoreticalEvaluation(final int iteration,
127                                                                 final int evaluation,
128                                                                 final SpacecraftState[] states) {
129 
130         final SpacecraftState state = states[0];
131 
132         // Phase derivatives are computed with respect to spacecraft state in inertial frame
133         // and station parameters
134         // ----------------------
135         //
136         // Parameters:
137         //  - 0..2 - Position of the spacecraft in inertial frame
138         //  - 3..5 - Velocity of the spacecraft in inertial frame
139         //  - 6..n - station parameters (ambiguity, clock offset, station offsets, pole, prime meridian...)
140         int nbParams = 6;
141         final Map<String, Integer> indices = new HashMap<>();
142         for (ParameterDriver driver : getParametersDrivers()) {
143             if (driver.isSelected()) {
144                 indices.put(driver.getName(), nbParams++);
145             }
146         }
147         final FieldVector3D<Gradient> zero = FieldVector3D.getZero(GradientField.getField(nbParams));
148 
149         // Coordinates of the spacecraft expressed as a gradient
150         final TimeStampedFieldPVCoordinates<Gradient> pvaDS = getCoordinates(state, 0, nbParams);
151 
152         // transform between station and inertial frame, expressed as a gradient
153         // The components of station's position in offset frame are the 3 last derivative parameters
154         final FieldTransform<Gradient> offsetToInertialDownlink =
155                         station.getOffsetToInertial(state.getFrame(), getDate(), nbParams, indices);
156         final FieldAbsoluteDate<Gradient> downlinkDateDS =
157                         offsetToInertialDownlink.getFieldDate();
158 
159         // Station position in inertial frame at end of the downlink leg
160         final TimeStampedFieldPVCoordinates<Gradient> stationDownlink =
161                         offsetToInertialDownlink.transformPVCoordinates(new TimeStampedFieldPVCoordinates<>(downlinkDateDS,
162                                                                                                             zero, zero, zero));
163 
164         // Compute propagation times
165         // (if state has already been set up to pre-compensate propagation delay,
166         //  we will have delta == tauD and transitState will be the same as state)
167 
168         // Downlink delay
169         final Gradient tauD = signalTimeOfFlight(pvaDS, stationDownlink.getPosition(), downlinkDateDS);
170 
171         // Transit state & Transit state (re)computed with gradients
172         final Gradient        delta        = downlinkDateDS.durationFrom(state.getDate());
173         final Gradient        deltaMTauD   = tauD.negate().add(delta);
174         final SpacecraftState transitState = state.shiftedBy(deltaMTauD.getValue());
175         final TimeStampedFieldPVCoordinates<Gradient> transitStateDS = pvaDS.shiftedBy(deltaMTauD);
176 
177         // prepare the evaluation
178         final EstimatedMeasurement<Phase> estimated =
179                         new EstimatedMeasurement<Phase>(this, iteration, evaluation,
180                                                         new SpacecraftState[] {
181                                                             transitState
182                                                         }, new TimeStampedPVCoordinates[] {
183                                                             transitStateDS.toTimeStampedPVCoordinates(),
184                                                             stationDownlink.toTimeStampedPVCoordinates()
185                                                         });
186 
187         // Clock offsets
188         final ObservableSatellite satellite = getSatellites().get(0);
189         final Gradient            dts       = satellite.getClockOffsetDriver().getValue(nbParams, indices);
190         final Gradient            dtg       = station.getClockOffsetDriver().getValue(nbParams, indices);
191 
192         // Phase value
193         final double   cOverLambda = Constants.SPEED_OF_LIGHT / wavelength;
194         final Gradient ambiguity   = ambiguityDriver.getValue(nbParams, indices);
195         final Gradient phase       = tauD.add(dtg).subtract(dts).multiply(cOverLambda).add(ambiguity);
196 
197         estimated.setEstimatedValue(phase.getValue());
198 
199         // Phase partial derivatives with respect to state
200         final double[] derivatives = phase.getGradient();
201         estimated.setStateDerivatives(0, Arrays.copyOfRange(derivatives, 0, 6));
202 
203         // set partial derivatives with respect to parameters
204         // (beware element at index 0 is the value, not a derivative)
205         for (final ParameterDriver driver : getParametersDrivers()) {
206             final Integer index = indices.get(driver.getName());
207             if (index != null) {
208                 estimated.setParameterDerivatives(driver, derivatives[index]);
209             }
210         }
211 
212         return estimated;
213 
214     }
215 
216 }