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.ArrayList;
20  import java.util.List;
21  import java.util.Map;
22  
23  import org.hipparchus.fitting.PolynomialCurveFitter;
24  import org.hipparchus.fitting.WeightedObservedPoint;
25  import org.hipparchus.util.FastMath;
26  import org.orekit.gnss.CombinedObservationData;
27  import org.orekit.gnss.CombinedObservationDataSet;
28  import org.orekit.gnss.Frequency;
29  import org.orekit.gnss.MeasurementType;
30  import org.orekit.gnss.ObservationDataSet;
31  import org.orekit.gnss.SatelliteSystem;
32  import org.orekit.time.AbsoluteDate;
33  
34  
35  /**
36   * Geometry free cycle slip detectors.
37   * The detector is based the algorithm given in <a
38   * href="https://gssc.esa.int/navipedia/index.php/Detector_based_in_carrier_phase_data:_The_geometry-free_combination">
39   * Detector based in carrier phase data: The geometry-free combination</a> by Zornoza and M. Hernández-Pajares. Within this class
40   * a second order polynomial is used to smooth the data. We consider a cycle-slip occurring if the current measurement is  too
41   * far from the one predicted with the polynomial.
42   * <p>
43   * For building the detector, one should give a threshold and a gap time limit.
44   * After construction of the detectors, one can have access to a List of CycleData. Each CycleDate represents
45   * a link between the station (define by the RINEX file) and a satellite at a specific frequency. For each cycle data,
46   * one has access to the begin and end of availability, and a sorted set which contains all the date at which
47   * cycle-slip have been detected
48   * </p>
49   * <p>
50   * @author David Soulard
51   * @author Bryan Cazabonne
52   * @since 10.2
53   */
54  public class GeometryFreeCycleSlipDetector extends AbstractCycleSlipDetector {
55  
56      /** Threshold above which cycle-slip occurs. */
57      private final double threshold;
58  
59      /**
60       * Constructor.
61       * @param dt time gap threshold between two consecutive measurement (if time between two consecutive measurement is greater than dt, a cycle slip is declared)
62       * @param threshold threshold above which cycle-slip occurs
63       * @param n number of measurement before starting
64       */
65      public GeometryFreeCycleSlipDetector(final double dt, final double threshold, final int n) {
66          super(dt, n);
67          this.threshold = threshold;
68      }
69  
70  
71      /** {@inheritDoc} */
72      @Override
73      protected void manageData(final ObservationDataSet observation) {
74  
75          // Extract observation data
76          final int             prn    = observation.getPrnNumber();
77          final AbsoluteDate    date   = observation.getDate();
78          final SatelliteSystem system = observation.getSatelliteSystem();
79  
80          // Geometry-free combination of measurements
81          final GeometryFreeCombination geometryFree = MeasurementCombinationFactory.getGeometryFreeCombination(system);
82          final CombinedObservationDataSet cods = geometryFree.combine(observation);
83  
84          // Initialize list of measurements
85          final List<CombinedObservationData> phasesGF = new ArrayList<>();
86  
87          // Loop on observation data to fill lists
88          for (final CombinedObservationData cod : cods.getObservationData()) {
89              if (!Double.isNaN(cod.getValue()) && cod.getMeasurementType() == MeasurementType.CARRIER_PHASE) {
90                  phasesGF.add(cod);
91              }
92          }
93  
94          // Loop on Geometry-free phase measurements
95          for (CombinedObservationData cod : phasesGF) {
96              final String nameSat = setName(prn, observation.getSatelliteSystem());
97              // Check for cycle-slip detection
98              final Frequency frequency = cod.getUsedObservationData().get(0).getObservationType().getFrequency(system);
99              final boolean slip = cycleSlipDetection(nameSat, date, cod.getValue(), frequency);
100             if (!slip) {
101                 // Update cycle slip data
102                 cycleSlipDataSet(nameSat, date, cod.getValue(), cod.getUsedObservationData().get(0).getObservationType().getFrequency(system));
103             }
104         }
105 
106     }
107 
108     /**
109      * Compute if there is a cycle slip at an specific date.
110      * @param nameSat name of the satellite, on the pre-defined format (e.g.: GPS - 07 for satellite 7 of GPS constellation)
111      * @param currentDate the date at which we check if a cycle-slip occurs
112      * @param valueGF geometry free measurement
113      * @param frequency frequency used
114      * @return true if a cycle slip has been detected.
115      */
116     private boolean cycleSlipDetection(final String nameSat, final AbsoluteDate currentDate,
117                                        final double valueGF, final Frequency frequency) {
118 
119         // Access the cycle slip results to know if a cycle-slip already occurred
120         final List<CycleSlipDetectorResults>         data  = getResults();
121         final List<Map<Frequency, DataForDetection>> stuff = getStuffReference();
122 
123         // If a cycle-slip already occurred
124         if (data != null) {
125 
126             // Loop on cycle-slip results
127             for (CycleSlipDetectorResults resultGF : data) {
128 
129                 // Found the right cycle data
130                 if (resultGF.getSatelliteName().compareTo(nameSat) == 0 && resultGF.getCycleSlipMap().containsKey(frequency)) {
131                     final Map<Frequency, DataForDetection> values = stuff.get(data.indexOf(resultGF));
132                     final DataForDetection dataForDetection = values.get(frequency);
133 
134                     // Check the time gap condition
135                     final double deltaT = FastMath.abs(currentDate.durationFrom(dataForDetection.getFiguresReference()[dataForDetection.getWrite()].getDate()));
136                     if (deltaT > getMaxTimeBeetween2Measurement()) {
137                         resultGF.addCycleSlipDate(frequency, currentDate);
138                         dataForDetection.resetFigures(new SlipComputationData[getMinMeasurementNumber()], valueGF, currentDate);
139                         resultGF.setDate(frequency, currentDate);
140                         return true;
141                     }
142 
143                     // Compute the fitting polynomial if there are enough measurement since last cycle-slip
144                     if (dataForDetection.getCanBeComputed() >= getMinMeasurementNumber()) {
145                         final List<WeightedObservedPoint> xy = new ArrayList<>();
146                         for (int i = 0; i < getMinMeasurementNumber(); i++) {
147                             final SlipComputationData current = dataForDetection.getFiguresReference()[i];
148                             xy.add(new WeightedObservedPoint(1.0, current.getDate().durationFrom(currentDate),
149                                                              current.getValue()));
150                         }
151 
152                         final PolynomialCurveFitter fitting = PolynomialCurveFitter.create(2);
153                         // Check if there is a cycle_slip
154                         if (FastMath.abs(fitting.fit(xy)[0] - valueGF) > threshold) {
155                             resultGF.addCycleSlipDate(frequency, currentDate);
156                             dataForDetection.resetFigures(new SlipComputationData[getMinMeasurementNumber()], valueGF, currentDate);
157                             resultGF.setDate(frequency, currentDate);
158                             return true;
159                         }
160 
161                     } else {
162                         break;
163                     }
164 
165                 }
166 
167             }
168 
169         }
170 
171         // No cycle-slip
172         return false;
173     }
174 
175 }