1   /* Copyright 2002-2022 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.io.File;
20  import java.io.FileInputStream;
21  import java.io.IOException;
22  import java.net.URISyntaxException;
23  import java.util.Arrays;
24  import java.util.List;
25  
26  import org.junit.Assert;
27  import org.junit.Before;
28  import org.junit.Test;
29  import org.orekit.Utils;
30  import org.orekit.data.DataFilter;
31  import org.orekit.data.GzipFilter;
32  import org.orekit.data.DataSource;
33  import org.orekit.data.UnixCompressFilter;
34  import org.orekit.gnss.Frequency;
35  import org.orekit.gnss.HatanakaCompressFilter;
36  import org.orekit.gnss.ObservationDataSet;
37  import org.orekit.gnss.RinexObservationLoader;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.time.TimeScalesFactory;
40  
41  
42  public class PhaseMinusCodeCycleSlipDetectorTest {
43      
44      @Before
45      public void setUp() {
46          Utils.setDataRoot("regular-data");
47      }
48  
49      @Test
50      public void testTheNumberOFCycleFind() throws URISyntaxException, IOException {
51          
52          final String inputPath = GeometryFreeCycleSlipDetectorTest.class.getClassLoader().getResource("gnss/cycleSlip/seat0440.16d.Z").toURI().getPath();
53          final File input  = new File(inputPath);
54          String fileName = "seat0440.16d.Z";
55          DataSource nd = new DataSource(fileName,
56                                       () -> new FileInputStream(new File(input.getParentFile(), fileName)));
57          for (final DataFilter filter : Arrays.asList(new GzipFilter(),
58                                                       new UnixCompressFilter(),
59                                                       new HatanakaCompressFilter())) {
60              nd = filter.filter(nd);
61          }
62          final RinexObservationLoader loader = new RinexObservationLoader(nd);
63          final List<ObservationDataSet> obserDataSets = loader.getObservationDataSets();
64          PhaseMinusCodeCycleSlipDetector slipDetectors =
65              new PhaseMinusCodeCycleSlipDetector(90, 10, 20, 3); 
66          final List<CycleSlipDetectorResults> results = slipDetectors.detect(obserDataSets);
67          for(CycleSlipDetectorResults d: results) {
68              switch(getPrn(d)) {
69                  case 1: 
70                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,  2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
71                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,33 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
72  
73                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,  2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
74                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,40 ,00.0000000, TimeScalesFactory.getTAI())),1e-9);
75                      break;
76                      
77                  case 5: 
78                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,  2, 13,  2, 44, 30.0000000, TimeScalesFactory.getTAI())),1e-9);
79                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
80  
81                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,  2, 13,  2, 44, 30.0000000, TimeScalesFactory.getTAI())),1e-9);
82                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
83                      break;
84                      
85                  case 6: 
86                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016, 2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
87                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,  2, 13,  4, 28, 30.0000000, TimeScalesFactory.getTAI())),1e-9);
88  
89                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016, 2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
90                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,  2, 13,  4, 31, 0.0000000, TimeScalesFactory.getTAI())),1e-9);
91                      break;
92  
93                  case 7: 
94                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,  2, 13,  4, 13,  30.0000000, TimeScalesFactory.getTAI())),1e-9);
95                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
96  
97                      Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,  2, 13,  4, 11,  00.0000000, TimeScalesFactory.getTAI())),1e-9);
98                      Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
99                      break;
100 
101                 case 9: 
102                     Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016, 2, 13,  2, 32,  00.0000000, TimeScalesFactory.getTAI())),1e-9);
103                     Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
104 
105                     Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016, 2, 13,  2, 31,  30.0000000, TimeScalesFactory.getTAI())),1e-9);
106                     Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
107                     break;
108 
109                 case 11: 
110                     Assert.assertEquals(19.0, d.getEndDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,  2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
111                     Assert.assertEquals(19.0, d.getBeginDate(Frequency.G01).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
112 
113                     Assert.assertEquals(19.0, d.getEndDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,  2, 13,  5,  0,  0.0000000, TimeScalesFactory.getTAI())),1e-9);
114                     Assert.assertEquals(19.0, d.getBeginDate(Frequency.G02).durationFrom(new AbsoluteDate(2016,2, 13  ,2  ,31 ,30.0000000, TimeScalesFactory.getTAI())),1e-9);
115                     break;
116 
117                 default:   break;
118             }   
119         }  
120     }
121 
122     //Test to verify that the cycle-slips because of data gap are computed
123     @Test
124     public void testTimeCycleSlip() throws URISyntaxException, IOException {
125         
126         final String inputPath = GeometryFreeCycleSlipDetectorTest.class.getClassLoader().getResource("gnss/cycleSlip/WithCycleSlip.16o").toURI().getPath();
127         final File input  = new File(inputPath);
128         String fileName = "WithCycleSlip.16o";
129         DataSource nd = new DataSource(fileName,
130                                      () -> new FileInputStream(new File(input.getParentFile(), fileName)));
131         for (final DataFilter filter : Arrays.asList(new GzipFilter(),
132                                                      new UnixCompressFilter(),
133                                                      new HatanakaCompressFilter())) {
134             nd = filter.filter(nd);
135         }
136         final RinexObservationLoader loader = new RinexObservationLoader(nd);
137         final List<ObservationDataSet> obserDataSets = loader.getObservationDataSets();
138         PhaseMinusCodeCycleSlipDetector slipDetectors =
139             new PhaseMinusCodeCycleSlipDetector(90, 1e15, 20, 3);
140         final List<CycleSlipDetectorResults> results = slipDetectors.detect(obserDataSets);
141         for(CycleSlipDetectorResults d: results) {
142             switch(getPrn(d)) {
143                 case 1: 
144                     //The date have been created  manually within the file
145                     AbsoluteDate[] dateCycleSlipL1 = new AbsoluteDate[] {
146                         new AbsoluteDate(2016,02,13,04,37,43.000 , TimeScalesFactory.getUTC()),
147                         new AbsoluteDate(2016,02,13,04,45,13.000 , TimeScalesFactory.getUTC()),
148                         new AbsoluteDate(2016,02,13,04,54,13.000 , TimeScalesFactory.getUTC())};
149                     int i1 = 0;
150                     for (AbsoluteDate dateL1: d.getCycleSlipMap().get(Frequency.G01)) {
151                         Assert.assertEquals(0,  dateL1.compareTo(dateCycleSlipL1[i1]));
152                         i1++;
153                     }
154                     //The dates have been created manually within the file
155                     AbsoluteDate[] dateCycleSlipL2 = new AbsoluteDate[] {
156                         new AbsoluteDate(2016,02,13,04,38,13.000 , TimeScalesFactory.getUTC()),
157                         new AbsoluteDate(2016,02,13,04,41,13.000 , TimeScalesFactory.getUTC()),
158                         new AbsoluteDate(2016,02,13,04,45,43.000 , TimeScalesFactory.getUTC()),
159                         new AbsoluteDate(2016,02,13,04,54,13.000 , TimeScalesFactory.getUTC())};
160                     int i2 = 0;
161                     for(AbsoluteDate dateL2: d.getCycleSlipMap().get(Frequency.G02)) {
162                         Assert.assertEquals(0,  dateL2.compareTo(dateCycleSlipL2[i2]));
163                         i2++;
164                     }
165                 default:    break;
166                                 
167                    
168             }
169         }
170     }
171 
172     //Test to check the detectors find the cycle slip added to data on purpose
173     @Test
174     public void testCycleSlipDetection() throws URISyntaxException, IOException {
175         final String inputPath = GeometryFreeCycleSlipDetectorTest.class.getClassLoader().getResource("gnss/cycleSlip/WithCycleSlip.16o").toURI().getPath();
176         final File input  = new File(inputPath);
177         String fileName = "WithoutCycleSlip.16o";
178         DataSource nd = new DataSource(fileName,
179                                      () -> new FileInputStream(new File(input.getParentFile(), fileName)));
180         for (final DataFilter filter : Arrays.asList(new GzipFilter(),
181                                                      new UnixCompressFilter(),
182                                                      new HatanakaCompressFilter())) {
183             nd = filter.filter(nd);
184         }
185         final RinexObservationLoader loader = new RinexObservationLoader(nd);
186         final List<ObservationDataSet> obserDataSets = loader.getObservationDataSets();
187         final double dt = 31; //great time gap threshold to don't detect cycle-slip because of time gap
188         final int N = 25;
189         final int m = 2;
190         //Test on L2
191         final double[] thresholdL2 = new double[] {
192             0.4,2};
193         //The date have been computed with an excel spreadsheet.
194         final AbsoluteDate[] dateL2 = new AbsoluteDate[] {
195           new AbsoluteDate(2016, 2, 13 ,1, 43, 13.000, TimeScalesFactory.getUTC()),
196         };
197         
198         for(int i = 0; i<thresholdL2.length; i++) {
199             PhaseMinusCodeCycleSlipDetector slipDetectorsL2 =
200                             new PhaseMinusCodeCycleSlipDetector(dt, thresholdL2[i], N, m);
201             final List<CycleSlipDetectorResults> resultsL2 = slipDetectorsL2.detect(obserDataSets);
202             for(CycleSlipDetectorResults d : resultsL2) {
203                 if(i == 0) {
204                     final List<AbsoluteDate> computedDateOnL2 = d.getCycleSlipMap().get(Frequency.G02);
205                     Assert.assertEquals(3, computedDateOnL2.size());
206                     Assert.assertEquals(0.0, computedDateOnL2.get(0).durationFrom(dateL2[i]), 0.0);
207                 } else {
208                     final List<AbsoluteDate> computedDateOnL2 = d.getCycleSlipMap().get(Frequency.G02);
209                     Assert.assertEquals(0, computedDateOnL2.size());
210                 }
211 
212             }
213         }
214         /////////////////////////////////////////////////////////////////////////////////////
215         //Test on L1
216         final double[] thresholdL1 = new double[] {
217             0.1,0.30,2};
218         //The date have been computed with an excel spreadsheet.
219         final AbsoluteDate[] dateL1 = new AbsoluteDate[] {
220             new AbsoluteDate(2016, 2, 13 ,1, 43, 13.000, TimeScalesFactory.getUTC()),
221             new AbsoluteDate(2016, 02, 13, 01, 55, 43.000, TimeScalesFactory.getUTC()),
222             new AbsoluteDate(2016, 02, 13, 02, 8, 13.000, TimeScalesFactory.getUTC())
223           };
224         for (int i=0; i<thresholdL1.length; i++) {
225             PhaseMinusCodeCycleSlipDetector slipDetectorsL1 =
226                             new PhaseMinusCodeCycleSlipDetector(dt, thresholdL1[i], N, m+1);
227             final List<CycleSlipDetectorResults> resultsL1 = slipDetectorsL1.detect(obserDataSets);
228             for(CycleSlipDetectorResults d : resultsL1) {
229                 if(i == 0) {
230                     final List<AbsoluteDate> computedDateOnL1 = d.getCycleSlipMap().get(Frequency.G01);
231                     Assert.assertEquals(3, computedDateOnL1.size());
232                     int i1 = 0;
233                     for(AbsoluteDate date: computedDateOnL1) {
234                         Assert.assertEquals(0.0, date.durationFrom(dateL1[i1]), 1e-9);
235                         i1++;
236                     }
237                 } else if (i == 1) {
238                     final List<AbsoluteDate> computedDateOnL1 = d.getCycleSlipMap().get(Frequency.G01);
239                     Assert.assertEquals(2, computedDateOnL1.size());
240                 } else {
241                     final List<AbsoluteDate> computedDateOnL1 = d.getCycleSlipMap().get(Frequency.G01);
242                     Assert.assertEquals(0, computedDateOnL1.size());
243                 }
244             }
245         }
246     }
247     
248     /** Getter on the PRN of the satellite. */
249     private int getPrn(final CycleSlipDetectorResults d) {
250         
251         if(d.getSatelliteName().substring(6).compareTo("1")==0) {return 1;};
252         if(d.getSatelliteName().substring(6).compareTo("2")==0) {return 2;};
253         if(d.getSatelliteName().substring(6).compareTo("3")==0) {return 3;};
254         if(d.getSatelliteName().substring(6).compareTo("4")==0) {return 4;};
255         if(d.getSatelliteName().substring(6).compareTo("5")==0) {return 5;};
256         if(d.getSatelliteName().substring(6).compareTo("6")==0) {return 6;};
257         if(d.getSatelliteName().substring(6).compareTo("7")==0) {return 7;};
258         if(d.getSatelliteName().substring(6).compareTo("8")==0) {return 8;};
259         if(d.getSatelliteName().substring(6).compareTo("9")==0) {return 9;};
260         if(d.getSatelliteName().substring(6).compareTo("10")==0) {return 10;};
261         if(d.getSatelliteName().substring(6).compareTo("11")==0) {return 11;};
262         if(d.getSatelliteName().substring(6).compareTo("12")==0) {return 12;};
263         if(d.getSatelliteName().substring(6).compareTo("13")==0) {return 13;};
264         if(d.getSatelliteName().substring(6).compareTo("14")==0) {return 14;};
265         if(d.getSatelliteName().substring(6).compareTo("15")==0) {return 15;};
266         if(d.getSatelliteName().substring(6).compareTo("16")==0) {return 16;};
267         if(d.getSatelliteName().substring(6).compareTo("17")==0) {return 17;};
268         if(d.getSatelliteName().substring(6).compareTo("18")==0) {return 18;};
269         if(d.getSatelliteName().substring(6).compareTo("19")==0) {return 19;};
270         if(d.getSatelliteName().substring(6).compareTo("20")==0) {return 20;};
271         if(d.getSatelliteName().substring(6).compareTo("21")==0) {return 21;};
272         if(d.getSatelliteName().substring(6).compareTo("22")==0) {return 22;};
273         if(d.getSatelliteName().substring(6).compareTo("23")==0) {return 23;};
274         if(d.getSatelliteName().substring(6).compareTo("24")==0) {return 24;};
275         if(d.getSatelliteName().substring(6).compareTo("25")==0) {return 25;};
276         if(d.getSatelliteName().substring(6).compareTo("26")==0) {return 26;};
277         if(d.getSatelliteName().substring(6).compareTo("27")==0) {return 27;};
278         if(d.getSatelliteName().substring(6).compareTo("28")==0) {return 28;};
279         if(d.getSatelliteName().substring(6).compareTo("29")==0) {return 29;};
280         if(d.getSatelliteName().substring(6).compareTo("30")==0) {return 30;};
281         if(d.getSatelliteName().substring(6).compareTo("31")==0) {return 31;} else {return 32;}
282               
283     }
284 }