1 //==============================================================================
2 //
3 // This file is part of GPSTk, the GPS Toolkit.
4 //
5 // The GPSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GPSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GPSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2020, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38
39 /**
40 * @file DoubleDifference.cpp
41 * Form double differences and buffer them, for program DDBase.
42 */
43
44 //------------------------------------------------------------------------------------
45 // TD DoubleDifference.cpp make small limit on DD buff size an input parameter
46 // TD DoubleDifference.cpp do we allow 'gaps' in ref sat's data?
47
48 //------------------------------------------------------------------------------------
49 // system includes
50 #include "TimeString.hpp"
51 // GPSTk
52
53 // DDBase
54 #include "DDBase.hpp"
55
56 //------------------------------------------------------------------------------------
57 using namespace std;
58 using namespace gpstk;
59 using namespace StringUtils;
60
61 //------------------------------------------------------------------------------------
62 // prototypes -- this module only
63 void ComputeSingleDifferences(string baseline, map<SDid,RawData>& SDmap);
64 int ComputeDoubleDifferences(map<SDid,RawData>& SDmap);
65
66 //------------------------------------------------------------------------------------
67 // other prototypes
68 // ElevationMask.cpp
69 bool ElevationMask(double elevation, double azimuth);
70
71 //------------------------------------------------------------------------------------
DoubleDifference(void)72 int DoubleDifference(void)
73 {
74 try {
75 int j,k;
76 size_t i,n;
77 // map to hold all buffered single differences for one baseline
78 map<SDid,RawData> SDmap;
79
80 if(CI.Verbose) oflog << "BEGIN DoubleDifference()"
81 << " at total time " << fixed << setprecision(3)
82 << double(clock()-totaltime)/double(CLOCKS_PER_SEC) << " seconds."
83 << endl;
84
85 // clear any existing DDs
86 DDDataMap.clear();
87
88 // loop over baselines
89 for(n=0; n<Baselines.size(); n++) {
90
91 // ----------------------------------------------------------
92 // for this baseline, compute all SDs, then DDs, and buffer them
93 if(CI.Verbose) oflog << "DoubleDifference() for baseline "
94 << Baselines[n] << endl;
95
96 // clear the SD map
97 SDmap.clear();
98
99 // ----------------------------------------------------------
100 // compute all single differences for this baseline
101 // give it same ordering as Baseline
102 ComputeSingleDifferences(Baselines[n],SDmap);
103
104 // loop over SD data, edit small ones and dump summary
105 if(CI.Verbose) oflog << "Single difference summary for baseline "
106 << Baselines[n] << endl;
107
108 vector<SDid> Remove; // these will be small dataset to delete later
109
110 map<SDid,RawData>::const_iterator kt;
111 for(k=1,kt=SDmap.begin(); kt != SDmap.end(); k++,kt++) {
112
113 if(CI.Verbose) {
114 oflog << " " << setw(2) << k << " " << kt->first
115 << " " << setw(5) << kt->second.count.size();
116 if(kt->second.count.size() > 0)
117 oflog << " " << setw(5) << kt->second.count.at(0) << " - "
118 << setw(5) << kt->second.count.at(kt->second.count.size()-1);
119 else
120 oflog << " na - na";
121
122 // gaps - (count : number of pts)
123 if(kt->second.count.size() > 0) { // gcc needs this ...
124 for(i=0; i<kt->second.count.size()-1; i++) {
125 j = kt->second.count.at(i+1) - kt->second.count.at(i);
126 if(j > 1) oflog
127 << " (" << kt->second.count.at(i)+1 << ":" << j-1 << ")";
128 }
129 }
130 }
131
132 // ignore small datasets
133 if(kt->second.count.size() < 10) { // TD make input parameter
134 Remove.push_back(kt->first);
135 if(CI.Verbose) oflog << " **Rejected";
136 }
137
138 if(CI.Verbose) oflog << endl;
139
140 } // end summary loop
141
142 // delete marked SD buffers
143 for(i=0; i<Remove.size(); i++) SDmap.erase(Remove[i]);
144
145 // ----------------------------------------------------------
146 // now compute double differences - according to timetable
147 if(ComputeDoubleDifferences(SDmap)) return 1;
148
149 // check that there are non-zero double differences
150
151 } // end loop over baselines
152
153 return 0;
154 }
155 catch(Exception& e) { GPSTK_RETHROW(e); }
156 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
157 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
158 } // end DoubleDifference()
159
160 //------------------------------------------------------------------------------------
161 // Compute all single differences 'site1' - 'site2', using the RawDataBuffers in
162 // Stations[site], and store the results in the given map<SDid,RawData>.
ComputeSingleDifferences(string baseline,map<SDid,RawData> & SDmap)163 void ComputeSingleDifferences(string baseline, map<SDid,RawData>& SDmap)
164 {
165 try {
166 int beg,end;
167 size_t i,j;
168 GSatID sat;
169
170 // decompose the baseline
171 string site1=word(baseline,0,'-');
172 string site2=word(baseline,1,'-');
173
174 // find the beginning and ending *counts* of good data for this baseline
175 if(QueryTimeTable(baseline,beg,end)) {
176 oflog << "ERROR - baseline " << baseline
177 << " not found in timetable. No single differences computed." << endl;
178 return;
179 }
180
181 // find satellites in common
182 map<GSatID,RawData>::const_iterator it1,it2;
183
184 // loop over satellites at first site
185 for(it1 = Stations[site1].RawDataBuffers.begin();
186 it1 != Stations[site1].RawDataBuffers.end(); it1++) {
187
188 sat = it1->first;
189 // it1->second is RawData={ L1,L2,P1,P2,elev,az,count buffers = vector<> }
190
191 // does this sat have data at the other station?
192 it2 = Stations[site2].RawDataBuffers.find(sat);
193 if(it2 == Stations[site2].RawDataBuffers.end()) continue; // no
194
195 // compute single differences for this satellite
196 // here is where you define the ordering of sites: first(1) - second(2)
197 SDid sdid(site1,site2,sat);
198 RawData sddata;
199
200 // loop over epochs, finding common data. start and stop the loop
201 // at times determined by the timetable, NOT by the raw data buffers.
202 i = j = 0;
203 while(i < it1->second.count.size() && j < it2->second.count.size()) {
204
205 // impose limits from timetable
206 if(it1->second.count[i] > end) break;
207 else if(it2->second.count[j] > end) break;
208 else if(it1->second.count[i] < beg) i++;
209 else if(it2->second.count[j] < beg) j++;
210 // i and j are the same count (epoch)
211 else if(it1->second.count[i] == it2->second.count[j]) {
212 // reject data below MinElevation here
213 //if(it1->second.elev[i] > CI.MinElevation &&
214 //it2->second.elev[j] > CI.MinElevation) {
215 if(ElevationMask(it1->second.elev[i],it1->second.az[i]) &&
216 ElevationMask(it2->second.elev[j],it2->second.az[j])) {
217
218 // buffer the differences
219 sddata.count.push_back(it1->second.count[i]);
220 sddata.L1.push_back(it1->second.L1[i] - it2->second.L1[j]);
221 sddata.L2.push_back(it1->second.L2[i] - it2->second.L2[j]);
222 sddata.P1.push_back(it1->second.P1[i] - it2->second.P1[j]);
223 sddata.P2.push_back(it1->second.P2[i] - it2->second.P2[j]);
224 sddata.ER.push_back(it1->second.ER[i] - it2->second.ER[j]);
225 sddata.elev.push_back(it1->second.elev[i]);
226
227 } // end if elevation is ok
228
229 // next epoch
230 i++;
231 j++;
232 }
233 // i is behind j in time(count)
234 else if(it1->second.count[i] < it2->second.count[j])
235 i++;
236 // i is ahead of j in time(count)
237 else
238 j++;
239
240 } // end while
241
242 // save it in the map
243 SDmap[sdid] = sddata;
244
245 } // end loop over satellites at first site
246
247 }
248 catch(Exception& e) { GPSTK_RETHROW(e); }
249 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
250 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
251 }
252
253 //------------------------------------------------------------------------------------
254 // Assume SDmap is all for the same baseline
ComputeDoubleDifferences(map<SDid,RawData> & SDmap)255 int ComputeDoubleDifferences(map<SDid,RawData>& SDmap)
256 {
257 try {
258 bool frst,ok;
259 int indx,count = 0,ddsign;
260 long nn1,nn2;
261 double ddL1,ddL2,ddER,ddP1,ddP2,dd,db1,db2;
262 CommonTime tt,ttnext; // ttnext is the time of the next reference satellite switch
263 //SDid sid,ref; // SDid of the current satellite and reference satellite
264 map<SDid,int> Inext; // index in count (all) buffers which is to be processed next
265 map<SDid,RawData>::const_iterator it;
266
267 if(SDmap.size() == 0) return 0;
268
269 // initialize the 'next index' map to zero
270 // find the smallest (earliest) count
271 frst = true;
272 for(it=SDmap.begin(); it != SDmap.end(); it++) {
273 int jj=0;
274 Inext[it->first] = jj;
275 if(frst || it->second.count[0] < count) {
276 count = it->second.count[0];
277 frst = false;
278 }
279 }
280
281 // ref will be the SDid for the reference satellite
282 SDid ref = SDmap.begin()->first; // ref.sat is TBD by timetable
283
284 // loop over epochs in the SDs
285 ttnext = CommonTime::BEGINNING_OF_TIME;
286 while(1) {
287 // time at this count
288 tt = FirstEpoch + count * CI.DataInterval;
289
290 // get the reference satellite at this time
291 if(tt > ttnext) {
292 ttnext = tt;
293 if(QueryTimeTable(ref, ttnext)) { // error - timetable failed
294 oflog << "DD: Error - failed to find reference from timetable at "
295 << printTime(tt,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << " count "
296 << count << " for baseline " << ref.site1 << "-" << ref.site2 << endl;
297 return 1;
298 }
299 if(CI.Verbose) oflog << "DD: reference is set to " << ref << " at "
300 << printTime(tt,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g")
301 << " count " << count << endl;
302 }
303
304 // does reference satellite have data at this count?
305 if(SDmap[ref].count[Inext[ref]] != count) {
306 oflog << "Error - failed to find reference data " << ref << " at "
307 << printTime(tt,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << endl;
308 // TD return here, or just skip the epoch?
309 // question is do we allow 'holes' in ref sat's data?
310 return 1;
311 }
312
313 // compute DDs
314 for(it=SDmap.begin(); it != SDmap.end(); it++) {
315
316 // sid is the SDid for the current satellite
317 SDid sid = it->first;
318 indx = Inext[sid];
319
320 // end of buffer has been reached
321 if(Inext[sid] >= int(SDmap[sid].count.size())) continue;
322
323 if(sid == ref) continue; // ignore the reference
324
325 // no data for this satellite at this count
326 if(SDmap[sid].count[indx] != count) continue;
327
328 // compute DD phases and DD nominal range
329 ddL1 = wl1 * (SDmap[sid].L1[indx] - SDmap[ref].L1[Inext[ref]]);
330 ddL2 = wl2 * (SDmap[sid].L2[indx] - SDmap[ref].L2[Inext[ref]]);
331 ddP1 = SDmap[sid].P1[indx] - SDmap[ref].P1[Inext[ref]];
332 ddP2 = SDmap[sid].P2[indx] - SDmap[ref].P2[Inext[ref]];
333 ddER = (SDmap[sid].ER[indx] - SDmap[ref].ER[Inext[ref]]);
334
335 // get the appropriate DDData from the map, or create a new one
336 map<DDid,DDData>::iterator jt;
337 DDid ddid((ref.ssite == 1 ? ref.site1 : ref.site2),
338 (ref.ssite == 1 ? ref.site2 : ref.site1),sid.sat,ref.sat);
339 if(DDDataMap.find(ddid) == DDDataMap.end()) {
340 // create a new DDData
341 DDData tddb;
342 dd = (-ddL1+ddER)/wl1;
343 nn1 = int(dd + (dd > 0 ? 0.5 : -0.5));
344 tddb.L1bias = wl1 * nn1;
345 dd = (-ddL2+ddER)/wl2;
346 nn2 = int(dd + (dd > 0 ? 0.5 : -0.5));
347 tddb.L2bias = wl2 * nn2;
348 oflog << " Phase bias (initial) on " << ddid
349 << " at " << setw(4) << count << " "
350 << printTime(tt,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g");
351 if(CI.Frequency != 2) oflog << " L1: " << setw(10) << nn1;
352 if(CI.Frequency != 1) oflog << " L2: " << setw(10) << nn2;
353 oflog << endl;
354 //tddb.lastresetcount = count;
355 tddb.resets.push_back(tddb.count.size()); // always one at beginning
356 tddb.prevL1 = (ddL1-ddER)+tddb.L1bias;
357 tddb.prevL2 = (ddL2-ddER)+tddb.L2bias;
358 DDDataMap[ddid] = tddb;
359 }
360
361 // get the current DDData structure, and relative sign
362 jt = DDDataMap.find(ddid); // never fail...
363 ddsign = DDid::compare(ddid,jt->first);
364 DDData& ddb=jt->second;
365 ok = true; // if ok, buffer this DDData = ddb
366
367 // reset bias?
368 db1 = (ddsign*(ddL1-ddER) + ddb.L1bias - ddb.prevL1)/wl1;
369 db2 = (ddsign*(ddL2-ddER) + ddb.L2bias - ddb.prevL2)/wl2;
370 if((CI.Frequency != 2 && fabs(db1) > CI.PhaseBiasReset) ||
371 (CI.Frequency != 1 && fabs(db2) > CI.PhaseBiasReset)) {
372 long ndb1 = long(db1 + (db1 > 0 ? 0.5 : -0.5));
373 long ndb2 = long(db2 + (db2 > 0 ? 0.5 : -0.5));
374 oflog << " Phase bias (reset ) on " << ddid
375 << " at " << setw(4) << count << " "
376 << printTime(tt,"%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g");
377 if(CI.Frequency != 2) oflog << " L1: " << setw(10) << ndb1;
378 if(CI.Frequency != 1) oflog << " L2: " << setw(10) << ndb2;
379 oflog << endl;
380 ddb.L1bias -= wl1 * ndb1;
381 ddb.L2bias -= wl2 * ndb2;
382 //ddb.lastresetcount = count;
383 ddb.resets.push_back(ddb.count.size());
384 }
385
386 // remove the bias from the data
387 ddL1 += ddsign * ddb.L1bias;
388 ddL2 += ddsign * ddb.L2bias;
389
390 // save for next time
391 ddb.prevL1 = ddsign*(ddL1-ddER);
392 ddb.prevL2 = ddsign*(ddL2-ddER);
393
394 // buffer the debiased DDs
395 ddb.DDL1.push_back(ddsign*ddL1);
396 ddb.DDL2.push_back(ddsign*ddL2);
397 ddb.DDP1.push_back(ddsign*ddP1);
398 ddb.DDP2.push_back(ddsign*ddP2);
399 ddb.DDER.push_back(ddsign*ddER);
400 ddb.count.push_back(count);
401
402 // increment Inext
403 Inext[sid]++;
404
405 } // end loop over SD's in SDmap
406
407 Inext[ref]++;
408
409 // find the next count
410 // quit when all Inext are at end
411 frst = true;
412 ok = false;
413 for(it=SDmap.begin(); it != SDmap.end(); it++) {
414 if(Inext[it->first] < int(SDmap[it->first].count.size())) {
415 if(frst || it->second.count[Inext[it->first]] < count) {
416 count = it->second.count[Inext[it->first]];
417 frst = false;
418 }
419 ok = true;
420 }
421 }
422
423 if(!ok) break;
424
425 } // end while loop over epochs
426
427 return 0;
428 }
429 catch(Exception& e) { GPSTK_RETHROW(e); }
430 catch(std::exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
431 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
432 }
433
434 //------------------------------------------------------------------------------------
435 //------------------------------------------------------------------------------------
436