gdc.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk 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 GNSSTk 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 GNSSTk; 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-2022, 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 
47 #include "gdc.hpp"
48 #include "GNSSconstants.hpp"
49 #include "logstream.hpp"
50 #include "stl_helpers.hpp"
51 
52 using namespace std;
53 
54 namespace gnsstk
55 {
56 
57  //---------------------------------------------------------------------------------
58  // static data
59  //---------------------------------------------------------------------------------
60  // class gdc: string giving version of gnsstk Discontinuity Corrector
61  const string gdc::GDCVersion = string("9.0 5/20/17");
62 
63  /* ----------------------- flags and bitmaps
64  values for flags[]; NB flags[] is either good (0) or bad (non-zero)
65  not to be confused with Arc::marks or SatPass flags TD bitmap? */
66  const unsigned gdc::OK = 0; // good; NB SatPass::OK == 1
67  const unsigned gdc::BAD = 1; // bad in SatPass; NB SatPass::BAD == 0
68  const unsigned gdc::WLOUTLIER = 2; // called outlier by WL filter
69  const unsigned gdc::GFOUTLIER = 3; // called outlier by GF filter
70  const unsigned gdc::WLSHORT = 4; // data with Arc.ngood < MinPts
71  const unsigned gdc::GFSHORT = 5; // data with Arc.ngood < MinPts
72  const unsigned gdc::ISOLATED =
73  6; // final check - isolated good points (<MinPts)
74 
75  // Bitmap values used by Arc::mark - see create_mark_string_map() routine.
76  const unsigned Arc::BEG = 1;
77  const unsigned Arc::WLSLIP = 2;
78  const unsigned Arc::GFSLIP = 4;
79  const unsigned Arc::WLFIX = 8;
80  const unsigned Arc::GFFIX = 16;
81  const unsigned Arc::WLMARK = 32;
82  const unsigned Arc::GFMARK = 64;
83  const unsigned Arc::REJ = 128;
84  // NB add additional values to the Arc::create_mark_string_map() routine.
85 
86  // ------------------------ mere conveniences
87  const unsigned gdc::WL = 0; // keep these 0,1 for use as array indexes
88  const unsigned gdc::GF = 1;
89  const vector<unsigned> gdc::SLIP = gdc::create_vector_SLIP();
90  const vector<unsigned> gdc::FIX = gdc::create_vector_FIX();
91  const vector<string> gdc::LAB = gdc::create_vector_LAB();
92 
93  // static const map giving string descriptors for each Arc::mark value
94  const map<unsigned, string> Arc::markStr = Arc::create_mark_string_map();
95 
96  // used to access CFG - also see setcfg(a,b,c) below
97  #define cfg(a) cfg_func(#a) // #a 'string-izes' a "a"
98 
99  //---------------------------------------------------------------------------------
100  //---------------------------------------------------------------------------------
101  /* See doc in .hpp
102  what it does
103  Initialize() define globals and wavelengths, check GLOn and compute if
104  nec. FillDataVectors define first Arc(BEG), fill 5 parallel arrays
105  xdata=time-begin, (dataWL, dataGF, dataIF), flag=OK or BAD
106  units cycles(WL,GF,NL), remove biases WLbias, GFbias, IFbias,
107  get N12bias internally data (WL,GF,IF,GR) has units
108  wl(WL,GF,NL,(meters))
109  but dump all in meters.
110  ProcessOneCombo(WL)
111  GrossProcessing
112  filterFirstDiff - filter using first diff
113  mergeFilterResultsIntoArcs
114  outlier or short seg -> set flag for all segment;
115  slip: if already a slip of this type, "unfix";
116  else add an Arc and mark "slip" ** define slip, not
117  Nslip **
118  getArcStats
119  findLargeGaps
120  fixSlips fix from here out using ** step, ADD to Nslip **
121  FineProcessing
122  filterWindow
123  mergeFilterResultsIntoArcs
124  getArcStats
125  fixSlips
126  end ProcessOneCombo
127 
128  ProcessOneCombo(GF)
129 
130  final "double check" -- tbd
131 
132  applyFixesToSatPass - only routine to change SatPass - apply all fixes
133  and outputs to SatPass and editing commands.
134  NB
135  Arc "floats above" data vectors and does not modify them at all
136  Arcs don't "know" about each other - type BEG,slip,rej -- see top of
137  gdc.cpp flag vector also does not change data vectors, bitmaps -- see top
138  of gdc.cpp flag is either good (0) or not (!0); flags NOT related to
139  SatPass or Arc::marks Arcs handles BEG, SLIP/FIX, and GAP(new BEG); flags
140  handle bad data w/in Arc NB BAD != SatPass::BAD, but ONLY SatPass::BAD on
141  input produces flag = BAD
142 
143  convert SatPass to data arrays and call DC(arrays) */
144  int gdc::DiscontinuityCorrector(SatPass& SP, std::string& retMsg,
145  std::vector<std::string>& cmds, int GLOn)
146  {
147  try
148  {
149  sat = SP.getSat();
150  isGLO = (sat.system == SatelliteSystem::Glonass);
151  GLOchan = GLOn;
152  // if GLONASS frequency channel not given, try to find it
153  if (isGLO && GLOchan == -99)
154  {
155  if (!SP.getGLOchannel(GLOchan, retMsg))
156  {
157  retMsg =
158  " Error - unable to compute GLO channel - fail: " + retMsg;
159  return -9;
160  }
161  LOG(VERBOSE) << "# Compute GLO channel = " << GLOchan << " "
162  << retMsg;
163  }
164 
165  // get obstypes for this pass
166  vector<string> obstypes(SP.getObsTypes());
167  string L1("L1"), L2("L2"), P1("P1"), P2("P2");
168  // useCA? no, assume caller knows what he is doing and only gave you C
169  // || P
170  if (vectorindex(obstypes, string("P1")) == -1)
171  {
172  P1 = string("C1");
173  }
174  if (vectorindex(obstypes, string("P2")) == -1)
175  {
176  P2 = string("C2");
177  }
178 
179  outfmt = SP.getOutputFormat();
180  Epoch beg = SP.getFirstTime();
181 
182  vector<double> L1_in, L2_in, P1_in, P2_in, dt_in;
183  vector<int> flags_in;
184 
185  // loop over the pass - MUST keep flags_in, dt_in, arrays all parallel
186  for (unsigned int i = 0; i < SP.size(); i++)
187  {
188 
189  // save the seconds since beg
190  dt_in.push_back(SP.time(i) - beg);
191 
192  // test for good data
193  // must consistently mark bad data in SP with SatPass::BAD
194  if (!(SP.spdvector[i].flag & SatPass::OK) ||
195  SP.data(i, L1) == 0.0 || SP.data(i, L2) == 0.0 ||
196  SP.data(i, P1) == 0.0 || SP.data(i, P2) == 0.0)
197  {
198  flags_in.push_back(0); // 0 bad - as in SatPass
199  L1_in.push_back(0.0);
200  L2_in.push_back(0.0);
201  P1_in.push_back(0.0);
202  P2_in.push_back(0.0);
203  continue;
204  }
205 
206  // good data
207  flags_in.push_back(1); // 1 good data - as in SatPass
208  L1_in.push_back(SP.data(i, L1));
209  L2_in.push_back(SP.data(i, L2));
210  P1_in.push_back(SP.data(i, P1));
211  P2_in.push_back(SP.data(i, P2));
212  }
213 
214  // save the first GDC output line from SatPass (SPS)
215  {
216  ostringstream oss;
217  oss << "GDC " << setw(3) << unique + 1 << " SPS " << SP;
218  SPSstr = oss.str();
219  }
220 
221  int iret = DiscontinuityCorrector(sat, SP.getDT(), beg, L1_in, L2_in,
222  P1_in, P2_in, dt_in, flags_in,
223  retMsg, cmds, GLOchan, outfmt);
224  if (iret)
225  {
226  return iret;
227  }
228 
229  // apply fixes to SatPass
230  if (cfg(doFix))
231  {
232  applyFixesToSatPass(SP); // ,breaks,marks);
233  }
234 
235  return 0;
236  }
237  catch (Exception& e)
238  {
239  GNSSTK_RETHROW(e);
240  }
241  } // end int gdc::DiscontinuityCorrector(SatPass& SP, string& retMsg, int GLOn)
242 
243  //---------------------------------------------------------------------------------
244  /* Call to DC without SatPass.
245  Flags on input must be either 1(OK) or 0(BAD), as in SatPass */
246  int gdc::DiscontinuityCorrector(
247  const RinexSatID& sat_in, const double& nominalDT, const Epoch& beginTime,
248  std::vector<double> dataL1, std::vector<double> dataL2,
249  std::vector<double> dataP1, std::vector<double> dataP2,
250  std::vector<double> dt_in, std::vector<int> flags_in, std::string& retMsg,
251  std::vector<std::string>& cmds, int GLOn, std::string outfmt_in)
252  {
253  try
254  {
255  int iret;
256  unsigned int i;
257 
258  // TD check that arrays are all full length
259 
260  sat = sat_in;
261  dt = nominalDT;
262  beginT = beginTime;
263  beginT += dt_in[0];
264  outfmt = outfmt_in;
265 
266  isGLO = (sat.system == SatelliteSystem::Glonass);
267  GLOchan = GLOn;
268  if (isGLO && GLOchan == -99)
269  {
270  return -9;
271  }
272 
273  // NB wl1,alpha,beta are used only in this routine...
274  wl1 = getWavelength(sat.system, 1, GLOchan); // GLOchan ignored by GPS
275  wl2 = getWavelength(sat.system, 2, GLOchan);
276  alpha = getAlpha(sat.system, 1, 2);
277  beta = getBeta(sat.system, 1, 2);
278  wlWL =
279  wl2 * (beta + 1.0) / alpha; // wl(WL) = 86cm GPS, depends on GLOchan
280  wlGF = wl2 - wl1; // wl(GF) = wl1-wl2 = 5.376cm GPS, or f(GLOchan)
281  // wl(NL) = 10.7cm GPS, used for IF
282 
283  // fill data vectors from input -------------------------------
284  Arc arc(0, 0, 0, Arc::BEG);
285  xdata.clear();
286  flags.clear();
287  dataWL.clear();
288  dataGF.clear();
289 
290  // loop over the pass - MUST keep xdata, flags, dataWL and dataGF
291  // parallel
292  double d, dtlast;
293  arc.ngood = 0;
294  for (i = 0; i < dt_in.size(); i++)
295  {
296  // save the seconds since beginT
297  xdata.push_back(dt_in[i]);
298 
299  // test for good data
300  // caller must consistently mark bad data with SatPass::BAD(0)
301  if (!(flags_in[i] & SatPass::OK) || dataL1[i] == 0.0 ||
302  dataL2[i] == 0.0 || dataP1[i] == 0.0 || dataP2[i] == 0.0)
303  {
304  flags.push_back(BAD); // bad data from SatPass
305  dataWL.push_back(0.0);
306  dataGF.push_back(0.0);
307  continue;
308  }
309 
310  // good data
311  dtlast = dt_in[i];
312  flags.push_back(OK); // 0 good data
313  arc.ngood++;
314 
315  // WLC = (WLphase - NLrange) in units of WLwl
316  d = ((beta * wl1 * dataL1[i] - wl2 * dataL2[i]) / (beta - 1.0) -
317  (beta * dataP1[i] + dataP2[i]) / (beta + 1.0)) /
318  wlWL;
319  if (arc.ngood == 1)
320  {
321  WLbias = d;
322  }
323  dataWL.push_back(d - WLbias);
324 
325  // LGF = wl1*L1 - wl2*L2 in units of GFwl
326  d = (wl1 * dataL1[i] - wl2 * dataL2[i]) / wlGF;
327  if (arc.ngood == 1)
328  {
329  GFbias = d;
330  }
331  dataGF.push_back(d - GFbias);
332 
333  // initial phase biases - mainly just for output
334  if (arc.ngood == 1)
335  {
336  N1bias = static_cast<long long>(dataP1[i] / wl1 - dataL1[i]);
337  N2bias = static_cast<long long>(dataP2[i] / wl2 - dataL2[i]);
338  }
339  }
340 
341  // fill the Arc as much as possible
342  Arcs.clear();
343  arc.index = 0; // always since we use data[0,size-1]
344  arc.npts = xdata.size();
345  Arcs[arc.index] = arc;
346 
347  // Begin GDC processing ----------------------------------------
348  unique++;
349  // generate strings for output
350  ostringstream oss;
351  oss << "GDC " << setw(3) << unique;
352  tag = oss.str();
353  if (SPSstr.empty())
354  {
355  Epoch endT(beginT);
356  endT += dtlast;
357  oss.str("");
358  oss << tag << " SPS " << setw(4) << arc.npts // ntot
359  << " " << sat << " " << setw(4) << arc.ngood // sat ngood
360  << " 0 " // status flag
361  << printTime(beginT, outfmt) << " " << printTime(endT, outfmt)
362  << " " << setprecision(1) << fixed << dt << " L1 L2 P1 P2";
363  SPSstr = oss.str();
364  }
365  LOG(INFO) << SPSstr;
366  SPSstr = string(); // clear for next call
367 
368  // dump data with tag RAW
369  if (cfg(RAW))
370  {
371  dumpData(LOGstrm, tag + " RAW");
372  }
373 
374  // check that segment is long enough
375  if (flags.size() < 2 * cfg(width))
376  {
377  for (i = 0; i < flags.size(); i++)
378  { // cf. flagBadData()
379  if (flags[i] == OK)
380  {
381  flags[i] = BAD;
382  }
383  }
384  LOG(INFO) << tag
385  << " Pass is too short to analyze: " << flags.size()
386  << " < 2 * window width = " << 2 * cfg(width);
387  }
388 
389  while (1)
390  {
391  // Process WL --------------------------------------
392  iret = ProcessOneCombo(WL);
393  if (iret < 0)
394  {
395  break;
396  }
397 
398  // Process GF --------------------------------------
399  iret = ProcessOneCombo(GF);
400  if (iret < 0)
401  {
402  break; // TD not return?
403  }
404 
405  // Check value of slips found
406 
407  break; // mandatory
408  } // end while(1)
409 
410  // final check
411  iret = FinalCheck();
412  if (iret < 0)
413  {
414  return iret;
415  } // never?
416 
417  // build the return message
418  retMsg = returnMessage();
419  if (cfg_func("verbose"))
420  {
421  DumpArcs("#" + tag + " FIN", "");
422  }
423 
424  // generate editing commands
425  if (cfg(doCmds))
426  {
427  generateCmds(cmds);
428  }
429 
430  // are there ever going to be breaks?
431  // for(i=0; i<breaks.size(); i++)
432  // LOG(WARNING) << " Warning - GDC breaks pass at " << breaks[i];
433 
434  return 0;
435  }
436  catch (Exception& e)
437  {
438  GNSSTK_RETHROW(e);
439  }
440  } // end int gdc::DiscontinuityCorrector(dataL1, dataL2,...)
441 
442  //---------------------------------------------------------------------------------
443  // NB return value == nslips is not used
444  int gdc::ProcessOneCombo(const unsigned which)
445  {
446  try
447  {
448  int iret, nslips(0);
449 
450  /* -------------------------------------------------------------------------
451  first look for gross slips using 1st differences, then compute
452  stats, look for gaps, and fix the slips (WLG GFG) */
453  iret = GrossProcessing(which);
454  if (iret < 0)
455  {
456  return iret;
457  }
458  nslips += iret;
459 
460  /* -------------------------------------------------------------------------
461  now look for small slips using window filter, then compute stats and
462  fix the slips (WLW GFW) */
463  iret = FineProcessing(which);
464  if (iret < 0)
465  {
466  return iret;
467  }
468  nslips += iret;
469 
470  return nslips;
471  }
472  catch (Exception& e)
473  {
474  GNSSTK_RETHROW(e);
475  }
476  } // end ProcessOneCombo()
477 
478  //---------------------------------------------------------------------------------
479  /* process one combo (WL or GF) using 1st differences; called by
480  ProcessOneCombo
481  @return return value of filter() if negative, otherwise number of slips
482  found */
483  int gdc::GrossProcessing(const unsigned which)
484  {
485  try
486  {
487  int i, nslips(0), iret;
488  double limit;
489  string label;
490  vector<FilterHit<double>> filterResults;
491 
492  /* filter using first difference, for gross slips and outliers
493  WL1 GF1 */
494  label = LAB[which] + "1";
495  limit = cfg_func(LAB[which] + "grossStep");
496  iret = filterFirstDiff(which, label, limit, filterResults);
497  if (iret < 0)
498  {
499  return iret;
500  }
501  nslips += iret;
502 
503  // dump filter hits
504  if (cfg_func("debug") > -1)
505  {
506  DumpHits(filterResults, "#" + tag, label, 2);
507  }
508 
509  /* merge 1st difference filter results with Arcs; returns number of new
510  arcs NB i unused */
511  i = mergeFilterResultsIntoArcs(filterResults, which);
512 
513  /* recompute stats in each segment
514  not until window filter - gross slip can use Arc.info.step =
515  FilterHit.step */
516  getArcStats(which);
517 
518  // dump Arcs
519  if (cfg_func(label))
520  {
521  DumpArcs("#" + tag, label, 2);
522  }
523 
524  /* look for gaps > MaxGap, end Arc there, add Arc(BEG) where data
525  resumes */
526  findLargeGaps();
527 
528  /* fix gross slips
529  remove slips that are "size 0" -- do this in vector<FilterHit> ? no
530  NB i unused */
531  i = fixSlips(which);
532 
533  // dump data (WLG GFG)
534  label = LAB[which] + "G";
535  if (cfg_func(label))
536  {
537  dumpData(LOGstrm, tag + " " + label);
538  }
539  return nslips;
540  }
541  catch (Exception& e)
542  {
543  GNSSTK_RETHROW(e);
544  }
545  } // end GrossProcessing()
546 
547  //---------------------------------------------------------------------------------
548  /* process one combo (WL or GF) using window filter; called by
549  ProcessOneCombo
550  @return return value of filter() if negative, otherwise number of slips
551  found */
552  int gdc::FineProcessing(const unsigned which)
553  {
554  try
555  {
556  int i, nslips(0), iret;
557  double limit;
558  string label;
559  vector<FilterHit<double>> filterResults;
560  // filter using the window filter
561  label = LAB[which] + "W"; // WLW or GFW
562  limit = cfg_func(LAB[which] + "fineStep");
563  iret = filterWindow(which, label, limit, filterResults);
564  if (iret < 0)
565  {
566  /* a segment is too small...
567  DumpArcs("#After FilterWindow ",label,2); */
568  return iret;
569  }
570  nslips += iret; // iret >= 1 -- counts BOD
571 
572  // dump filter hits
573  if (cfg_func("debug") > -1)
574  {
575  DumpHits(filterResults, "#" + tag, label, 2);
576  }
577 
578  // merge window filter results with Arcs
579  i = mergeFilterResultsIntoArcs(filterResults, which);
580 
581  /* Recompute stats in each segment.
582  NB Filters define FilterHit.step using their analysis; that step is
583  then copied over to Arc.xxinfo.step in mergeFilterResultsIntoArcs().
584  The first difference step is used to fix gross slips. The window
585  filter step (same as in the window algorithm) can also be used, or
586  you could try to re-compute step using more data. */
587  getArcStats(which);
588 
589  /* fix small slips using stats
590  NB i unused */
591  i = fixSlips(which);
592 
593  // dump Arcs
594  if (cfg_func(label))
595  {
596  DumpArcs("#" + tag, label, 2);
597  }
598 
599  // dump data WLF GFF
600  label = LAB[which] + "F";
601  if (cfg_func(label))
602  {
603  dumpData(LOGstrm, tag + " " + label);
604  }
605 
606  return nslips;
607  }
608  catch (Exception& e)
609  {
610  GNSSTK_RETHROW(e);
611  }
612  } // end FineProcessing()
613 
614  //---------------------------------------------------------------------------------
615  /* filter using first differences, to find gross slips and outliers
616  param which is either WL or GF
617  param label string to be passed to dump e.g. "GF1"
618  param limit pass to filter.setLimit()
619  return return value of filter() if negative (failure),
620  otherwise return number of FilterHit found (>= 1) */
621  int gdc::filterFirstDiff(const unsigned which, const string label,
622  double limit, vector<FilterHit<double>>& hits)
623  {
624  try
625  {
626  // configure first difference filter
627  FirstDiffFilter<double> fdf(xdata, (which == GF ? dataGF : dataWL),
628  flags);
629  fdf.setw(cfg(oswidth));
630  fdf.setprecision(cfg(osprec));
631  fdf.setLimit(limit);
632 
633  // run it
634  int iret = fdf.filter();
635  if (iret < 0)
636  {
637  return iret;
638  }
639 
640  // analyze results
641  iret = fdf.analyze();
642 
643  // compute stats on each segment, then get results to return
644  for (int i = 0; i < fdf.results.size(); i++)
645  fdf.getStats(fdf.results[i]);
646 
647  // NB must do this after getStats()
648  hits = fdf.getResults();
649 
650  // dump filter results - will use stats from getStats() WL1 GF1
651  // fdf.setDumpNoAnal(cfg(debug)>-1);
652  if (cfg_func(label))
653  {
654  fdf.dump(LOGstrm, tag + " " + label);
655  }
656 
657  return iret;
658  }
659  catch (Exception& e)
660  {
661  GNSSTK_RETHROW(e);
662  }
663  } // end gdc::filterFirstDiff()
664 
665  //---------------------------------------------------------------------------------
666  /* filter using window filter
667  param which is either WL or GF
668  param label string to be passed to dump e.g. "GFG"
669  param limit pass to filter.setLimit()
670  param return vector< FilterHit<T> > hits containing all outliers and slips
671  return return value of filter() if negative, otherwise number of slips
672  found */
673  int gdc::filterWindow(const unsigned which, const string label,
674  double limit, vector<FilterHit<double>>& hits)
675  {
676  try
677  {
678  unsigned int i;
679 
680  // configure window filter
681  WindowFilter<double> wf(xdata, (which == GF ? dataGF : dataWL), flags);
682  wf.setWidth(cfg(width));
683  wf.setw(cfg(oswidth));
684  wf.setprecision(cfg(osprec));
685  wf.setMinStep(limit);
686  wf.setTwoSample(which == GF);
687  // wf.setDebug(true); // TEMP
688 
689  // run it
690  int iret = wf.filter();
691  if (iret == -2)
692  {
693  LOG(ERROR) << " Call to GF window filter without time data!";
694  GNSSTK_THROW(
695  Exception("Call to GF window filter without time data"));
696  }
697  else if (iret == -1 || iret == -3)
698  { // segment is too small
699  return iret;
700  }
701 
702  // analyze results
703  iret = wf.analyze();
704 
705  // compute stats on each segment, then get results to return
706  for (int i = 0; i < wf.results.size(); i++)
707  {
708  wf.getStats(wf.results[i]);
709  }
710 
711  // NB this must be after getStats()
712  hits = wf.getResults();
713 
714  // dump filter results - will use stats from getStats()
715  wf.setDumpAnalMsg(cfg(debug) > -1 || cfg(verbose) != 0);
716  if (cfg_func(label))
717  {
718  wf.dump(LOGstrm, tag + " " + label);
719  }
720 
721  /*LOG(INFO) << " There are " << wf.maybes.size() << " maybes";
722  for(i=0; i<wf.maybes.size(); i++) {
723  if(wf.maybes[i].score < 66) continue;
724  if(wf.maybes[i].score < 85) continue;
725  LOG(INFO) << "#" << tag << " " << sat
726  << " poss(" << wf.maybes[i].score << "%)"
727  << " " << LAB[which] << " slip: step "
728  << fixed << setprecision(2) << setw(6) << wf.maybes[i].step <<
729  " wl"
730  << " indx " << wf.maybes[i].index
731  << " " << printTime(xtime(wf.maybes[i].index),outfmt)
732  << (cfg(debug)>-1 ? wf.maybes[i].msg : "");
733  } */
734 
735  return iret;
736  }
737  catch (std::exception& e)
738  {
739  LOG(ERROR) << "std exception " << e.what();
740  GNSSTK_THROW(Exception(string("std exception") + e.what()));
741  }
742  catch (Exception& e)
743  {
744  GNSSTK_RETHROW(e);
745  }
746  } // end int gdc::filterWindow()
747 
748  //---------------------------------------------------------------------------------
749  /* merge filter results (vector<FilterHit>) into the Arcs list, and set
750  flags[]. The merge will mark outliers, add new Arc's where there are
751  slips, and call fixUpArcs(), if necessary. Test with cases where there is
752  huge data rejection in GF, after WL slips, etc. param hits
753  vector<FilterHit> which is results of filter param which is either WL or
754  GF return the number of new Arcs in Arcs */
755  int gdc::mergeFilterResultsIntoArcs(vector<FilterHit<double>>& hits,
756  const unsigned which)
757  {
758  try
759  {
760  // is this necessary? ever used?
761  if (Arcs.size() == 0)
762  {
763  GNSSTK_THROW(Exception("No Arcs found"));
764  }
765  if (hits.size() == 0)
766  {
767  GNSSTK_THROW(Exception("No Filter results found"));
768  }
769 
770  bool fixup(false);
771  int i, narcs(0);
772  unsigned int minpts(cfg(MinPts));
773  double lostslip(0.0);
774 
775  // flag data BAD for new outliers and small segments
776  for (i = 0; i < hits.size(); i++)
777  {
778  LOG(DEBUG) << "#" << tag << " merge " << LAB[which]
779  << " hit into Arc[" << hits[i].index << "] "
780  << hits[i].asString();
781 
782  // hits[i].type can be BOD slip outlier other(never used)
783  if (hits[i].type == FilterHit<double>::BOD)
784  {
785  ; // nothing to do
786  }
787 
788  else if (hits[i].type == FilterHit<double>::outlier)
789  {
790  unsigned int flag = (which == WL ? WLOUTLIER : GFOUTLIER);
791  // mark all the data in this hit
792  flagBadData(hits[i], flag);
793  fixup = true;
794  }
795 
796  // else if(hits[i].type == FilterHit<double>::other) {
797  // never used
798  //}
799 
800  else if (hits[i].type == FilterHit<double>::slip)
801  { // slip
802  // if too short, mark it and don't make an Arc; however
803  // accumulate the slip magnitude for following slips
804  if (hits[i].ngood < minpts)
805  { // too short
806  // mark all the data in this hit
807  flagBadData(hits[i], (which == WL ? WLSHORT : GFSHORT));
808  fixup = true;
809 
810  // save the slip, to add to later slips
811  lostslip += hits[i].step;
812 
813  continue;
814  }
815 
816  // find the Arc in which this hit lies
817  map<int, Arc>::iterator ait = Arcs.begin();
818  findArc(hits[i].index, ait);
819 
820  // is there already an Arc here?
821  if (hits[i].index == ait->second.index)
822  {
823  // already an Arc at this index
824  if (ait->second.mark & SLIP[which])
825  { // already a slip(which) here
826  if (ait->second.mark & FIX[which])
827  { // and its been fixed
828  // happens when gross and fine slip
829  // just remove the fix mark...
830  ait->second.mark ^= FIX[which];
831  }
832  else
833  {
834  GNSSTK_THROW( // marked SLIP but not fixed
835  Exception(
836  "Already marked but unfixed should not happen"));
837  }
838  }
839  else // no slip(which) here
840  {
841  ait->second.mark |= SLIP[which]; // so mark it SLIP
842  }
843 
844  // no need for fixup
845  }
846 
847  // no Arc at this point - add one
848  else
849  {
850  if (addArc(hits[i].index, SLIP[which]))
851  {
852  narcs++;
853  fixup = true;
854  }
855  }
856 
857  // copy hits[i].step into the Arc; this will be used in
858  // fixSlips()
859  (which == WL ? Arcs[hits[i].index].WLinfo.step
860  : Arcs[hits[i].index].GFinfo.step) =
861  lostslip + hits[i].step;
862  (which == WL ? Arcs[hits[i].index].WLinfo.sigma
863  : Arcs[hits[i].index].GFinfo.sigma) =
864  lostslip + hits[i].sigma;
865  lostslip = 0.0;
866 
867  } // end if slip
868 
869  // NB note there is a continue stmt above
870 
871  } // end loop over hits
872 
873  if (fixup)
874  {
875  fixUpArcs(); // recompute points for all Arcs
876  }
877 
878  return narcs;
879  }
880  catch (Exception& e)
881  {
882  GNSSTK_RETHROW(e);
883  }
884  } // end int gdc::mergeFilterResultsIntoArcs()
885 
886  //---------------------------------------------------------------------------------
887  /* Flag bad data in the flags[] array, using a filter hit object. Don't alter
888  Arcs NB FixUpArcs() must be called after this routine to recompute ngood.
889  Note that flags[] is changed ONLY if flags[] currently OK....TD use
890  bitmap?
891  @param hit the FilterHit, which describes a segment to be marked
892  @param flagvalue the value to which flags[] is set for each outlier. */
893  void gdc::flagBadData(const FilterHit<double>& hit, const unsigned flagvalue)
894  {
895  // loop over all the data in this segment (hit)
896  for (unsigned int i = hit.index; i < hit.index + hit.npts; i++)
897  {
898  if (flags[i] == OK)
899  {
900  flags[i] = flagvalue; // TD let flags be a bitmap too?
901  // don't modify Arc.ngood b/c fixUpArcs() will have to be called
902  // anyway
903  }
904  }
905  } // end void gdc::flagBadData()
906 
907  //---------------------------------------------------------------------------------
908  /* add a new Arc to Arcs at index, using the given value for mark. If there
909  is already an Arc at index, instead just assign the mark using &= (its a
910  bitmap). param index the index into data[] at which to add the Arc param
911  mark the value to assign to the new Arc's mark return true if new Arc
912  created -> fixUpArcs() should be called. */
913  bool gdc::addArc(const int index, const unsigned mark)
914  {
915  // find the Arc containing the given index
916  map<int, Arc>::iterator ait = Arcs.begin();
917  findArc(index, ait);
918 
919  // create a new segment
920  Arc B;
921  B.mark = mark;
922  B.index = index;
923  B.npts = ait->second.npts - (index - ait->second.index);
924  B.ngood = 0;
925 
926  // add it
927  Arcs[B.index] = B;
928  // if(cfg(debug)>-1) LOG(INFO) << "# Add new Arc " << B.asString();
929 
930  // modify the existing segment
931  ait->second.npts = index - ait->second.index;
932  ait->second.ngood = 0;
933 
934  // replace it
935  B = ait->second; // B is now the original Arc
936  Arcs.erase(ait);
937  Arcs[B.index] = B;
938 
939  return true; // have to recompute ngood's in fixUpArcs()
940 
941  } // end bool gdc::addArc()
942 
943  //---------------------------------------------------------------------------------
944  /* modify Arcs: recompute npts and ngood, remove empty Arcs
945  TD also remove Arcs with slip == 0 */
946  void gdc::fixUpArcs()
947  {
948  try
949  {
950  map<int, Arc>::iterator ait;
951  vector<int> tmpArcs;
952 
953  // recompute all the Arcs npts and ngood
954  recomputeArcs();
955 
956  // ensure that Arcs[a.index] = a;
957  map<int, Arc> oldArcs(Arcs);
958  Arcs.clear();
959  for (ait = oldArcs.begin(); ait != oldArcs.end(); ++ait)
960  Arcs[ait->second.index] = ait->second;
961  }
962  catch (Exception& e)
963  {
964  GNSSTK_RETHROW(e);
965  }
966  } // end void gdc::fixUpArcs()
967 
968  //---------------------------------------------------------------------------------
969  /* recompute the npts and ngood for each Arc using the indexes in the map:
970  Arcs */
971  void gdc::recomputeArcs()
972  {
973  try
974  {
975  map<int, Arc>::iterator ait, cit;
976 
977  // loop over Arcs, recomputing npts
978  cit = ait = Arcs.begin();
979  while (++cit != Arcs.end())
980  {
981  ait->second.npts = cit->second.index - ait->second.index;
982  ++ait;
983  }
984  ait->second.npts = xdata.size() - ait->second.index; // last one
985 
986  // loop over Arcs, recomputing ngood
987  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
988  {
989  computeNgood(ait->second); // ignore return value
990  // LOG(INFO) << "GDC recomputeArcs " << ait->second.asString();
991  }
992  }
993  catch (Exception& e)
994  {
995  GNSSTK_RETHROW(e);
996  }
997  }
998 
999  //---------------------------------------------------------------------------------
1000  /* compute stats for 'which' data (WL or GF but not both) for the given Arc.
1001  NB this is sneaky and goes across fixed slips....
1002  NB do not confuse this with *Filter::getStats()
1003  NB there is a getArcStats(which) which loops over all Arcs calling this
1004  for each. param iterator pointing to element of Arcs that contains the Arc
1005  param which is either WL or GF */
1006  void gdc::getArcStats(map<int, Arc>::iterator& ait, const unsigned which)
1007  {
1008  StatsFilterBase<double> *ptrStats = nullptr;
1009  try
1010  {
1011  bool isWL(which == WL);
1012  int i, index, npts;
1013  map<int, Arc>::iterator cit(ait);
1014 
1015  if (isWL)
1016  {
1017  ptrStats = new OneSampleStatsFilter<double>();
1018  }
1019  else
1020  {
1021  ptrStats = new TwoSampleStatsFilter<double>();
1022  }
1023  i = index = cit->second.index;
1024  npts = cit->second.npts;
1025 
1026  // loop over continuous data in the arc
1027  while (cit != Arcs.end() && i < xdata.size())
1028  {
1029  // add to stats (xdata is ignored in OneSampleStats)
1030  // don't include bad data, unless this is a REJ arc...
1031  if (flags[i] == OK || cit->second.mark == Arc::REJ)
1032  {
1033  ptrStats->Add(xdata[i], (isWL ? dataWL[i] : dataGF[i]));
1034  }
1035  i++;
1036  if (i == index + npts)
1037  { // reached end of Arc
1038  cit++; // go to the next one..
1039  if (cit == Arcs.end()) // ..unless there isn't one
1040  {
1041  break;
1042  }
1043  if (!(cit->second.mark & SLIP[which])) // ..or its not a slip
1044  {
1045  break;
1046  }
1047  if (!(cit->second.mark & FIX[which])) // ..or its not been fixed
1048  {
1049  break;
1050  }
1051  index = cit->second.index;
1052  npts = cit->second.npts;
1053  }
1054  }
1055 
1056  // store results (N,ave/aveY,sig/sigYX) in the original Arc
1057  (isWL ? ait->second.WLinfo : ait->second.GFinfo).n = ptrStats->N();
1058  (isWL ? ait->second.WLinfo : ait->second.GFinfo).ave =
1059  ptrStats->Average();
1060  (isWL ? ait->second.WLinfo : ait->second.GFinfo).sig =
1061  ptrStats->StdDev();
1062  }
1063  catch (Exception& e)
1064  {
1065  delete ptrStats;
1066  GNSSTK_RETHROW(e);
1067  }
1068  delete ptrStats;
1069  } // end void gdc::getArcStats()
1070 
1071  //---------------------------------------------------------------------------------
1072  /* find gaps within Arc in Arcs; if gap is larger than MaxGap, break the Arc
1073  into two, adding a BEG Arc (unless its at the very end of the data). param
1074  which is either WL or GF */
1075  void gdc::findLargeGaps()
1076  {
1077  try
1078  {
1079  int limit(cfg(MaxGap));
1080  map<int, int> allgaps, gaps;
1081  map<int, int>::const_iterator git;
1082  map<int, Arc>::iterator ait;
1083 
1084  // loop over Arcs
1085  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1086  {
1087  // find all the gaps
1088  gaps = findGaps(ait->second);
1089  if (gaps.size() == 0)
1090  {
1091  continue;
1092  }
1093 
1094  allgaps.insert(gaps.begin(), gaps.end());
1095  }
1096 
1097  /* process the gaps
1098  must do it this way (not within Arcs loop), in case one Arc gets
1099  split twice */
1100  bool fixup(false);
1101  for (git = gaps.begin(); git != gaps.end(); ++git)
1102  {
1103  if (git->second <= limit)
1104  {
1105  continue; // skip small gaps
1106  }
1107 
1108  // find the Arc it belongs to
1109  ait = Arcs.begin();
1110  findArc(git->first, ait);
1111 
1112  if (git->first == ait->second.index)
1113  {
1114  continue; // skip 'gap' at beginning of Arc
1115  }
1116  if (git->first + git->second ==
1117  ait->second.index + ait->second.npts)
1118  {
1119  continue; // skip 'gap' at end of Arc
1120  }
1121 
1122  // Arc at ait must be split // we don't need fixUp
1123  addArc(git->first + git->second, Arc::BEG);
1124 
1125  // must recompute ngood, but only for one Arc .. oh well
1126  ait = Arcs.begin(); // iterator is corrupted by addArc
1127  findArc(git->first + git->second, ait);
1128  computeNgood(ait->second);
1129  }
1130  if (fixup)
1131  {
1132  fixUpArcs(); // recompute points for all Arcs
1133  }
1134  }
1135  catch (Exception& e)
1136  {
1137  GNSSTK_RETHROW(e);
1138  }
1139  } // end void gdc::findLargeGaps()
1140 
1141  //---------------------------------------------------------------------------------
1142  /* find gaps within the given Arc, including those at the very beginning
1143  (index==Arc.index) and at the very end
1144  (index+nptsgap==Arc.index+Arc.npts). return map with key=index of
1145  beginning of gap, value=number of points in gap. */
1146  map<int, int> gdc::findGaps(const Arc& arc)
1147  {
1148  try
1149  {
1150  unsigned int i, count,
1151  index; // count consecutive bad pts, starting at index
1152  map<int, int> gaps;
1153  for (count = 0, i = arc.index; i < arc.index + arc.npts; ++i)
1154  {
1155  if (flags[i] == OK)
1156  { // good
1157  if (count > 0)
1158  { // is there a gap ending here?
1159  gaps[index] = count;
1160  count = 0;
1161  }
1162  }
1163  else
1164  { // bad - add to count
1165  if (count == 0)
1166  {
1167  index = i;
1168  }
1169  count++;
1170  }
1171  }
1172  if (count > 0)
1173  { // is there a gap at the very end?
1174  gaps[index] = count;
1175  }
1176 
1177  return gaps;
1178  }
1179  catch (Exception& e)
1180  {
1181  GNSSTK_RETHROW(e);
1182  }
1183  } // end map<int,int> gdc::findGaps(const Arc& arc) throw(Exception)
1184 
1185  //---------------------------------------------------------------------------------
1186  /* fix slips between Arcs, using info.step (NOT info.Nslip), which is defined
1187  by the filter in results(FilterHit). Compute an integer from step and ADD
1188  it to Nslip. Thus Nslip has the total slip, while step has only the last
1189  estimate, used to fix. In the case of the FirstDifferenceFilter this step
1190  is only an approximate fix; for the WindowFilter the step is defined by
1191  stats on the two segments (one-sample for WL and two-sample for GF). param
1192  which is either WL or GF return the number of slips fixed. */
1193  int gdc::fixSlips(const unsigned which)
1194  {
1195  try
1196  {
1197  int i, nslips(0);
1198  long N;
1199  double step, istep;
1200  static const double GFfactor(wl2 / wlGF);
1201  static const double IFfactor(
1202  isGLO ? 3.5 : 3.52941176470588); // TD what is this?
1203 
1204  // loop over Arcs using iterator ait //, with dummy copy cit
1205  map<int, Arc>::iterator ait; //, cit;
1206  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1207  {
1208  if ((ait->second.mark & SLIP[which]) == 0) // its not a slip
1209  {
1210  continue;
1211  }
1212  if ((ait->second.mark & FIX[which]) != 0) // its been fixed
1213  {
1214  continue;
1215  }
1216 
1217  LOG(DEBUG) << "#" << tag << " fix slip for Arc[" << ait->first
1218  << "] " << ait->second.asString();
1219 
1220  // first get the step in wavelengths for this slip (not the total)
1221  step = (which == WL ? ait->second.WLinfo.step
1222  : ait->second.GFinfo.step);
1223 
1224  // get the integer-wavelengths step, and use this for fixing
1225  N = long(step + (step > 0.0 ? 0.5 : -0.5));
1226  istep = static_cast<double>(N);
1227 
1228  /* accumulate Nslip, which is the total slip
1229  NB b/c gross and fine both contribute, and want the total slip at
1230  the end */
1231  (which == WL ? ait->second.WLinfo.Nslip
1232  : ait->second.GFinfo.Nslip) += N;
1233 
1234  // mark it fixed, TD but if N==0 shouldn't you remove the SLIP?
1235  ait->second.mark |= FIX[which];
1236 
1237  /* if it's a non-zero step, modify the data from here all the way
1238  out */
1239  if (N)
1240  {
1241  nslips++; // count it
1242  if (which == WL)
1243  {
1244  for (i = ait->second.index; i < xdata.size(); i++)
1245  {
1246  dataWL[i] -= istep;
1247  dataGF[i] -= GFfactor * istep;
1248  }
1249  }
1250  else
1251  {
1252  for (i = ait->second.index; i < xdata.size(); i++)
1253  {
1254  dataGF[i] -= istep;
1255  }
1256  }
1257 
1258  // save 'learn' message here, with time tag, which, step, etc
1259 
1260  } // end if N
1261  }
1262 
1263  return nslips;
1264  }
1265  catch (Exception& e)
1266  {
1267  GNSSTK_RETHROW(e);
1268  }
1269  } // end int gdc::fixSlips(const unsigned which) throw(Exception)
1270 
1271  // do a final check on the pass. Look for isolated good points (< MinPts good
1272  // points surrounded by N(?) bad points on each side.
1273  int gdc::FinalCheck()
1274  {
1275  try
1276  {
1277  bool fixup(false);
1278  unsigned int i, k;
1279  int iret(0), j, currstate, ngood, nbad;
1280  vector<int> gbs;
1281 
1282  // look for segments of N good (+) or bad (-) points
1283  currstate = ngood = nbad = 0;
1284  for (i = 0; i < xdata.size(); i++)
1285  {
1286  if (flags[i] == 0)
1287  {
1288  if (currstate == -1)
1289  { // end of bad segment
1290  gbs.push_back(-nbad);
1291  nbad = ngood = 0;
1292  }
1293  currstate = 1;
1294  ngood++;
1295  }
1296  else
1297  {
1298  if (currstate == 1)
1299  { // end of good segment
1300  gbs.push_back(ngood);
1301  nbad = ngood = 0;
1302  }
1303  currstate = -1;
1304  nbad++;
1305  }
1306  }
1307 
1308  // look for isolated good segments
1309  int min = cfg(MinPts);
1310  if (min > 10)
1311  {
1312  min = 10;
1313  }
1314  int gap = cfg(MaxGap);
1315  if (gap > 10)
1316  {
1317  gap = 10;
1318  }
1319 
1320  // only do for >3 segments
1321  // k is current index, used for marking
1322  if (gbs.size() > 3)
1323  {
1324  for (k = 0, i = 0; i < gbs.size(); i++)
1325  {
1326  if (gbs[i] > 0 && gbs[i] < min)
1327  { // current segment is good and short
1328  // can now assume gbs[i-1] and gbs[i+1] are < 0 (bad)
1329 
1330  // prev segment is first or last
1331  bool prev1(i == 1);
1332  // next segment is last
1333  bool next1(i == gbs.size() - 2);
1334  // prev segment is big and bad, or not there
1335  bool prevbb(i == 0 || -gbs[i - 1] > gap);
1336  // next segment is big and bad, or not there
1337  bool nextbb(i == gbs.size() - 1 || -gbs[i + 1] > gap);
1338  // prev is 2nd or 3rd, and 1st or 2nd is small
1339  bool prev23((i == 2 && gbs[0] < gap) ||
1340  (i == 3 && gbs[1] < gap));
1341  // next is 2nd or 3rd from end, and 1st or 2nd is small
1342  bool next23((i == gbs.size() - 1 && -gbs[i - 1] < gap) ||
1343  (i == gbs.size() - 2 && -gbs[i + 1] < gap) ||
1344  (i == gbs.size() - 3 && -gbs[i + 1] < gap));
1345  bool ibeg(i == 0 && -gbs[1] > gap);
1346  bool iend(i == gbs.size() - 1 && -gbs[i - 1] > gap);
1347 
1348  if ((prevbb && nextbb) || prev1 || next1 || prev23 ||
1349  next23 || ibeg || iend)
1350  {
1351  fixup = true;
1352  // oss2 << " Mark:";
1353  for (j = 0; j < gbs[i]; j++)
1354  {
1355  flags[k + j] = ISOLATED;
1356  // oss2 << " " << k+j;
1357  }
1358  }
1359  } // end if long and bad
1360 
1361  // keep k pointed to start of this segment
1362  k += ::abs(gbs[i]);
1363  }
1364  }
1365 
1366  if (fixup)
1367  {
1368  fixUpArcs(); // recompute points for all Arcs
1369  }
1370 
1371  if (cfg_func("FIN"))
1372  {
1373  dumpData(LOGstrm, tag + " FIN");
1374  }
1375 
1376  return iret;
1377  }
1378  catch (Exception& e)
1379  {
1380  GNSSTK_RETHROW(e);
1381  }
1382  }
1383 
1384  // end of processing routines
1385 
1386  //---------------------------------------------------------------------------------
1387  void gdc::dumpData(ostream& os, const string msg)
1388  {
1389  size_t i, k;
1390  map<int, Arc>::const_iterator kt = Arcs.begin();
1391 
1392  for (k = 0, i = 0; i < xdata.size(); i++)
1393  {
1394  string arcmsg;
1395  if (kt != Arcs.end() && kt->second.index == i)
1396  {
1397  arcmsg = string(" ") + kt->second.asString();
1398  kt++;
1399  }
1400 
1401  os << msg << " " << sat << " " << printTime(xtime(i), outfmt) << fixed
1402  << setprecision(3) << " " << setw(9)
1403  << (xtime(i) - beginT) // xdata[i]
1404  << " " << setw(2) << flags[i] << fixed << setprecision(4) << " "
1405  << setw(9) << dataWL[i] * wlWL << " " << setw(9) << dataGF[i] * wlGF
1406  << arcmsg << endl;
1407  }
1408  } // end void gdc::dumpData(ostream& os, const string msg)
1409 
1410  //---------------------------------------------------------------------------------
1411  void gdc::DumpHits(const vector<FilterHit<double>>& filterResults,
1412  const string& tag, const string& label, int prec)
1413  {
1414  if (prec == -1)
1415  {
1416  prec = cfg(osprec);
1417  }
1418  for (int i = 0; i < filterResults.size(); i++)
1419  {
1420  LOG(INFO) << tag << " " << label << " Hit" << i + 1 << "["
1421  << filterResults[i].index << "] "
1422  << filterResults[i].asStatsString(prec);
1423  }
1424  }
1425 
1426  //---------------------------------------------------------------------------------
1427  void gdc::DumpArcs(const string& tag, const string& label, int prec)
1428  {
1429  if (prec == -1)
1430  {
1431  prec = cfg(osprec);
1432  }
1433  int i(1);
1434  map<int, Arc>::const_iterator ait;
1435  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1436  {
1437  LOG(INFO) << tag << (label.size() ? " " + label : "") << " Arc" << i++
1438  << "[" << ait->first << "] " << ait->second.asString(prec);
1439  }
1440  } // end void gdc::DumpArcs()
1441 
1442  //---------------------------------------------------------------------------------
1443  // build the string that is returned by the discontinuity corrector
1444  string gdc::returnMessage(int prec, int wid)
1445  {
1446  int minpts(cfg(MinPts));
1447  string retmsg;
1448  ostringstream oss, oss2;
1449  map<int, Arc>::iterator ait;
1450 
1451  if (prec == -1)
1452  {
1453  prec = cfg(osprec);
1454  }
1455  if (wid == -1)
1456  {
1457  wid = cfg(oswidth);
1458  }
1459  oss << fixed << setprecision(prec);
1460  oss2 << fixed << setprecision(prec);
1461 
1462  // Find segs (>MinPts) w/ ngood=0, call "REJ" Arcs; recompute stats
1463  // each Arc can potentially be broken into three: REJ,Arc,REJ
1464  map<int, Arc> newArcs;
1465  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1466  {
1467  int n, i, ib(ait->second.index),
1468  ie(ait->second.index + ait->second.npts);
1469 
1470  // first part of Arc
1471  for (n = 0, i = ib; i < ie; i++)
1472  {
1473  if (flags[i] == 0)
1474  {
1475  break;
1476  }
1477  else
1478  {
1479  n++;
1480  }
1481  }
1482 
1483  // if the entire Arc is bad data, just relabel it REJ
1484  if (n == ait->second.npts)
1485  {
1486  ait->second.mark = Arc::REJ;
1487  continue;
1488  }
1489 
1490  // if first part of Arc is >minpts of bad data, call it a REJ Arc
1491  if (n > minpts)
1492  {
1493  Arc A;
1494  A.mark = Arc::REJ;
1495  A.index = 0;
1496  A.npts = n;
1497  A.ngood = 0;
1498  newArcs[A.index] = A;
1499  ait->second.index = ib + n;
1500  }
1501 
1502  // last part of Arc
1503  for (n = 0, i = ie - 1; i >= ib; i--)
1504  {
1505  if (flags[i] == 0)
1506  {
1507  break;
1508  }
1509  else
1510  {
1511  n++;
1512  }
1513  }
1514  if (n > minpts)
1515  {
1516  Arc A;
1517  A.mark = Arc::REJ;
1518  A.index = ie - n;
1519  A.npts = n;
1520  A.ngood = 0;
1521  newArcs[A.index] = A;
1522  }
1523  }
1524 
1525  if (newArcs.size())
1526  {
1527  // add existing ones to newArcs, changing the key to (possibly new)
1528  // index
1529  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1530  {
1531  newArcs[ait->second.index] = ait->second;
1532  }
1533  // now remake Arcs
1534  Arcs.clear();
1535  for (ait = newArcs.begin(); ait != newArcs.end(); ++ait)
1536  {
1537  Arcs[ait->second.index] = ait->second;
1538  }
1539  recomputeArcs();
1540  }
1541  // recompute stats for all arcs // TD better to do it before this
1542  getArcStats(WL);
1543  getArcStats(GF);
1544 
1545  // TD add a SUM line giving number of slips, fixes, etc.
1546 
1547  // loop over the Arcs, converting each to a line of the message
1548  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1549  {
1550  const Arc A(ait->second);
1551  const unsigned mark(A.mark);
1552 
1553  oss.str("");
1554  oss2.str("");
1555  if (mark & Arc::BEG)
1556  {
1557  oss << "BEG";
1558  }
1559  else if ((mark & Arc::WLSLIP) || (mark & Arc::GFSLIP))
1560  {
1561  oss << "FIX";
1562  oss2 << " n(WL,GF) " << A.WLinfo.Nslip << "," << A.GFinfo.Nslip;
1563  }
1564  else if (mark & Arc::REJ)
1565  {
1566  oss << "REJ";
1567  // NB REJ'd arcs can still hold slips
1568  if (A.WLinfo.Nslip != 0 || A.GFinfo.Nslip != 0)
1569  {
1570  oss2 << " n(WL,GF) " << A.WLinfo.Nslip << "," << A.GFinfo.Nslip;
1571  }
1572  }
1573  oss << " " << setw(4) << A.index << " "
1574  << printTime(xtime(A.index), outfmt) << " " << setw(4) << A.npts
1575  << " " << setw(4) << A.ngood;
1576  if (A.WLinfo.n > 0)
1577  {
1578  oss << " WL " << setw(4) << A.WLinfo.n << " " << setw(wid)
1579  << A.WLinfo.ave << " +- " << setw(wid) << A.WLinfo.sig;
1580  }
1581  if (A.GFinfo.n > 0)
1582  {
1583  oss << " GF " << setw(4) << A.GFinfo.n << " " << setw(wid)
1584  << A.GFinfo.ave << " +- " << setw(wid) << A.GFinfo.sig;
1585  }
1586  oss << oss2.str();
1587 
1588  retmsg += oss.str() + "\n";
1589  }
1590  StringUtils::stripTrailing(retmsg, '\n');
1591 
1592  return retmsg;
1593  }
1594 
1595  //---------------------------------------------------------------------------------
1596  /* apply the results to fix the input SatPass cf. cfg(doFix)
1597  param SP SatPass object containing the input data.
1598  param breaks vector of indexes where SatPass SP must be broken into two
1599  param marks vector of indexes in SatPass SP where breaks are suspected */
1600  void gdc::applyFixesToSatPass(SatPass& SP)
1601  {
1602  try
1603  {
1604  unsigned int i, j;
1605  long long nGF, nWL, nL1, nL2;
1606  Epoch ttag, tbeg, tend;
1607  map<int, Arc>::iterator ait;
1608 
1609  double dL1, dL2; // just double versions of nL1, nL2
1610  const string L1("L1"), L2("L2");
1611 
1612  nL1 = nL2 = 0; // running total slip
1613  ait = Arcs.begin();
1614  for (i = 0; i < SP.size(); i++)
1615  {
1616  if (ait != Arcs.end() && i == ait->first)
1617  { // at new arc
1618  Arc &arc(ait->second);
1619 
1620  LOG(DEBUG) << "#" << tag << " applyFix with Arc[" << ait->first
1621  << "] " << ait->second.asString();
1622 
1623  if ((arc.mark & (Arc::WLSLIP | Arc::GFSLIP)) ||
1624  (arc.mark & Arc::REJ))
1625  {
1626  // redefine biases nL1 nL2
1627  nGF = arc.GFinfo.Nslip;
1628  nWL = arc.WLinfo.Nslip;
1629  // real slips do accumulate here
1630  nL1 -= nGF; // b/c Ngf(corrected) = -N1
1631  nL2 -= (nGF + nWL); // b/c Nwl = N1-N2
1632  dL1 = static_cast<double>(nL1);
1633  dL2 = static_cast<double>(nL2);
1634  }
1635 
1636  /* if((arc.mark & Arc::WLMARK) || (arc.mark & Arc::GFMARK)) {
1637  marks.push_back(i);
1638  NO continue;
1639  } */
1640 
1641  if (arc.mark & Arc::REJ)
1642  {
1643  // reject all the data in this Arc
1644  for (j = 0; j < arc.npts; j++)
1645  {
1646  SP.setFlag(i + j, SatPass::BAD);
1647  if (cfg(UserFlag))
1648  {
1649  SP.setUserFlag(i + j, cfg(UserFlag));
1650  }
1651  }
1652  }
1653 
1654  if (arc.mark & Arc::BEG)
1655  {
1656  if (i != 0)
1657  {
1658  LOG(WARNING) << " Warning - GDC breaks pass at index " << i
1659  << " sat " << sat << " time "
1660  << printTime(xtime(i), outfmt);
1661  }
1662  }
1663 
1664  // increment ait, prep for next arc
1665  ++ait;
1666  }
1667 
1668  if (flags[i] == BAD) // nothing to do - SatPass set before call
1669  {
1670  continue;
1671  }
1672 
1673  if (flags[i] != OK)
1674  {
1675  SP.setFlag(i, SatPass::BAD);
1676  if (cfg(UserFlag))
1677  {
1678  SP.setUserFlag(i, cfg(UserFlag));
1679  }
1680  }
1681  else
1682  {
1683  if (nL1)
1684  {
1685  SP.data(i, L1) -= dL1;
1686  }
1687  if (nL2)
1688  {
1689  SP.data(i, L2) -= dL2;
1690  }
1691  }
1692  } // end loop over data in SP
1693  }
1694  catch (Exception& e)
1695  {
1696  GNSSTK_RETHROW(e);
1697  }
1698  }
1699 
1700  //---------------------------------------------------------------------------------
1701  /* apply the results to generate editing commands; cfg(doCmds)
1702  Use tk-RinEdit form for commands (--IF name, etc) since EditRinex also
1703  takes.
1704  @param cmds vector of strings giving editing commands for RINEX
1705  editor. */
1706  void gdc::generateCmds(vector<string>& cmds)
1707  {
1708  try
1709  {
1710  unsigned int i, j, k;
1711  long long nGF, nWL, nL1, nL2;
1712  Epoch ttag, tbeg, tend;
1713  map<int, Arc>::iterator ait;
1714  static const string L1(cfg(doRINEX3) ? "L1C" : "L1"),
1715  L2(cfg(doRINEX3) ? "L2W" : "L2");
1716 
1717  // generate commands
1718  ostringstream oss;
1719  oss << "--BD+ " << sat << "," << L1 << ","
1720  << printTime(beginT, "%Y,%m,%d,%H,%M,%S,") << N1bias
1721  << printTime(beginT, " # initial L1 bias at %F,%.3g");
1722  cmds.push_back(oss.str());
1723  oss.str("");
1724  oss << "--BD+ " << sat << "," << L2 << ","
1725  << printTime(beginT, "%Y,%m,%d,%H,%M,%S,") << N2bias
1726  << printTime(beginT, " # initial L2 bias at %F,%.3g");
1727  cmds.push_back(oss.str());
1728  oss.str("");
1729 
1730  for (ait = Arcs.begin(); ait != Arcs.end(); ++ait)
1731  {
1732  Arc &arc(ait->second);
1733 
1734  // apply slips -- REJ can store a slip - see karr0880.10o pass 7
1735  if ((arc.mark & (Arc::WLSLIP | Arc::GFSLIP)) ||
1736  (arc.mark & Arc::REJ))
1737  {
1738  nGF = arc.GFinfo.Nslip;
1739  nWL = arc.WLinfo.Nslip;
1740  // slips don't accumulate here, but editing commands do
1741  nL1 = -nGF; // b/c Ngf(corrected) = -N1
1742  nL2 = -nGF - nWL; // b/c Nwl = N1-N2
1743  ttag = xtime(arc.index);
1744  if (nL1)
1745  {
1746  oss << "--BD+ " << sat << "," << L1 << ","
1747  << printTime(ttag, "%Y,%m,%d,%H,%M,%S,") << -nL1
1748  << printTime(ttag, " # L1 slip at %F,%.3g");
1749  cmds.push_back(oss.str());
1750  oss.str("");
1751  }
1752  if (nL2)
1753  {
1754  oss << "--BD+ " << sat << "," << L2 << ","
1755  << printTime(ttag, "%Y,%m,%d,%H,%M,%S,") << -nL2
1756  << printTime(ttag, " # L2 slip at %F,%.3g");
1757  cmds.push_back(oss.str());
1758  oss.str("");
1759  }
1760  }
1761 
1762  // delete entire segment
1763  if (arc.mark & Arc::REJ)
1764  {
1765  tbeg = xtime(arc.index);
1766  tend = xtime(arc.index + arc.npts - 1);
1767  tend += dt; // NB DD- means stop here, don't do this one
1768  if (arc.npts == 1)
1769  {
1770  oss << "--DD " << sat << "," << L1 << ","
1771  << printTime(
1772  tbeg,
1773  "%Y,%m,%d,%H,%M,%S # delete outlier at %F,%.3g");
1774  cmds.push_back(oss.str());
1775  oss.str("");
1776  oss << "--DD " << sat << "," << L2 << ","
1777  << printTime(
1778  tbeg,
1779  "%Y,%m,%d,%H,%M,%S # delete outlier at %F,%.3g");
1780  cmds.push_back(oss.str());
1781  oss.str("");
1782  }
1783  else
1784  {
1785  oss << "--DD+ " << sat << "," << L1 << ","
1786  << printTime(tbeg, "%Y,%m,%d,%H,%M,%S # from %F,%.3g")
1787  << " - delete entire segment = " << arc.npts << " epochs";
1788  cmds.push_back(oss.str());
1789  oss.str("");
1790  oss << "--DD- " << sat << "," << L1 << ","
1791  << printTime(tend, "%Y,%m,%d,%H,%M,%S # to %F,%.3g");
1792  cmds.push_back(oss.str());
1793  oss.str("");
1794  oss << "--DD+ " << sat << "," << L2 << ","
1795  << printTime(tbeg, "%Y,%m,%d,%H,%M,%S # from %F,%.3g");
1796  cmds.push_back(oss.str());
1797  oss.str("");
1798  oss << "--DD- " << sat << "," << L2 << ","
1799  << printTime(tend, "%Y,%m,%d,%H,%M,%S # to %F,%.3g");
1800  cmds.push_back(oss.str());
1801  oss.str("");
1802  }
1803 
1804  continue;
1805  }
1806 
1807  // if there are no outliers, done
1808  if (arc.ngood == arc.npts)
1809  {
1810  continue;
1811  }
1812 
1813  // now run over the data in this Arc looking for outliers
1814  bool bad = false;
1815  j = 0;
1816  k = arc.index + arc.npts;
1817  for (i = arc.index; i < k; i++)
1818  {
1819  if (!bad)
1820  {
1821  if (flags[i] == OK)
1822  {
1823  continue;
1824  }
1825  j = i;
1826  bad = true;
1827  }
1828 
1829  if (bad && (flags[i] == OK || i == k - 1))
1830  {
1831  if ((flags[i] == OK && i == j + 1) || (i == k - 1 && i == j))
1832  {
1833  // isolated outlier
1834  ttag = xtime(j);
1835  oss << "--DD " << sat << "," << L1 << ","
1836  << printTime(
1837  ttag,
1838  "%Y,%m,%d,%H,%M,%S # delete outlier at %F,%.3g");
1839  cmds.push_back(oss.str());
1840  oss.str("");
1841  oss << "--DD " << sat << "," << L2 << ","
1842  << printTime(
1843  ttag,
1844  "%Y,%m,%d,%H,%M,%S # delete outlier at %F,%.3g");
1845  cmds.push_back(oss.str());
1846  oss.str("");
1847  }
1848  else
1849  { // more than one outlier
1850  ttag = xtime(j);
1851  oss << "--DD+ " << sat << "," << L1 << ","
1852  << printTime(ttag, "%Y,%m,%d,%H,%M,%S # delete "
1853  "outliers starting at %F,%.3g");
1854  cmds.push_back(oss.str());
1855  oss.str("");
1856  oss << "--DD+ " << sat << "," << L2 << ","
1857  << printTime(ttag, "%Y,%m,%d,%H,%M,%S # delete "
1858  "outliers starting at %F,%.3g");
1859  cmds.push_back(oss.str());
1860  oss.str("");
1861 
1862  ttag = xtime(i == k - 1 ? i : i - 1);
1863  ttag += dt;
1864  oss << "--DD- " << sat << "," << L1 << ","
1865  << printTime(ttag, "%Y,%m,%d,%H,%M,%S # end deleting "
1866  "outliers at %F,%.3g");
1867  cmds.push_back(oss.str());
1868  oss.str("");
1869  oss << "--DD- " << sat << "," << L2 << ","
1870  << printTime(ttag, "%Y,%m,%d,%H,%M,%S # end deleting "
1871  "outliers at %F,%.3g");
1872  cmds.push_back(oss.str());
1873  oss.str("");
1874  }
1875  bad = false;
1876  }
1877  }
1878 
1879  } // end loop over Arcs
1880  }
1881  catch (Exception& e)
1882  {
1883  GNSSTK_RETHROW(e);
1884  }
1885  }
1886 
1887  //---------------------------------------------------------------------------------
1888  // configuration
1889  //---------------------------------------------------------------------------------
1890  /* Set a parameter in the configuration; the input string 'cmd' is of the
1891  form
1892  '[--DC]<id><s><value>' : separator s is one of ':=,'; leading --DC is
1893  optional. */
1894  bool gdc::setParameter(std::string cmd)
1895  {
1896  try
1897  {
1898  if (cmd.empty())
1899  {
1900  return false;
1901  }
1902  // remove leading --DC
1903  while (cmd[0] == '-')
1904  cmd.erase(0, 1);
1905  if (cmd.substr(0, 2) == "DC")
1906  {
1907  cmd.erase(0, 2);
1908  }
1909 
1910  string label, value;
1911  string::size_type pos = cmd.find_first_of(",=:");
1912  if (pos == string::npos)
1913  {
1914  label = cmd;
1915  }
1916  else
1917  {
1918  label = cmd.substr(0, pos);
1919  value = cmd;
1920  value.erase(0, pos + 1);
1921  }
1922 
1923  return setParameter(label, StringUtils::asDouble(value));
1924  }
1925  catch (Exception& e)
1926  {
1927  GNSSTK_RETHROW(e);
1928  }
1929  } // end bool gdc::setParameter(string cmd) throw(Exception)
1930 
1931  //---------------------------------------------------------------------------------
1932  /* Set a parameter in the configuration using the label and the value,
1933  for booleans use (T,F)=(non-zero,zero). */
1934  bool gdc::setParameter(std::string label, double value)
1935  {
1936  if (CFG.find(label) == CFG.end())
1937  {
1938  return false; // GNSSTK_THROW(Exception("Unknown configuration label " +
1939  }
1940  // label));
1941 
1942  CFG[label] = value;
1943 
1944  // if debug is turned on, turn on some/all of output as well
1945  if (CFG["debug"] > -1)
1946  {
1947  LOG(DEBUG) << "Set GDC " << label << " to " << value;
1948  }
1949 
1950  /* Turn output on/off if debug is being set
1951  NB this does NOT set LOG(DEBUG) */
1952  if (label == string("debug"))
1953  {
1954  // first return to default
1955  CFG["RAW"] = CFG["WLF"] = CFG["GFF"] = CFG["WL1"] = CFG["WLG"] =
1956  CFG["WLW"] = CFG["WLF"] = CFG["GF1"] = CFG["GFG"] = CFG["GFW"] =
1957  CFG["GFF"] = CFG["FIN"] = 0;
1958  if (value > -1)
1959  {
1960  CFG["verbose"] = 1; // debug implies verbose
1961  CFG["WLF"] = 1; // WL after fixing
1962  CFG["GFF"] = 1; // GF after fixing
1963  CFG["FIN"] = 1; // after final check
1964  LOG(INFO) << "GDC:debug sets GDC to output fixed data WLF GFF";
1965  if (CFG["debug"] > 0)
1966  {
1967  CFG["RAW"] = 1; // data (WL,GF) before any processing
1968  LOG(INFO) << "GDC:debug sets GDC to output RAW data";
1969  }
1970  if (CFG["debug"] > 1)
1971  {
1972  CFG["WL1"] = 1; // results of 1st diff filter on WL
1973  CFG["WLG"] = 1; // WL after fixing gross slips (after fdif)
1974  CFG["GF1"] = 1; // results of 1st diff filter on GF
1975  CFG["GFG"] = 1; // GF after fixing gross slips (after fdif)
1976  LOG(INFO)
1977  << "GDC:debug sets GDC to output 1st diff fixes WL1, GF1";
1978  LOG(INFO) << "GDC:debug sets GDC to output gross fixes WLG, GFG";
1979  }
1980  if (CFG["debug"] > 2)
1981  {
1982  CFG["WLW"] = 1; // results of window filter on WL
1983  CFG["GFW"] = 1; // results of window filter on GF
1984  LOG(INFO)
1985  << "GDC:debug sets GDC to output window filters WLW, GFW";
1986  }
1987  }
1988  }
1989  if (label == string("verbose"))
1990  {
1991  if (value != 0)
1992  {
1993  CFG["verbose"] = 1;
1994  }
1995  else
1996  {
1997  CFG["verbose"] = 0;
1998  }
1999  }
2000 
2001  return true;
2002  } // end bool gdc::setParameter(string label, double value)
2003 
2004  //---------------------------------------------------------------------------------
2005  /* Print help page, including descriptions and current values of all
2006  the parameters, to the ostream. */
2007  void gdc::DisplayParameterUsage(ostream& os, string tag, bool advanced)
2008  {
2009  try
2010  {
2011  static const unsigned name_val_width(18), adv_name_val_width(18);
2012  os << tag << "GNSSTk Discontinuity Corrector (GDC) v." << GDCVersion
2013  << " configuration:" << endl;
2014 
2015  string name;
2016  // map<string,double>::const_iterator it;
2017  map<int, string>::const_iterator it;
2018  for (it = CFGlist.begin(); it != CFGlist.end(); it++)
2019  {
2020  name = it->second;
2021  if (CFGdesc[name][0] == '*') // advanced options
2022  {
2023  continue;
2024  }
2025  ostringstream stst;
2026  stst << name // label
2027  << "=" << CFG[name]; // value
2028  os << tag << " "
2029  << StringUtils::leftJustify(stst.str(), name_val_width) << " : "
2030  << CFGdesc[name] // description
2031  << endl;
2032  }
2033  if (advanced)
2034  {
2035  os << tag << " Advanced options :\n";
2036  for (it = CFGlist.begin(); it != CFGlist.end(); it++)
2037  {
2038  name = it->second;
2039  if (CFGdesc[name][0] != '*') // ordinary options
2040  {
2041  continue;
2042  }
2043  ostringstream stst;
2044  stst << name // label
2045  << "=" << CFG[name]; // value
2046  os << tag << " "
2047  << StringUtils::leftJustify(stst.str(), adv_name_val_width)
2048  << " : " << CFGdesc[name].substr(2) // description
2049  << endl;
2050  }
2051  }
2052  }
2053  catch (Exception& e)
2054  {
2055  GNSSTK_RETHROW(e);
2056  }
2057  } // end void gdc::DisplayParameterUsage(ostream& os, bool advanced)
2058 
2059 //------------------------------------------------------------------------------------
2060 /* use to define configuration members with value and descriptive string
2061  NB #a expands to string("a"); note that the parameters are just strings ==
2062  keys. */
2063 #define setcfg(a, b, c) \
2064  { \
2065  CFG[#a] = b; \
2066  CFGdesc[#a] = c; \
2067  CFGlist[CFGindex] = #a; \
2068  CFGindex++; \
2069  }
2070 
2071  /* initialize with default values - only place setcfg is used
2072  Notes:
2073  Don't make grossStep too small, or noise will -> many 'slips' in tiny
2074  segments Note that units of step limits are wavelengths: internally wl,
2075  output meters. */
2076  void gdc::init()
2077  {
2078  try
2079  {
2080  unique = 0; // unique number for each call of DiscCorr(), for output
2081  CFGindex = 0; // just a count of configuration members
2082 
2083  // name, value, "description" NB "* description" makes it
2084  // 'advanced'
2085  setcfg(MaxGap, 10, "maximum allowed gap within a segment (points)");
2086  setcfg(MinPts, 10,
2087  "minimum number of good points in phase segment (points)");
2088  setcfg(width, 20, "* sliding window width (points)");
2089  // WL
2090  setcfg(WLgrossStep, 6.0, "WL gross slip detection threshold (WLwl)");
2091  setcfg(WLfineStep, 0.7, "WL fine slip detection threshold (WLwl)");
2092  // GF
2093  setcfg(GFgrossStep, 6.0, "GF gross slip detection threshold (GFwl)");
2094  setcfg(GFfineStep, 0.7, "GF fine slip detection threshold (GFwl)");
2095  // I/O
2096  setcfg(oswidth, 7, "output stream width (chars)");
2097  setcfg(osprec, 3, "output stream precision (chars)");
2098  setcfg(debug, -1,
2099  "level of diagnostic output, from -1(none) to 3(all)");
2100  setcfg(verbose, 0, "output analysis message in window filter");
2101 
2102  // types of labeled output
2103  setcfg(RAW, 0,
2104  "* output data (WL,GF) before any processing (m) [0=don't]");
2105 
2106  setcfg(WL1, 0,
2107  "* output results of 1st diff filter on WL (wl) [0=don't]");
2108  setcfg(WLG, 0, "* output WL after fixing gross slips (m) [0=don't]");
2109  setcfg(WLW, 0,
2110  "* output results of window filter on WL (wl) [0=don't]");
2111  setcfg(WLF, 0, "* output WL after fixing (m) [0=don't]");
2112 
2113  setcfg(GF1, 0,
2114  "* output results of 1st diff filter on GF (wl) [0=don't]");
2115  setcfg(GFG, 0, "* output GF after fixing gross slips (m) [0=don't]");
2116  setcfg(GFW, 0,
2117  "* output results of window filter on GF (wl) [0=don't]");
2118  setcfg(GFF, 0, "* output GF after fixing (m) [0=don't]");
2119  setcfg(FIN, 0, "* output WL/GF after final check [0=don't]");
2120 
2121  // options to fix input SatPass, and/or generate editcmds
2122  setcfg(doFix, 0, "apply fixes to input L1 and L2 SatPass arrays");
2123  setcfg(doCmds, 0, "generate editing commands");
2124  setcfg(doRINEX3, 1, "* editing commands use L1C L2W instead of L1 L2");
2125 
2126  // when rejecting data, set SatPass UserFlag using this value
2127  setcfg(UserFlag, 0, "* call SatPass::setUserFlag(value) for rejects");
2128  }
2129  catch (Exception& e)
2130  {
2131  GNSSTK_RETHROW(e);
2132  }
2133  } // end void gdc::init()
2134 
2135 //------------------------------------------------------------------------------------
2136 #undef setcfg
2137 #undef cfg
2138 
2139 } // end namespace gnsstk
gnsstk::WindowFilter::setWidth
void setWidth(int w)
Definition: WindowFilter.hpp:367
gnsstk::StatsFilterBase::N
virtual unsigned int N() const =0
return the sample size
gnsstk::FilterHit::index
unsigned int index
index in the data array(s) at which this event occurs
Definition: StatsFilterHit.hpp:91
gdc.hpp
gnsstk::Arc::Arcinfo::n
int n
number of points in stats(ave,sig) - may be > npts in Arc
Definition: gdc.hpp:138
gnsstk::DiscontinuityCorrector
int DiscontinuityCorrector(SatPass &SP, GDCconfiguration &config, std::vector< std::string > &EditCmds, std::string &retMsg, int GLOn=-99)
Definition: DiscCorr.cpp:698
gnsstk::SatPass::data
double & data(unsigned int i, const std::string &type)
Definition: SatPass.cpp:882
gnsstk::FirstDiffFilter::getResults
std::vector< FilterHit< T > > getResults()
return a copy of the filter hit results
Definition: FirstDiffFilter.hpp:156
gnsstk::StatsFilterBase::StdDev
virtual T StdDev() const =0
return computed standard deviation, in 2-sample stats this is SigmaYX()
gnsstk::WindowFilter::dump
void dump(std::ostream &s, std::string tag=std::string())
Definition: WindowFilter.hpp:1341
L1
gnsstk::Matrix< double > L1
Definition: Matrix_LUDecomp_T.cpp:46
gnsstk::FirstDiffFilter::dump
void dump(std::ostream &s, const std::string &tag=std::string())
Definition: FirstDiffFilter.hpp:682
gnsstk::WindowFilter::setw
void setw(int w)
Definition: WindowFilter.hpp:403
stl_helpers.hpp
gnsstk::Arc::index
int index
index in data arrays of beginning ("break")
Definition: gdc.hpp:150
gnsstk::SatPass::getFirstTime
Epoch getFirstTime() const
Definition: SatPass.cpp:1012
gnsstk::Arc::GFinfo
Arcinfo GFinfo
Definition: gdc.hpp:155
logstream.hpp
gnsstk::WindowFilter::getStats
void getStats(FilterHit< T > &fe, bool noedge=true)
Definition: WindowFilter.hpp:1445
gnsstk::WindowFilter::setDumpAnalMsg
void setDumpAnalMsg(bool b)
Definition: WindowFilter.hpp:394
gnsstk::SatPass
Definition: SatPass.hpp:71
gnsstk::FirstDiffFilter::filter
int filter(const size_t index=0, int npts=-1)
Definition: FirstDiffFilter.hpp:266
gnsstk::WindowFilter::filter
int filter(const size_t index=0, int npts=-1)
Definition: WindowFilter.hpp:506
gnsstk::Arc::Arcinfo::Nslip
int Nslip
Definition: gdc.hpp:134
gnsstk::SatPass::getOutputFormat
std::string getOutputFormat()
Definition: SatPass.hpp:417
gnsstk::Arc::npts
unsigned npts
number of points in the segment -> last index
Definition: gdc.hpp:151
gnsstk::getWavelength
double getWavelength(SatelliteSystem sys, int rinexBandNum, int gloChan=0) noexcept
Definition: FreqConsts.hpp:191
GNSSconstants.hpp
gnsstk::WARNING
@ WARNING
Definition: logstream.hpp:57
LOGstrm
#define LOGstrm
Definition: logstream.hpp:322
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
P1
gnsstk::Matrix< double > P1
Definition: Matrix_LUDecomp_T.cpp:49
gnsstk::TwoSampleStatsFilter
A StatsFilter class for two-sample statistics that inherits StatsFilterBase.
Definition: WindowFilter.hpp:231
gnsstk::getAlpha
double getAlpha(SatelliteSystem sys, int na, int nb) noexcept
Definition: FreqConsts.hpp:286
gnsstk::Arc::ngood
unsigned ngood
number of good data points in the segment
Definition: gdc.hpp:152
gnsstk::Exception
Definition: Exception.hpp:151
gnsstk::StringUtils::stripTrailing
std::string & stripTrailing(std::string &s, const std::string &aString, std::string::size_type num=std::string::npos)
Definition: StringUtils.hpp:1453
gnsstk::SatPass::getObsTypes
std::vector< std::string > getObsTypes()
Definition: SatPass.hpp:449
P2
gnsstk::Matrix< double > P2
Definition: Matrix_LUDecomp_T.cpp:49
gnsstk::Arc::mark
unsigned mark
bitmap identifying actions that created or modified Arc
Definition: gdc.hpp:149
gnsstk::FirstDiffFilter
Definition: FirstDiffFilter.hpp:111
debug
#define debug
Definition: Rinex3ClockHeader.cpp:51
gnsstk::Arc::Arcinfo::sig
double sig
Definition: gdc.hpp:140
gnsstk::StatsFilterBase::Add
virtual void Add(const T &x, const T &y)=0
Add data to the statistics; in 1-sample stats the x is ignored.
gnsstk::WindowFilter
Definition: WindowFilter.hpp:302
gnsstk::OneSampleStatsFilter
A StatsFilter class for one-sample statistics that inherits StatsFilterBase.
Definition: WindowFilter.hpp:177
gnsstk::ERROR
@ ERROR
Definition: logstream.hpp:57
gnsstk::FirstDiffFilter::setprecision
void setprecision(int p)
Definition: FirstDiffFilter.hpp:152
gnsstk::FirstDiffFilter::results
std::vector< FilterHit< T > > results
vector of FilterHit, generated by analyze(), also for use in dump()
Definition: FirstDiffFilter.hpp:250
gnsstk::SatPass::getDT
double getDT() const
Definition: SatPass.hpp:515
gnsstk::min
T min(const SparseMatrix< T > &SM)
Maximum element - return 0 if empty.
Definition: SparseMatrix.hpp:858
gnsstk::FilterHit::npts
unsigned int npts
number data points in segment (= a delta index)
Definition: StatsFilterHit.hpp:93
gnsstk::GalHealthStatus::OK
@ OK
Signal OK.
wl1
double wl1
Definition: DiscCorr.cpp:656
gnsstk::VERBOSE
@ VERBOSE
Definition: logstream.hpp:57
gnsstk::FirstDiffFilter::setLimit
void setLimit(T val)
get and set
Definition: FirstDiffFilter.hpp:148
gnsstk::Arc
Definition: gdc.hpp:82
gnsstk::beta
double beta(double x, double y)
Definition: SpecialFuncs.cpp:204
gnsstk::StringUtils::asDouble
double asDouble(const std::string &s)
Definition: StringUtils.hpp:705
wl2
double wl2
Definition: DiscCorr.cpp:656
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
example4.pos
pos
Definition: example4.py:125
gnsstk::SatPass::getGLOchannel
bool getGLOchannel(int &n, std::string &msg)
Definition: SatPass.cpp:309
gnsstk::WindowFilter::setprecision
void setprecision(int p)
Definition: WindowFilter.hpp:404
gnsstk::WindowFilter::setTwoSample
void setTwoSample(bool b)
Definition: WindowFilter.hpp:369
gnsstk::SatPass::spdvector
std::vector< SatPassData > spdvector
ALL data in the pass, stored in SatPassData objects, in time order.
Definition: SatPass.hpp:185
gnsstk::getBeta
double getBeta(SatelliteSystem sys, int na, int nb) noexcept
Definition: FreqConsts.hpp:271
gnsstk::SatPass::setUserFlag
void setUserFlag(unsigned int i, unsigned int inflag)
Definition: SatPass.cpp:964
LOG
#define LOG(level)
define the macro that is used to write to the log stream
Definition: logstream.hpp:315
gnsstk::FirstDiffFilter::analyze
int analyze()
Definition: FirstDiffFilter.hpp:334
gnsstk::printTime
std::string printTime(const CommonTime &t, const std::string &fmt)
Definition: TimeString.cpp:64
gnsstk::RinexSatID
Definition: RinexSatID.hpp:63
std
Definition: Angle.hpp:142
cfg
#define cfg(a)
Definition: gdc.cpp:97
gnsstk::INFO
@ INFO
Definition: logstream.hpp:57
gnsstk::Epoch
Definition: Epoch.hpp:123
gnsstk::StringUtils::leftJustify
std::string & leftJustify(std::string &s, const std::string::size_type length, const char pad=' ')
Definition: StringUtils.hpp:1582
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::DEBUG
@ DEBUG
Definition: logstream.hpp:57
gnsstk::SatPass::time
Epoch time(unsigned int i) const
Definition: SatPass.cpp:1097
gnsstk::Arc::WLinfo
Arcinfo WLinfo
Arcinfo for each datatype.
Definition: gdc.hpp:155
gnsstk::WindowFilter::getResults
std::vector< FilterHit< T > > getResults()
get results vector of FilterHit
Definition: WindowFilter.hpp:407
gnsstk::WindowFilter::setMinStep
void setMinStep(T val)
Definition: WindowFilter.hpp:380
gnsstk::FilterHit
Definition: StatsFilterHit.hpp:68
gnsstk::StatsFilterBase
Definition: WindowFilter.hpp:129
setcfg
#define setcfg(a, b, c)
Definition: gdc.cpp:2063
GLOn
int GLOn
Definition: DiscCorr.cpp:655
L2
gnsstk::Matrix< double > L2
Definition: Matrix_LUDecomp_T.cpp:46
gnsstk::WindowFilter::results
std::vector< FilterHit< T > > results
Definition: WindowFilter.hpp:496
gnsstk::SatPass::getSat
RinexSatID getSat() const
Definition: SatPass.hpp:509
gnsstk::SatPass::setFlag
void setFlag(unsigned int i, unsigned short flag)
Definition: SatPass.cpp:942
gnsstk::StatsFilterBase::Average
virtual T Average() const =0
return the average; in 2-sample stats this is AverageY()
gnsstk::SatPass::size
unsigned int size() const
Definition: SatPass.hpp:527
gnsstk::Arc::Arcinfo::ave
double ave
average value of the data in the Arc (data units)
Definition: gdc.hpp:139
gnsstk::vectorindex
int vectorindex(const std::vector< T > &vec, const T &value)
Definition: stl_helpers.hpp:123
gnsstk::FirstDiffFilter::setw
void setw(int w)
get and set for dump
Definition: FirstDiffFilter.hpp:151
gnsstk::FirstDiffFilter::getStats
void getStats(FilterHit< T > &fe)
Definition: FirstDiffFilter.hpp:730
gnsstk::WindowFilter::analyze
int analyze()
Definition: WindowFilter.hpp:945

gnsstk
Author(s):
autogenerated on Wed Oct 25 2023 02:40:39