DiscCorr.cpp

Go to the documentation of this file.
00001 #pragma ident "$Id: DiscCorr.cpp 2741 2011-06-22 16:37:02Z nwu $"
00002 
00003 //============================================================================
00004 //
00005 //  This file is part of GPSTk, the GPS Toolkit.
00006 //
00007 //  The GPSTk is free software; you can redistribute it and/or modify
00008 //  it under the terms of the GNU Lesser General Public License as published
00009 //  by the Free Software Foundation; either version 2.1 of the License, or
00010 //  any later version.
00011 //
00012 //  The GPSTk is distributed in the hope that it will be useful,
00013 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 //  GNU Lesser General Public License for more details.
00016 //
00017 //  You should have received a copy of the GNU Lesser General Public
00018 //  License along with GPSTk; if not, write to the Free Software Foundation,
00019 //  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020 //
00021 //  Copyright 2004, The University of Texas at Austin
00022 //
00023 //============================================================================
00024 
00034 //------------------------------------------------------------------------------------
00035 #include <string>
00036 #include <iostream>
00037 #include <sstream>
00038 #include <vector>
00039 #include <deque>
00040 #include <list>
00041 #include <algorithm>
00042 
00043 #include "StringUtils.hpp"
00044 #include "Stats.hpp"
00045 #include "PolyFit.hpp"
00046 #include "icd_200_constants.hpp"    // PI,C_GPS_M,OSC_FREQ,L1_MULT,L2_MULT
00047 #include "RobustStats.hpp"
00048 
00049 #include "DiscCorr.hpp"
00050 
00051 using namespace std;
00052 using namespace gpstk;
00053 using namespace StringUtils;
00054 
00055 //------------------------------------------------------------------------------------
00056 // class GDCconfiguration
00057 //------------------------------------------------------------------------------------
00058 // class GDCconfiguration: string giving version of gpstk Discontinuity Corrector
00059 string GDCconfiguration::GDCVersion = string("5.3 7/14/2008");
00060 
00061 // class GDCconfiguration: member functions
00062 //------------------------------------------------------------------------------------
00063 // Set a parameter in the configuration; the input string 'cmd' is of the form
00064 // '[--DC]<id><s><value>' : separator s is one of ':=,' and leading --DC is optional.
00065 void GDCconfiguration::setParameter(string cmd) throw(Exception)
00066 {
00067 try {
00068    if(cmd.empty()) return;
00069       // remove leading --DC
00070    while(cmd[0] == '-') cmd.erase(0,1);
00071    if(cmd.substr(0,2) == "DC") cmd.erase(0,2);
00072 
00073    string label, value;
00074    string::size_type pos=cmd.find_first_of(",=:");
00075    if(pos == string::npos) {
00076       label = cmd;
00077    }
00078    else {
00079       label = cmd.substr(0,pos);
00080       value = cmd;
00081       value.erase(0,pos+1);
00082    }
00083 
00084    setParameter(label, asDouble(value));
00085 }
00086 catch(Exception& e) { GPSTK_RETHROW(e); }
00087 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00088 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00089 }
00090 
00091 //------------------------------------------------------------------------------------
00092 // Set a parameter in the configuration using the label and the value,
00093 // for booleans use (T,F)=(non-zero,zero).
00094 void GDCconfiguration::setParameter(string label, double value) throw(Exception)
00095 {
00096 try {
00097    if(CFG.find(label) == CFG.end())
00098       ; // throw
00099    else {
00100       // log is not defined yet
00101       if(CFG["Debug"] > 0) *(p_oflog) << "GDCconfiguration::setParameter sets "
00102          << label << " to " << value << endl;
00103       CFG[label] = value;
00104    }
00105 }
00106 catch(Exception& e) { GPSTK_RETHROW(e); }
00107 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00108 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00109 }
00110 
00111 //------------------------------------------------------------------------------------
00112 // Print help page, including descriptions and current values of all
00113 // the parameters, to the ostream.
00114 void GDCconfiguration::DisplayParameterUsage(ostream& os, bool advanced)
00115    throw(Exception)
00116 {
00117 try {
00118    os << "GPSTk Discontinuity Corrector (GDC) v." << GDCVersion
00119       << " configuration:"
00120       //<< "  [ pass setParameter() a string '<label><sep><value>';"
00121       //<< " <sep> is one of ,=: ]"
00122       << endl;
00123 
00124    map<string,double>::const_iterator it;
00125    for(it=CFG.begin(); it != CFG.end(); it++) {
00126       if(CFGdescription[it->first][0] == '*')      // advanced options
00127          continue;
00128       ostringstream stst;
00129       stst << it->first                            // label
00130          << "=" << it->second;                     // value
00131       os << " " << leftJustify(stst.str(),18)
00132          << " : " << CFGdescription[it->first]     // description
00133          << endl;
00134    }
00135    if(advanced) {
00136       os << "   Advanced options:" << endl;
00137    for(it=CFG.begin(); it != CFG.end(); it++) {
00138       if(CFGdescription[it->first][0] != '*')      // ordinary options
00139          continue;
00140       ostringstream stst;
00141       stst << it->first                            // label
00142          << "=" << it->second;                     // value
00143       os << " " << leftJustify(stst.str(),25)
00144          << " : " << CFGdescription[it->first].substr(2)  // description
00145          << endl;
00146    }
00147    }
00148 }
00149 catch(Exception& e) { GPSTK_RETHROW(e); }
00150 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00151 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00152 }
00153 
00154 
00155 //------------------------------------------------------------------------------------
00156 #define setcfg(a,b,c) { CFG[#a]=b; CFGdescription[#a]=c; }
00157 // initialize with default values
00158 void GDCconfiguration::initialize(void)
00159 {
00160 try {
00161    p_oflog = &cout;
00162 
00163    // use cfg(DT) NOT dt -  dt is part of SatPass...
00164    setcfg(DT, -1, "nominal timestep of data (seconds) [required - no default!]");
00165    setcfg(Debug, 0, "level of diagnostic output to log, from none(0) to extreme(7)");
00166    setcfg(useCA, 0, "use C/A code pseudorange (C1) ()");
00167    setcfg(MaxGap, 180, "maximum allowed time gap within a segment (seconds)");
00168    setcfg(MinPts, 13, "minimum number of good points in phase segment ()");
00169    setcfg(WLSigma, 1.5, "expected WL sigma (WL cycle) [NB = ~0.83*p-range noise(m)]");
00170    setcfg(GFVariation, 16,                    // about 300 5.4-cm wavelengths
00171       "expected maximum variation in GF phase in time DT (meters)");
00172    // output
00173    setcfg(OutputGPSTime, 0,
00174       "if 0: Y,M,D,H,M,S  else: W,SoW (GPS) in editing commands");
00175    setcfg(OutputDeletes, 1,
00176       "if non-zero, include delete commands in the output cmd list");
00177 
00178    // -------------------------------------------------------------------------
00179    // advanced options - ordinary user will most likely NOT change
00180    setcfg(RawBiasLimit, 100, "* change in raw R-Ph that triggers bias reset (m)");
00181    // WL editing
00182    setcfg(WLNSigmaDelete, 2, "* delete segments with sig(WL) > this * WLSigma ()");
00183    // 040809 increased 10->20: ~colored variation was so large test >~ limit
00184    // 050109 change 20->10+50/dt; implement in preprocess()
00185    setcfg(WLWindowWidth, 50, "* sliding window width for WL slip detection = 10+this/dt) (points)");
00186    setcfg(WLNWindows, 2.5,
00187       "* minimum segment size for WL small slip search (WLWindowWidth)");
00188    setcfg(WLobviousLimit, 3,
00189       "* minimum delta(WL) that produces an obvious slip (WLSigma)");
00190    setcfg(WLNSigmaStrip, 3.5, "* delete points with WL > this * computed sigma ()");
00191    setcfg(WLNptsOutlierStats, 200,
00192       "* maximum segment size to use robust outlier detection (pts)");
00193    setcfg(WLRobustWeightLimit, 0.35,
00194       "* minimum good weight in robust outlier detection (0<wt<=1)");
00195    // WL small slips
00196    setcfg(WLSlipEdge, 3,
00197       "* minimum separating WL slips and end of segment, else edit (pts)");
00198    //050109 .67->1.0
00199    setcfg(WLSlipSize, 1.0, "* minimum WL slip size (WL wavelengths)");
00200    setcfg(WLSlipExcess, 0.1,
00201       "* minimum amount WL slip must exceed noise (WL wavelengths)");
00202    //050109 1.2->2.5
00203    setcfg(WLSlipSeparation, 2.5, "* minimum excess/noise ratio of WL slip ()");
00204    // GF small slips
00205    setcfg(GFSlipWidth, 5,
00206       "* minimum segment length for GF small slip detection (pts)");
00207    setcfg(GFSlipEdge, 3,
00208       "* minimum separating GF slips and end of segment, else edit (pts)");
00209    setcfg(GFobviousLimit, 1,
00210       "* minimum delta(GF) that produces an obvious slip (GFVariation)");
00211    setcfg(GFSlipOutlier, 5, "* minimum GF outlier magnitude/noise ratio ()");
00212    setcfg(GFSlipSize, 0.8, "* minimum GF slip size (5.4cm wavelengths)");
00213    setcfg(GFSlipStepToNoise, 2, "* maximum GF slip step/noise ratio ()");
00214    setcfg(GFSlipToStep, 3, "* minimum GF slip magnitude/step ratio ()");
00215    setcfg(GFSlipToNoise, 3, "* minimum GF slip magnitude/noise ratio ()");
00216    // GF fix
00217    setcfg(GFFixNpts, 15,
00218       "* maximum number of points on each side to fix GF slips ()");
00219    setcfg(GFFixDegree, 3, "* degree of polynomial used to fix GF slips ()");
00220    setcfg(GFFixMaxRMS, 100,
00221       "* limit on RMS fit residuals to fix GF slips, else delete (5.4cm)");
00222 
00223 }
00224 catch(Exception& e) { GPSTK_RETHROW(e); }
00225 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00226 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00227 }
00228 
00229 //------------------------------------------------------------------------------------
00230 //------------------------------------------------------------------------------------
00231 // class Segment - used internally only.
00232 // An object to hold information about segments = periods of continuous phase.
00233 // Keep a linked list of these objects, subdivide whenever a discontinuity is
00234 // detected, and join whenever one is fixed.
00235 class Segment {
00236 public:
00237       // member data
00238    int nbeg,nend;          // array indexes of the first and last good points
00239                            // always maintain these so they point to good data
00240    int npts;               // number of good points in this Segment
00241    int nseg;               // segment number - used for data dumps only
00242    double bias1;           // bias subtracted from WLbias for WLStats - only
00243    Stats<double> WLStats;  // includes N,min,max,ave,sig
00244    double bias2;           // bias subtracted from GFP for polynomial fit - only
00245    PolyFit<double> PF;     // for fit to GF range
00246    double RMSROF;          // RMS residual of fit of polynomial (PF) to GFR
00247    bool WLsweep;           // WLstatSweep(this) was called, used by detectWLsmallSlips
00248 
00249       // member functions
00250    Segment(void) : nbeg(0),nend(0),npts(0),nseg(0),bias1(0.0),bias2(0.0)
00251       { WLStats.Reset(); WLsweep=false; }
00252    Segment(const Segment& s)
00253       { *this = s; }
00254    ~Segment(void) { }
00255    Segment& operator=(const Segment& s) {
00256       if(this==&s) return *this;
00257       nbeg = s.nbeg; nend = s.nend; npts = s.npts; nseg = s.nseg;
00258       bias1 = s.bias1; bias2 = s.bias2;
00259       WLStats = s.WLStats; PF = s.PF; RMSROF = s.RMSROF; WLsweep = s.WLsweep;
00260       return *this;
00261    }
00262 }; // end class Segment
00263 
00264 //------------------------------------------------------------------------------------
00265 //------------------------------------------------------------------------------------
00266 // class Slip - used internally only
00267 class Slip {
00268 public:
00269    int index;           // index in arrays where this slip occurs
00270    long NWL,N1;         // slip fixes for WL (N1-N2) and GF (=N1)
00271    string msg;          // string to be output after '#' on edit cmds
00272    explicit Slip(int in) : index(in),NWL(0),N1(0) { }
00273    bool operator<(const Slip &rhs) const { return index < rhs.index; }
00274 }; // end class Slip
00275 
00276 //------------------------------------------------------------------------------------
00277 //------------------------------------------------------------------------------------
00278 // class GDCPass - inherit SatPass and GDConfiguration - use internally only
00279 // this object is used to implement the DC algorithm.
00280 class GDCPass : public SatPass, public GDCconfiguration
00281 {
00282 public:
00283    static const unsigned short DETECT;           // both WL and GF
00284    static const unsigned short FIX;              // both WL and GF
00285    static const unsigned short WLDETECT;
00286    static const unsigned short WLFIX;
00287    static const unsigned short GFDETECT;
00288    static const unsigned short GFFIX;
00289 
00290    explicit GDCPass(SatPass& sp, const GDCconfiguration& gdc);
00291 
00292    //~GDCPass(void) { };
00293 
00295    int preprocess(void) throw(Exception);
00296 
00300    int linearCombinations(void) throw(Exception);
00301 
00303    int detectWLslips(void) throw(Exception);
00304 
00308    int detectObviousSlips(string which) throw(Exception);
00309 
00313    int firstDifferences(string which) throw(Exception);
00314 
00317    void WLcomputeStats(list<Segment>::iterator& it) throw(Exception);
00318 
00321    void WLsigmaStrip(list<Segment>::iterator& it) throw(Exception);
00322 
00326    int WLstatSweep(list<Segment>::iterator& it, int width) throw(Exception);
00327 
00332    int detectWLsmallSlips(void) throw(Exception);
00333 
00343    bool foundWLsmallSlip(list<Segment>::iterator& it, int i) throw(Exception);
00344 
00347    int fixAllSlips(string which) throw(Exception);
00348 
00351    void fixOneSlip(list<Segment>::iterator& kt, string which) throw(Exception);
00352 
00354    void WLslipFix(list<Segment>::iterator& left,
00355                   list<Segment>::iterator& right)
00356       throw(Exception);
00357 
00361    int prepareGFdata(void) throw(Exception);
00362 
00364    int detectGFslips(void) throw(Exception);
00365 
00368    int GFphaseResiduals(list<Segment>::iterator& it) throw(Exception);
00369 
00373    int detectGFsmallSlips(void) throw(Exception);
00374 
00381    bool foundGFoutlier(int i,int inew,Stats<double>& pastSt,Stats<double>& futureSt)
00382       throw(Exception);
00383 
00394    bool foundGFsmallSlip(int i, int nseg, int iend, int ibeg,
00395       deque<int>& pastIn, deque<int>& futureIn,
00396       Stats<double>& pastSt, Stats<double>& futureSt)
00397       throw(Exception);
00398 
00400    void GFslipFix(list<Segment>::iterator& left,
00401                   list<Segment>::iterator& right)
00402       throw(Exception);
00403 
00407    long EstimateGFslipFix(list<Segment>::iterator& left,
00408                           list<Segment>::iterator& right,
00409                           int nb, int ne, long n1)
00410       throw(Exception);
00411 
00413    int WLconsistencyCheck(void) throw(Exception);
00414 
00417    string finish(int iret, SatPass& svp, vector<string>& editCmds) throw(Exception);
00418 
00422    list<Segment>::iterator createSegment(list<Segment>::iterator sit,
00423                                          int ibeg,
00424                                          string msg=string())
00425       throw(Exception);
00426 
00432    string dumpSegments(string msg, int level=2, bool extra=false)
00433       throw(Exception);
00434 
00437    void deleteSegment(list<Segment>::iterator& it, string msg=string())
00438       throw(Exception);
00439 
00440 private:
00441 
00444    double cfg_func(string a) throw(Exception)
00445    {
00446       if(CFGdescription[a] == string()) {
00447          Exception e("cfg(UNKNOWN LABEL) : " + a);
00448          GPSTK_THROW(e);
00449       }
00450       return CFG[a];
00451    }
00452 
00455    list<Segment> SegList;
00456 
00458    list<Slip> SlipList;
00459 
00461    //vector<double> A1,A2;
00462 
00464    Stats<double> WLPassStats;
00465 
00467    Stats<double> GFPassStats;
00468 
00470    PolyFit<double> GFPassFit;
00471 
00473    map<string,int> learn;
00474 
00475 }; // end class GDCPass
00476 
00477 //------------------------------------------------------------------------------------
00478 // local data
00479 //------------------------------------------------------------------------------------
00480 // these are used only to associate a unique number in the log file with each pass
00481 int GDCUnique=0;     // unique number for each call
00482 int GDCUniqueFix;    // unique for each (WL,GF) fix
00483 
00484 // conveniences only...
00485 #define log *(p_oflog)
00486 #define cfg(a) cfg_func(#a)
00487 // gcc doesn't like const enum...
00488 enum obstypeenum {  L1=0, L2=1, P1=2, P2=3, A1=4, A2=5 };   // P1 will <=> C1 or P1
00489 vector<string> DCobstypes;       // above are indexes into both data and this vector
00490 
00491 // constants used in linear combinations
00492 const double CFF=C_GPS_M/OSC_FREQ;
00493 const double F1=L1_MULT;   // 154.0;
00494 const double F2=L2_MULT;   // 120.0;
00495    // wavelengths
00496 const double wl1=CFF/F1;                        // 19.0cm
00497 const double wl2=CFF/F2;                        // 24.4cm
00498 const double wlwl=CFF/(F1-F2);                  // 86.2cm, the widelane wavelength
00499 const double wl21=CFF*(1.0/F2 - 1.0/F1);        // 5.4cm, the 'GF' wavelength
00500    // for widelane R & Ph
00501 const double wl1r=F1/(F1+F2);
00502 const double wl2r=F2/(F1+F2);
00503 const double wl1p=wl1*F1/(F1-F2);
00504 const double wl2p=-wl2*F2/(F1-F2);
00505    // for geometry-free R and Ph
00506 const double gf1r=-1;
00507 const double gf2r=1;
00508 const double gf1p=wl1;
00509 const double gf2p=-wl2;
00510 
00511 //------------------------------------------------------------------------------------
00512 // Return values (used by all routines within this module):
00513 const int BadInput=-5;
00514 const int NoData=-4;
00515 const int FatalProblem=-3;
00516 //const int PrematureEnd=-2;
00517 const int Singular=-1;
00518 const int ReturnOK=0;
00519 // NoData:
00520 //    preprocess
00521 //    fixAllSlips       segment list is empty
00522 // FatalProblem:
00523 //    preprocess        DT is not set
00524 //    firstDifferences  A1.size is wrong
00525 // Singular:
00526 //    prepareGFdata     polynomial fit to GF range failed
00527 //    no - delete segment instead GFphaseResiduals  polynomial fit to GF range failed
00528 //
00529 // -----------------------------------------------------------------------------------
00530 // Roadmap:
00531 // DiscontinuityCorrector(SatPass& svp,GDCconfiguration& gdc,vector<string>& editCmds)
00532 //                             Dump   Arrays with units (on output)
00533 //                             tag    L1              L2              P1              P2              A1              A2
00534 // preprocess()                BEF    L1(c)           L2(c)           P1(m)           P2(m)           P-L1(m)         P-L2(m)
00535 // linearCombinations()        LCD    gfp+gfr         gfp             wlbias          -gfr
00536 // detectWLslips()
00537 //    detectObviousSlips(WL)
00538 //       firstDifferences(WL)  DWL                                                                     dP1=dWLb        0
00539 //    WLcomputeStats
00540 //    WLsigmaStrip                                                                                    wlbias          0
00541 //    WLstatSweep             (WLS is stats)                                                          test            limit
00542 //                                                                                                    =dave(fut-past) =sqrt(sum var)
00543 //    detectWLsmallSlips       WLD
00544 // fixAllSlips("WL")
00545 //    fixOneSlip("WL")
00546 //       WLslipFix             WLF
00547 // prepareGFdata()                    gfp+gfr(resid)  gfp(wl21)                       -gfr(wl21)
00548 // detectGFslips()
00549 //    detectObviousSlips(GF)
00550 //       firstDifferences(GF)  DGF                                                                     d(gfp+gfr res)   d(gfp)
00551 //    GFphaseResiduals                                                                                1stdiff(gfp-(fit gfr))
00552 //    detectGFsmallSlips       GFD
00553 // WLconsistencyCheck()
00554 // fixAllSlips("GF")
00555 //    fixOneSlip("GF")
00556 //       GFslipFix             GFF
00557 // finish (fix slips, recomp)  AFT    L1(c)           L2(c)           P1(m)           P2(m)           P-L1(m)         P-L2(m)
00558 //
00559 //------------------------------------------------------------------------------------
00560 // Constants used to mark slips, etc. using the SatPass flag:
00561 // These are from SatPass.cpp
00562 //const unsigned short SatPass::BAD = 0; // used by caller and DC to mark bad data
00563 //const unsigned short SatPass::OK  = 1; // good data, no discontinuity
00564 //const unsigned short SatPass::LL1 = 2; // good data, discontinuity on L1 only
00565 //const unsigned short SatPass::LL2 = 4; // good data, discontinuity on L2 only
00566 //const unsigned short SatPass::LL3 = 6; // good data, discontinuity on L1 and L2
00567 const unsigned short GDCPass::WLDETECT =   2;
00568 const unsigned short GDCPass::GFDETECT =   4;
00569 const unsigned short GDCPass::DETECT   =   6;  // = WLDETECT | GFDETECT
00570 const unsigned short GDCPass::WLFIX    =   8;
00571 const unsigned short GDCPass::GFFIX    =  16;
00572 const unsigned short GDCPass::FIX      =  24;  // = WLFIX | GFFIX
00573 
00574 // notes on the use of these flags:
00575 // if(flag & DETECT) is true for EITHER WL or GF or both
00576 // if(flag & FIX)  is true for EITHER WL or GF or both
00577 // if((flag & WLDETECT) && (flag & GFDETECT)) is true only for both WL and GF
00578 
00579 // NB typical slip will have flag = DETECT+OK+FIX = 31
00580 //    typical unfixed slip   flag = DETECT+OK     =  7
00581 
00582 // BAD is used either as flag == BAD (for bad data) or flag != BAD (for good data),
00583 // thus there are two gotcha's
00584 //   - if a point is marked, but is later set BAD, that info is lost
00585 //   - if a BAD point is marked, it becomes 'good'
00586 // To avoid this we have to use OK rather than BAD:
00587 // either !(flag & OK) or (flag ^ OK) for bad data, and (flag & OK) for good data
00588 
00589 //------------------------------------------------------------------------------------
00590 // The discontinuity corrector function
00591 //------------------------------------------------------------------------------------
00592 // yes you need the gpstk::
00593 int gpstk::DiscontinuityCorrector(SatPass& svp,
00594                                   GDCconfiguration& gdc,
00595                                   vector<string>& editCmds,
00596                                   string& retMessage)
00597    throw(Exception)
00598 {
00599 try {
00600    int i,j,iret;
00601    GDCUnique++;
00602 
00603    // --------------------------------------------------------------------------------
00604    // require obstypes L1,L2,C1/P1,P2, and add two auxiliary arrays
00605    DCobstypes.clear();
00606    DCobstypes.push_back("L1");
00607    DCobstypes.push_back("L2");
00608    DCobstypes.push_back((int(gdc.getParameter("useCA"))) == 0 ? "P1" : "C1");
00609    DCobstypes.push_back("P2");
00610    DCobstypes.push_back("A1");
00611    DCobstypes.push_back("A2");
00612 
00613    // --------------------------------------------------------------------------------
00614    // test input for (a) some data and (b) the required obs types L1,L2,C1/P1,P2
00615    vector<double> newdata(6);
00616    try {
00617       newdata[L1] = svp.data(0,DCobstypes[L1]);
00618       newdata[L2] = svp.data(0,DCobstypes[L2]);
00619       newdata[P1] = svp.data(0,DCobstypes[P1]);
00620       newdata[P2] = svp.data(0,DCobstypes[P2]);
00621    }
00622    catch (Exception& e) { // if obs type is not found in input
00623       return BadInput;
00624    }
00625 
00626    // --------------------------------------------------------------------------------
00627    // create a SatPass using DCobstypes, and fill from input
00628    SatPass nsvp(svp.getSat(), svp.getDT(), DCobstypes);
00629 
00630    // fill the new SatPass with the input data
00631    nsvp.status() = svp.status();
00632    vector<unsigned short> lli(6),ssi(6);
00633    for(i=0; i<svp.size(); i++) {
00634       for(j=0; j<6; j++) {
00635          newdata[j] = j < 4 ? svp.data(i,DCobstypes[j]) : 0.0;
00636          lli[j] = j < 4 ? svp.LLI(i,DCobstypes[j]) : 0;
00637          ssi[j] = j < 4 ? svp.SSI(i,DCobstypes[j]) : 0;
00638       }
00639       // return value must be 0
00640       nsvp.addData(svp.time(i), DCobstypes, newdata, lli, ssi, svp.getFlag(i));
00641    }
00642 
00643    // --------------------------------------------------------------------------------
00644    // create a GDCPass from the input SatPass (modified) and GDC configuration
00645    GDCPass gp(nsvp,gdc);
00646 
00647    // --------------------------------------------------------------------------------
00648    // implement the DC algorithm using the GDCPass
00649    // NB search for 'change the arrays' for places where arrays are re-defined
00650    // NB search for 'change the data' for places where the data is modified (! biases)
00651    // NB search for 'change the bias' for places where the bias is changed
00652    for(;;) {      // a convenience...
00653       // preparation
00654       if( (iret = gp.preprocess() )) break;
00655       if( (iret = gp.linearCombinations() )) break;
00656 
00657       // WL
00658       if( (iret = gp.detectWLslips() )) break;
00659       if( (iret = gp.fixAllSlips("WL") )) break;
00660 
00661       // GF
00662       if( (iret = gp.prepareGFdata() )) break;
00663       if( (iret = gp.detectGFslips() )) break;
00664       if( (iret = gp.WLconsistencyCheck() )) break;
00665       if( (iret = gp.fixAllSlips("GF") )) break;
00666 
00667       break;      // mandatory
00668    }
00669 
00670    // --------------------------------------------------------------------------------
00671    // generate editing commands for deleted (flagged) data and slips,
00672    // use editing command (slips and deletes) to modify the original SatPass data
00673    // and print ending summary
00674    retMessage = gp.finish(iret, svp, editCmds);
00675 
00676    return iret;
00677 }
00678 catch(Exception& e) { GPSTK_RETHROW(e); }
00679 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00680 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00681 }
00682 
00683 //------------------------------------------------------------------------------------
00684 //------------------------------------------------------------------------------------
00685 // class GDCPass member functions
00686 //------------------------------------------------------------------------------------
00687 GDCPass::GDCPass(SatPass& sp, const GDCconfiguration& gdc)
00688       : SatPass(sp.getSat(), sp.getDT(), sp.getObsTypes())
00689 {
00690    int i,j;
00691    Status = sp.status();
00692    dt = sp.getDT();
00693    sat = sp.getSat();
00694    vector<string> ot = sp.getObsTypes();
00695    for(i=0; i<ot.size(); i++) {
00696       labelForIndex[i] = ot[i];
00697       indexForLabel[labelForIndex[i]] = i;
00698    }
00699 
00700    vector<double> vdata;
00701    vector<unsigned short> lli,ssi;
00702    for(i=0; i<sp.size(); i++) {
00703       vdata.clear();
00704       lli.clear();
00705       ssi.clear();
00706       for(j=0; j<ot.size(); j++) {
00707          vdata.push_back(sp.data(i,ot[j]));
00708          lli.push_back(sp.LLI(i,ot[j]));
00709          ssi.push_back(sp.SSI(i,ot[j]));
00710       }
00711       addData(sp.time(i),ot,vdata,lli,ssi,sp.getFlag(i));
00712    }
00713 
00714    *((GDCconfiguration*)this) = gdc;
00715 
00716    learn.clear();
00717 }
00718 
00719 //------------------------------------------------------------------------------------
00720 int GDCPass::preprocess(void) throw(Exception)
00721 {
00722 try {
00723    int i,ilast,Ngood;
00724    double biasL1,biasL2,dbias;
00725    list<Segment>::iterator it;
00726 
00727    if(cfg(Debug) >= 2) {
00728       DayTime CurrentTime;
00729       //CurrentTime.setLocalTime();
00730       log << "======== Beg GPSTK Discontinuity Corrector " << GDCUnique
00731           << " ================================================" << endl;
00732       log << "GPSTK Discontinuity Corrector Ver. " << GDCVersion << " Run "
00733           << CurrentTime << endl;
00734    }
00735 
00736       // check input
00737    if(cfg(DT) <= 0) {
00738       log << "Error: data time interval is not set...Abort" << endl;
00739       return FatalProblem;
00740    }
00741 
00742    // 050109 some parameters should depend on DT
00743    CFG["WLWindowWidth"] = 10 + int(0.5+CFG["WLWindowWidth"]/CFG["DT"]);
00744    //log << "WLWindowWidth is now " << cfg(WLWindowWidth) << endl;
00745 
00746       // create the first segment
00747    SegList.clear();
00748    {
00749       Segment s;
00750       s.nseg = 1;
00751       SegList.push_back(s);
00752    }
00753    it = SegList.begin();
00754 
00755       // loop over points in the pass
00756       // editing obviously bad data and adding segments where necessary
00757    for(ilast=-1,i=0; i<size(); i++) {
00758 
00759       // ignore data the caller has marked BAD
00760       if(!(spdvector[i].flag & OK)) continue;
00761 
00762       // just in case the caller has set it to something else...
00763       spdvector[i].flag = OK;
00764 
00765          // look for obvious outliers
00766          // Don't do this - sometimes the pseudoranges get extreme values b/c the
00767          // clock is allowed to run off for long times - perfectly normal
00768       //if(spdvector[i].data[P1] < cfg(MinRange) ||
00769       //   spdvector[i].data[P1] > cfg(MaxRange) ||
00770       //   spdvector[i].data[P2] < cfg(MinRange) ||
00771       //   spdvector[i].data[P2] > cfg(MaxRange) )
00772       //{
00773       //   spdvector[i].flag = BAD;
00774       //   learn["points deleted: obvious outlier"]++;
00775       //   if(cfg(Debug) > 6)
00776       //      log << "Obvious outlier " << GDCUnique << " " << sat
00777       //         << " at " << i << " " << time(i).printf(outFormat) << endl;
00778       //   continue;
00779       //}
00780 
00781          // note first good point
00782       if(ilast == -1) ilast = it->nbeg = i;
00783 
00784          // is there a gap here? if yes, create a new segment
00785          // TD? do this here? why not allow any gap in the WL, and look for gaps
00786          // TD? only in the GF?
00787       if(cfg(DT)*(i-ilast) > cfg(MaxGap))
00788          it = createSegment(it,i,"initial gap");
00789 
00790          // count good points
00791       it->npts++;
00792       ilast = i;
00793    }
00794 
00795       // note last good point
00796    if(ilast == -1) ilast = it->nbeg;
00797    it->nend = ilast;
00798 
00799       // 'change the arrays' A1, A2 to be range minus phase for output
00800       // do the same at the end ("AFT")
00801       // loop over segments, counting the number of non-trivial ones
00802    for(Ngood=0,it=SegList.begin(); it != SegList.end(); it++) {
00803       biasL1 = biasL2 = 0.0;
00804 
00805          // loop over points in this segment
00806       for(i=it->nbeg; i<=it->nend; i++) {
00807          if(!(spdvector[i].flag & OK)) continue;
00808 
00809          dbias = fabs(spdvector[i].data[P1]-wl1*spdvector[i].data[L1]-biasL1);
00810          if(dbias > cfg(RawBiasLimit)) {
00811             if(cfg(Debug) >= 2) log << "BEFresetL1 " << GDCUnique
00812                << " " << sat << " " << time(i).printf(outFormat)
00813                << " " << fixed << setprecision(3) << biasL1
00814                << " " << spdvector[i].data[P1] - wl1 * spdvector[i].data[L1] << endl;
00815             biasL1 = spdvector[i].data[P1] - wl1 * spdvector[i].data[L1];
00816          }
00817 
00818          dbias = fabs(spdvector[i].data[P2]-wl2*spdvector[i].data[L2]-biasL2);
00819          if(dbias > cfg(RawBiasLimit)) {
00820             if(cfg(Debug) >= 2) log << "BEFresetL2 " << GDCUnique
00821                << " " << sat << " " << time(i).printf(outFormat)
00822                << " " << fixed << setprecision(3) << biasL2
00823                << " " << spdvector[i].data[P2] - wl2 * spdvector[i].data[L2] << endl;
00824             biasL2 = spdvector[i].data[P2] - wl2 * spdvector[i].data[L2];
00825          }
00826 
00827          spdvector[i].data[A1] =
00828             spdvector[i].data[P1] - wl1 * spdvector[i].data[L1] - biasL1;
00829          spdvector[i].data[A2] =
00830             spdvector[i].data[P2] - wl2 * spdvector[i].data[L2] - biasL2;
00831 
00832       }  // end loop over points in the segment
00833 
00834          // delete small segments
00835       if(it->npts < int(cfg(MinPts)))
00836          deleteSegment(it,"insufficient data in segment");
00837       else
00838          Ngood++;
00839    }
00840 
00841    if(cfg(Debug) >= 2) dumpSegments("BEF",2,true);
00842 
00843    if(Ngood == 0) return NoData;
00844    return ReturnOK;
00845 }
00846 catch(Exception& e) { GPSTK_RETHROW(e); }
00847 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00848 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00849 }
00850 
00851 //------------------------------------------------------------------------------------
00852 //------------------------------------------------------------------------------------
00853 int GDCPass::linearCombinations(void) throw(Exception)
00854 {
00855 try {
00856    int i;
00857    double wlr,wlp,wlbias,gfr,gfp;
00858    list<Segment>::iterator it;
00859 
00860    // iterate over segments
00861    for(it=SegList.begin(); it != SegList.end(); it++) {
00862       it->npts = 0;                       // re-compute npts here
00863 
00864       // loop over points in this segment
00865       for(i=it->nbeg; i<=it->nend; i++) {
00866          if(!(spdvector[i].flag & OK)) continue;
00867 
00868          // narrow lane range (m)
00869          wlr = wl1r * spdvector[i].data[P1] + wl2r * spdvector[i].data[P2];
00870          // wide lane phase (m)
00871          wlp = wl1p * spdvector[i].data[L1] + wl2p * spdvector[i].data[L2];
00872          // geometry-free range (m)
00873          gfr =        spdvector[i].data[P1] -        spdvector[i].data[P2];
00874          // geometry-free phase (m)
00875          gfp = gf1p * spdvector[i].data[L1] + gf2p * spdvector[i].data[L2];
00876          // wide lane bias (cycles)
00877          wlbias = (wlp-wlr)/wlwl;
00878 
00879          // change the bias
00880          if(it->npts == 0) {                             // first good point
00881             it->bias1 = wlbias;                          // WL bias (NWL)
00882             it->bias2 = gfp;                             // GFP bias
00883          }
00884 
00885          // change the arrays
00886          spdvector[i].data[L1] = gfp + gfr;              // only used in GF
00887          spdvector[i].data[L2] = gfp;
00888          spdvector[i].data[P1] = wlbias;
00889          spdvector[i].data[P2] = - gfr;
00890 
00891          it->npts++;
00892       }
00893    }
00894 
00895    if(cfg(Debug) >= 2) dumpSegments("LCD");
00896 
00897    return ReturnOK;
00898 }
00899 catch(Exception& e) { GPSTK_RETHROW(e); }
00900 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00901 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00902 }
00903 
00904 //------------------------------------------------------------------------------------
00905 //------------------------------------------------------------------------------------
00906 // detect slips in the wide lane bias
00907 int GDCPass::detectWLslips(void) throw(Exception)
00908 {
00909 try {
00910    int iret;
00911    list<Segment>::iterator it;
00912 
00913    // look for obvious slips. this will break one segment into many
00914    if( (iret = detectObviousSlips("WL"))) return iret;
00915 
00916    for(it=SegList.begin(); it != SegList.end(); it++) {
00917 
00918       // compute stats and delete segments that are too small
00919       WLcomputeStats(it);
00920 
00921       // sigma-strip the WL bias, and remove small segments
00922       if(it->npts > 0) WLsigmaStrip(it);
00923 
00924       // print this before deleting segments with large sigma
00925       if(cfg(Debug) >= 1 && it->npts >= int(cfg(MinPts)))
00926          log << "WLSIG " << GDCUnique << " " << sat
00927             << " " << it->nseg
00928             << " " << time(it->nbeg).printf(outFormat)
00929             << fixed << setprecision(3)
00930             << " " << it->WLStats.StdDev()
00931             << " " << it->WLStats.Average()
00932             << " " << it->WLStats.Minimum()
00933             << " " << it->WLStats.Maximum()
00934             << " " << it->npts
00935             << " " << it->nbeg << " - " << it->nend
00936             << " " << it->bias1
00937             << " " << it->bias2
00938             << endl;
00939 
00940       // delete segments if sigma is too high...
00941       if(it->WLStats.StdDev() > cfg(WLNSigmaDelete)*cfg(WLSigma))
00942          deleteSegment(it,"WL sigma too big");
00943 
00944       // if there are less than about 2.5*cfg(WLWindowWidth) good points, don't bother
00945       // to use the sliding window method to look for slips; otherwise
00946       // compute stats for each segment using the 'two-paned sliding stats window',
00947       // store results in the temporary arrays
00948       if(double(it->npts) >= cfg(WLNWindows) * int(cfg(WLWindowWidth))) {
00949          iret = WLstatSweep(it,int(cfg(WLWindowWidth)));
00950          if(iret) return iret;
00951       }
00952 
00953    }  // end loop over segments
00954 
00955    // use the temporary arrays filled by WLstatSweep to detect slips in the WL bias
00956    // recompute WLstats, and break up the segments where slips are found
00957    if( (iret = detectWLsmallSlips())) return iret;
00958 
00959    // delete all segments that are too small
00960    for(it=SegList.begin(); it != SegList.end(); it++) {
00961       if(it->npts < int(cfg(MinPts)))
00962          deleteSegment(it,"insufficient data in segment");
00963    }
00964 
00965    if(cfg(Debug) >= 4) dumpSegments("WLD");
00966 
00967    return ReturnOK;
00968 }
00969 catch(Exception& e) { GPSTK_RETHROW(e); }
00970 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
00971 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
00972 }
00973 
00974 //------------------------------------------------------------------------------------
00975 // detect obvious slips by computing the first difference (of either WL or GFP)
00976 // and looking for outliers. create new segments where there are slips
00977 // which is either 'WL' or 'GF'.
00978 int GDCPass::detectObviousSlips(string which) throw(Exception)
00979 {
00980 try {
00981       // TD determine from range noise // ~ 2*range noise/wl2
00982    const double WLobviousNwlLimit=cfg(WLobviousLimit)*cfg(WLSigma);
00983    const double GFobviousNwlLimit=cfg(GFobviousLimit)*cfg(GFVariation)/wl21;
00984    bool outlier;
00985    int i,j,iret,ibad=0,igood,nok;
00986    double limit,wlbias;
00987    list<Segment>::iterator it;
00988 
00989    // compute 1st differences of (WL bias, GFP-GFR) as 'which' is (WL,GF)
00990    iret = firstDifferences(which);
00991    if(iret) return iret;
00992 
00993    if(cfg(Debug) >= 5) dumpSegments(string("D")+which,2,true); // DWL DGF
00994 
00995    // scan the first differences, eliminate outliers and
00996    // break into segments where there are WL slips.
00997    limit = (which == string("WL") ? WLobviousNwlLimit : GFobviousNwlLimit);
00998    it = SegList.begin();
00999    nok = 0;
01000    igood = -1;
01001    outlier = false;
01002    for(i=0; i<size(); i++) {
01003       if(i < it->nbeg) {
01004          outlier = false;
01005          continue;
01006       }
01007       if(i > it->nend) {                  // change segments
01008          if(outlier) {
01009             if(spdvector[ibad].flag & OK) nok--;
01010             spdvector[ibad].flag = BAD;
01011             learn[string("points deleted: ") + which + string(" slip outlier")]++;
01012             outlier = false;
01013          }
01014          it->npts = nok;
01015          // update nbeg and nend
01016          while(it->nbeg < it->nend
01017             && it->nbeg < size()
01018             && !(spdvector[it->nbeg].flag & OK) ) it->nbeg++;
01019          while(it->nend > it->nbeg
01020             && it->nend > 0
01021             && !(spdvector[it->nend].flag & OK) ) it->nend--;
01022          it++;
01023          if(it == SegList.end())
01024             return ReturnOK;
01025          nok = 0;
01026       }
01027 
01028       if(!(spdvector[i].flag & OK))
01029          continue;
01030       nok++;                                   // nok = # good points in segment
01031 
01032       if(igood == -1) igood = i;               // igood is index of last good point
01033 
01034       if(fabs(spdvector[i].data[A1]) > limit) {// found an outlier (1st diff, cycles)
01035          outlier = true;
01036          ibad = i;                             // ibad is index of last bad point
01037       }
01038       else if(outlier) {                       // this point good, but not past one(s)
01039          for(j=igood+1; j<ibad; j++) {
01040             if(spdvector[j].flag & OK)
01041                nok--;
01042             if(spdvector[j].flag & DETECT)
01043                log << "Warning - found an obvious slip, "
01044                   << "but marking BAD a point already marked with slip "
01045                   << GDCUnique << " " << sat
01046                   << " " << time(j).printf(outFormat) << " " << j << endl;
01047             spdvector[j].flag = BAD;             // mark all points between as bad
01048             learn[string("points deleted: ") + which + string(" slip outlier")]++;
01049          }
01050 
01051             // create a new segment, starting at the last outlier
01052          it->npts = nok-2;
01053          // WL slip gross  OR  GF slip gross
01054          it = createSegment(it,ibad,which+string(" slip gross"));
01055 
01056             // mark it
01057          spdvector[ibad].flag |= (which == string("WL") ? WLDETECT : GFDETECT);
01058 
01059             // change the bias in the new segment
01060          if(which == "WL") {
01061             wlbias = spdvector[ibad].data[P1];
01062             it->bias1 = long(wlbias+(wlbias > 0 ? 0.5 : -0.5));   // WL bias (NWL)
01063          }
01064          if(which == "GF")
01065             it->bias2 = spdvector[ibad].data[L2];                 // GFP bias
01066 
01067             // prep for next point
01068          nok = 2;
01069          outlier = false;
01070          igood = ibad;
01071       }
01072       else
01073          igood = i;
01074    }
01075    it->npts = nok;
01076 
01077    return ReturnOK;
01078 }
01079 catch(Exception& e) { GPSTK_RETHROW(e); }
01080 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01081 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01082 }
01083 
01084 //------------------------------------------------------------------------------------
01085 // compute first differences of data array(s) for WL and GF gross slip detection.
01086 // for WL difference the WLbias (P1); for GF, the GFP (L2) and the GFR (P2)
01087 // Store results in A1, and for GF put the range difference in A2
01088 int GDCPass::firstDifferences(string which) throw(Exception)
01089 {
01090 try {
01091    //if(A1.size() != size()) return FatalProblem;
01092 
01093    int i,iprev=-1;
01094 
01095    for(i=0; i<size(); i++) {
01096       // ignore bad data
01097       if(!(spdvector[i].flag & OK)) {
01098          spdvector[i].data[A1] = spdvector[i].data[A2] = 0.0;
01099          continue;
01100       }
01101 
01102       // compute first differences - 'change the arrays' A1 and A2
01103       if(which == string("WL")) {
01104          if(iprev == -1)
01105             spdvector[i].data[A1] = 0.0;
01106          else
01107             spdvector[i].data[A1] =
01108                (spdvector[i].data[P1] - spdvector[iprev].data[P1]) /
01109                   (spdvector[i].ndt-spdvector[iprev].ndt);
01110       }
01111       else if(which == string("GF")) {
01112          if(iprev == -1)            // first difference not defined at first point
01113             spdvector[i].data[A1] = spdvector[i].data[A2] = 0.0;
01114          else {
01115             // compute first difference of L1 = raw residual GFP-GFR
01116             spdvector[i].data[A1] =
01117                (spdvector[i].data[L1] - spdvector[iprev].data[L1]);
01118                   // 040809 should this be divided by delta N?
01119                   // / (spdvector[i].ndt-spdvector[iprev].ndt);
01120             // compute first difference of L2 = GFP
01121             spdvector[i].data[A2] =
01122                (spdvector[i].data[L2] - spdvector[iprev].data[L2]);
01123                   // 040809 should this be divided by delta N?
01124                   // / (spdvector[i].ndt-spdvector[iprev].ndt);
01125          }
01126       }
01127 
01128       // go to next point
01129       iprev = i;
01130    }
01131 
01132    return ReturnOK;
01133 }
01134 catch(Exception& e) { GPSTK_RETHROW(e); }
01135 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01136 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01137 }
01138 
01139 //------------------------------------------------------------------------------------
01140 // for one segment, compute conventional statistics on the
01141 // WL bias and count the number of good points
01142 void GDCPass::WLcomputeStats(list<Segment>::iterator& it) throw(Exception)
01143 {
01144 try {
01145    // compute WLStats
01146    it->WLStats.Reset();
01147    it->npts = 0;
01148 
01149    // loop over data, adding to Stats, and counting good points
01150    for(int i=it->nbeg; i<=it->nend; i++) {
01151       if(!(spdvector[i].flag & OK)) continue;
01152       it->WLStats.Add(spdvector[i].data[P1] - it->bias1);
01153       it->npts++;
01154    }
01155 
01156    // eliminate segments with too few points
01157    if(it->npts < int(cfg(MinPts)))
01158       deleteSegment(it,"insufficient data in segment");
01159 }
01160 catch(Exception& e) { GPSTK_RETHROW(e); }
01161 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01162 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01163 }
01164 
01165 //------------------------------------------------------------------------------------
01166 // for one segment, compute conventional statistics on the
01167 // WL bias, remove small segments, and mark bad points that lie outside N*sigma
01168 void GDCPass::WLsigmaStrip(list<Segment>::iterator& it) throw(Exception)
01169 {
01170 try {
01171    bool outlier,haveslip;
01172    unsigned short slip;
01173    int i,slipindex;
01174    double wlbias,nsigma,ave;
01175 
01176    // use robust stats on small segments, for big ones stick with conventional
01177 
01178    if(it->npts < cfg(WLNptsOutlierStats)) {   // robust
01179       // 'change the arrays' A1 and A2; they will be overwritten by WLstatSweep
01180       // but then dumped as DSCWLF...
01181       // use temp vectors so they can be passed to Robust::MAD and Robust::MEstimate
01182       int j,k;
01183       double median,mad;
01184       vector<double> vecA1,vecA2;      // use these temp, for Robust:: calls
01185 
01186       // put wlbias in vecA1, but without gaps: let j index good points only from nbeg
01187       for(j=i=it->nbeg; i<=it->nend; i++) {
01188          if(!(spdvector[i].flag & OK)) continue;
01189          wlbias = spdvector[i].data[P1] - it->bias1;
01190          vecA1.push_back(wlbias);
01191          vecA2.push_back(0.0);
01192          j++;
01193       }
01194 
01195       mad = Robust::MAD(&(vecA1[0]),j-it->nbeg,median,true);
01196       nsigma = cfg(WLNSigmaStrip) * mad;
01197       ave = Robust::MEstimate(&(vecA1[0]),j-it->nbeg,median,mad,&(vecA2[0]));
01198 
01199       // change the array : A1 is wlbias, A2 (output) will contain the weights
01200       // copy temps out into A1 and A2
01201       for(k=0,i=it->nbeg; i<j; k++,i++) {
01202          spdvector[i].data[A1] = vecA1[k];
01203          spdvector[i].data[A2] = vecA2[k];
01204       }
01205 
01206       haveslip = false;
01207       for(j=i=it->nbeg; i<=it->nend; i++) {
01208          if(!(spdvector[i].flag & OK)) continue;
01209 
01210          wlbias = spdvector[i].data[P1] - it->bias1;
01211 
01212          // TD ? use weights at all? they remove a lot of points
01213          // TD add absolute limit?
01214          if(fabs(wlbias-ave) > nsigma ||
01215                spdvector[j].data[A2] < cfg(WLRobustWeightLimit))
01216             outlier = true;
01217          else
01218             outlier = false;
01219 
01220          // remove points by sigma stripping
01221          if(outlier) {
01222             if(spdvector[i].flag & DETECT || i == it->nbeg) {
01223                haveslip = true;
01224                slipindex = i;        // mark
01225                slip = spdvector[i].flag; // save to put on first good point
01226                //log << "Warning - marking a slip point BAD in WL sigma strip "
01227                //   << GDCUnique << " " << sat
01228                //   << " " << time(i).printf(outFormat) << " " << i << endl;
01229             }
01230             spdvector[i].flag = BAD;
01231             learn["points deleted: WL sigma stripping"]++;
01232             it->npts--;
01233             it->WLStats.Subtract(wlbias);
01234          }
01235          else if(haveslip) {
01236             spdvector[i].flag = slip;
01237             haveslip = false;
01238          }
01239 
01240          if(cfg(Debug) >= 6) {
01241             log << "DSCWLR " << GDCUnique << " " << sat
01242             << " " << it->nseg
01243             << " " << time(i).printf(outFormat)
01244             << fixed << setprecision(3)
01245             << " " << setw(3) << spdvector[i].flag
01246             << " " << setw(13) << spdvector[j].data[A1] // wlbias
01247             << " " << setw(13) << fabs(wlbias-ave)
01248             << " " << setw(5) << spdvector[j].data[A2]  // 0 <= weight <= 1
01249             << " " << setw(3) << i
01250             << (outlier ? " outlier" : "");
01251             if(i == it->nbeg) log
01252                << " " << setw(13) << it->bias1
01253                << " " << setw(13) << it->bias2;
01254             log << endl;
01255          }
01256 
01257          j++;
01258       }
01259 
01260    }
01261    else {                                             // conventional
01262 
01263       //nsigma = WLsigmaStrippingNsigmaLimit * it->WLStats.StdDev();
01264       nsigma = cfg(WLNSigmaStrip) * it->WLStats.StdDev();
01265 
01266       haveslip = false;
01267       ave = it->WLStats.Average();
01268       for(i=it->nbeg; i<=it->nend; i++) {
01269          if(!(spdvector[i].flag & OK)) continue;
01270 
01271          wlbias = spdvector[i].data[P1] - it->bias1;
01272 
01273          // remove points by sigma stripping
01274          if(fabs(wlbias-ave) > nsigma) { // TD add absolute limit?
01275             if(spdvector[i].flag & DETECT) {
01276                haveslip = true;
01277                slipindex = i;        // mark
01278                slip = spdvector[i].flag; // save to put on first good point
01279                //log << "Warning - marking a slip point BAD in WL sigma strip "
01280                //   << GDCUnique << " " << sat
01281                //   << " " << time(i).printf(outFormat) << " " << i << endl;
01282             }
01283             spdvector[i].flag = BAD;
01284             learn["points deleted: WL sigma stripping"]++;
01285             it->npts--;
01286             it->WLStats.Subtract(wlbias);
01287          }
01288          else if(haveslip) {
01289             spdvector[i].flag = slip;
01290             haveslip = false;
01291          }
01292 
01293       }  // loop over points in segment
01294    }
01295 
01296    // change nbeg, but don't change the bias
01297    if(haveslip) {
01298       it->nbeg = slipindex;
01299       //wlbias = spdvector[slipindex].data[P1];
01300       //it->bias1 = long(wlbias+(wlbias > 0 ? 0.5 : -0.5));
01301    }
01302 
01303    // again
01304    if(it->npts < int(cfg(MinPts)))
01305       deleteSegment(it,"WL sigma stripping");
01306    else {
01307       // update nbeg and nend // TD add limit 0 size()
01308       while(it->nbeg < it->nend && !(spdvector[it->nbeg].flag & OK)) it->nbeg++;
01309       while(it->nend > it->nbeg && !(spdvector[it->nend].flag & OK)) it->nend--;
01310    }
01311 
01312 }
01313 catch(Exception& e) { GPSTK_RETHROW(e); }
01314 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01315 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01316 }
01317 
01318 //------------------------------------------------------------------------------------
01319 // in the given segment, compute statistics on the WL bias using a
01320 // 'two-paned sliding window', each pane of width 'width' good points.
01321 // store the results in the parallel (to SatPass::data) arrays A1,A2.
01322 int GDCPass::WLstatSweep(list<Segment>::iterator& it, int width) throw(Exception)
01323 {
01324 try {
01325    int iminus,i,iplus;
01326    double wlbias,test,limit;
01327    Stats<double> pastStats, futureStats;
01328 
01329    // ignore empty segments
01330    if(it->npts == 0) return ReturnOK;
01331    it->WLsweep = true;
01332 
01333       // Cartoon of the 'two-pane moving window'
01334       // windows:  'past window'      'future window'
01335       // stats  :  --- pastStats----  ----futureStats--
01336       // data   : (x x x x x x x x x)(x x x x x x x x x) x ...
01337       //           |               |  |                  |
01338       // indexes: iminus          i-1 i                 iplus
01339 
01340    // start with the window 'squashed' to one point - the first one
01341    iminus = i = iplus = it->nbeg;
01342 
01343    // fill up the future window to size 'width', but don't go beyond the segment
01344    while(futureStats.N() < width && iplus <= it->nend) {
01345       if(spdvector[iplus].flag & OK) {                // add only good data
01346          futureStats.Add(spdvector[iplus].data[P1] - it->bias1);
01347       }
01348       iplus++;
01349    }
01350 
01351    // now loop over all points in the segment
01352    for(i=it->nbeg; i<= it->nend; i++) {
01353       if(!(spdvector[i].flag & OK))                      // add only good data
01354          continue;
01355 
01356       // compute test and limit
01357       test = 0;
01358       if(pastStats.N() > 0 && futureStats.N() > 0)
01359          test = fabs(futureStats.Average()-pastStats.Average());
01360       limit = ::sqrt(futureStats.Variance() + pastStats.Variance());
01361       // 'change the arrays' A1 and A2
01362       spdvector[i].data[A1] = test;
01363       spdvector[i].data[A2] = limit;
01364 
01365       wlbias = spdvector[i].data[P1] - it->bias1;        // debiased WLbias
01366 
01367       // dump the stats
01368       if(cfg(Debug) >= 6) log << "WLS " << GDCUnique
01369          << " " << sat << " " << it->nseg
01370          << " " << time(i).printf(outFormat)
01371          << fixed << setprecision(3)
01372          << " " << setw(3) << pastStats.N()
01373          << " " << setw(7) << pastStats.Average()
01374          << " " << setw(7) << pastStats.StdDev()
01375          << " " << setw(3) << futureStats.N()
01376          << " " << setw(7) << futureStats.Average()
01377          << " " << setw(7) << futureStats.StdDev()
01378          << " " << setw(9) << spdvector[i].data[A1]
01379          << " " << setw(9) << spdvector[i].data[A2]
01380          << " " << setw(9) << wlbias
01381          << " " << setw(3) << i
01382          << endl;
01383 
01384       // update stats :
01385       // move point i from future to past, ...
01386       futureStats.Subtract(wlbias);
01387       pastStats.Add(wlbias);
01388       // ... and move iplus up by one (good) point, ...
01389       while(futureStats.N() < width && iplus <= it->nend) {
01390          if(spdvector[iplus].flag & OK) {
01391             futureStats.Add(spdvector[iplus].data[P1] - it->bias1);
01392          }
01393          iplus++;
01394       }
01395       // ... and move iminus up by one good point
01396       while(pastStats.N() > width && iminus <= it->nend) {// <= nend not really nec.
01397          if(spdvector[iminus].flag & OK) {
01398             pastStats.Subtract(spdvector[iminus].data[P1] - it->bias1);
01399          }
01400          iminus++;
01401       }
01402 
01403    }  // end loop over i=all points in segment
01404 
01405    return  ReturnOK;
01406 }
01407 catch(Exception& e) { GPSTK_RETHROW(e); }
01408 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01409 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01410 }
01411 
01412 //------------------------------------------------------------------------------------
01413 // look for slips in the WL using the results of WLstatSweep
01414 // if slip is close to either end (< window width), just chop off the small segment
01415 // recompute WLstats; when a slip is found, create a new segment
01416 int GDCPass::detectWLsmallSlips(void) throw(Exception)
01417 {
01418 try {
01419    int i,k,nok;
01420    double wlbias;
01421    list<Segment>::iterator it;
01422 
01423    // find first segment for which WLstatSweep was called
01424    it = SegList.begin();
01425    while(!it->WLsweep) {
01426       it++;
01427       if(it == SegList.end()) return ReturnOK;
01428    }
01429    it->WLStats.Reset();
01430 
01431    // loop over the data arrays - all segments
01432    i = it->nbeg;
01433    nok = 0;
01434    int halfwidth = int(cfg(WLSlipEdge));
01435    while(i < size())
01436    {
01437       // must skip segments for which WLstatSweep was not called
01438       while(i > it->nend || !it->WLsweep) {
01439          if(i > it->nend) {
01440             it->npts = nok;
01441             nok = 0;
01442          }
01443          it++;
01444          if(it == SegList.end()) return ReturnOK;
01445          i = it->nbeg;
01446          if(it->WLsweep) {
01447             it->WLStats.Reset();
01448          }
01449       }
01450 
01451       if(spdvector[i].flag & OK) {
01452          nok++;                                 // nok = # good points in segment
01453 
01454          if(nok == 1) {                         // change the bias, as WLStats reset
01455             wlbias = spdvector[i].data[P1];
01456             it->bias1 = long(wlbias+(wlbias > 0 ? 0.5 : -0.5));
01457          }
01458 
01459          //  condition 3 - near ends of segment?
01460          if(nok < halfwidth || (it->npts - nok) < halfwidth ) {
01461             // failed test 3 - near ends of segment
01462             // consider chopping off this end of segment - large limit?
01463             // TD must do something here ...
01464             if(cfg(Debug) >= 6) log << "too near end " << GDCUnique
01465                << " " << i << " " << nok << " " << it->npts-nok
01466                << " " << time(i).printf(outFormat)
01467                << " " << spdvector[i].data[A1] << " " << spdvector[i].data[A2]
01468                << endl;
01469          }
01470          else if(foundWLsmallSlip(it,i)) { // met condition 3
01471             // create new segment
01472             // TD what if nok < MinPts? -- cf detectGFsmallSlips
01473             k = it->npts;
01474             it->npts = nok;
01475             //log << "create new segment at i = " << i << " " << nok << "pts" << endl;
01476             it = createSegment(it,i,"WL slip small");
01477 
01478             // mark it
01479             spdvector[i].flag |= WLDETECT;
01480 
01481             // prep for next segment
01482             // biases remain the same in the new segment
01483             it->npts = k - nok;
01484             nok = 0;
01485             it->WLStats.Reset();
01486             wlbias = spdvector[i].data[P1]; // change the bias, as WLStats reset
01487             it->bias1 = long(wlbias+(wlbias > 0 ? 0.5 : -0.5));
01488          }
01489 
01490          it->WLStats.Add(spdvector[i].data[P1] - it->bias1);
01491 
01492       } // end if good data
01493 
01494       i++;
01495 
01496    }  // end loop over points in the pass
01497    it->npts = nok;
01498 
01499    return ReturnOK;
01500 }
01501 catch(Exception& e) { GPSTK_RETHROW(e); }
01502 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01503 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01504 }
01505 
01506 //------------------------------------------------------------------------------------
01507 // determine if a slip has been found at index i, in segment it
01508 // A1 = test = fabs(futureStats.Average()-pastStats.Average()) ~ step in ave WL
01509 // A2 = limit = sqrt(futureStats.Variance()+pastStats.Variance()) ~ noise in WL
01510 // ALL CONDITIONs needed for a slip to be detected:
01511 // 1. test must be > ~0.67(WLSlipSize) * WLwl
01512 // 2. test-limit must be > (WLSlipExcess) ie test
01513 // 3. slip must be far (>1/2 window) from either end
01514 // 4. test must be at a local maximum within ~ window width
01515 // 5. limit must be at a local minimum (")
01516 // 6. (test-limit)/limit > 1.2(WLSlipSeparation)
01517 // large limit (esp near end of a pass) means too much noise
01518 bool GDCPass::foundWLsmallSlip(list<Segment>::iterator& it, int i)
01519    throw(Exception)
01520 {
01521 try {
01522    const int minMaxWidth=int(cfg(WLSlipEdge)); // test 4,5, = ~~1/2 WLWindowWidth
01523    int j,jp,jm,pass4,pass5,Pass;
01524 
01525    // A1 = test = fabs(futureStats.Average() - pastStats.Average());
01526    // A2 = limit = ::sqrt(futureStats.Variance() + pastStats.Variance());
01527    // all units WL cycles
01528    double test = spdvector[i].data[A1];
01529    double lim = spdvector[i].data[A2];
01530 
01531    // 050109 if Debug=6, print only possible slips, if 7 print all
01532    bool isSlip=false;
01533    ostringstream oss;
01534    for(;;) {
01535 
01536       // Debug print - NB '>' is always pass, '<=' is fail....
01537       if(cfg(Debug) >= 6) oss << "WLslip " << GDCUnique
01538          << " " << sat << " " << setw(2) << it->nseg
01539          << " " << setw(3) << i
01540          << " " << time(i).printf(outFormat)
01541          //<< " " << it->npts << "pt"
01542          << fixed << setprecision(2)
01543          << " test=" << test << " lim=" << lim
01544          << " (1)" << spdvector[i].data[A1]
01545          << (spdvector[i].data[A1] > cfg(WLSlipSize) ? ">" : "<=")
01546          << cfg(WLSlipSize)
01547          << " (2)" << spdvector[i].data[A1]-spdvector[i].data[A2]
01548          << (spdvector[i].data[A1]-spdvector[i].data[A2]>cfg(WLSlipExcess)?">":"<=")
01549          << cfg(WLSlipExcess); // no endl
01550 
01551       // CONDITION 1  ||  CONDITION 2
01552       if(test <= cfg(WLSlipSize) || test-lim <= cfg(WLSlipExcess)) break;
01553 
01554       // CONDITIONs 4 and 5
01555       //         x
01556       //        x x
01557       //       x   x
01558       //      x     x
01559       //     x       x
01560       //    x         x
01561       // ------------------------
01562       //      jp=012345
01563       // do for WLSlipEdge pts on each side of point; best score is pass=2*WLSlipEdge
01564       double slope=(test-lim)/(8.0*minMaxWidth);
01565       j = pass4 = pass5 = Pass = 0;
01566       jp = jm = i;
01567       do {
01568          // find next good point in future
01569          do { jp++; } while(jp < it->nend && !(spdvector[jp].flag & OK));
01570          if(jp >= it->nend) break;
01571             // CONDITION 4: test(A1) is a local maximum
01572          if(spdvector[i].data[A1]-spdvector[jp].data[A1] > j*slope) pass4++;
01573             // CONDITION 5: limit(A2) is a local minimum
01574          if(spdvector[i].data[A2]-spdvector[jp].data[A2] < -j*slope) pass5++;
01575 
01576          // find next good point in past
01577          do { jm--; } while(jm > it->nbeg && !(spdvector[jm].flag & OK));
01578          if(jm <= it->nbeg) break;
01579             // CONDITION 4: test(A1) is a local maximum
01580          if(spdvector[i].data[A1]-spdvector[jm].data[A1] > j*slope) pass4++;
01581             // CONDITION 5: limit(A2) is a local minimum
01582          if(spdvector[i].data[A2]-spdvector[jm].data[A2] < -j*slope) pass5++;
01583 
01584       } while(++j < minMaxWidth);
01585 
01586       // perfect = 2*minMaxWidth; allow 1 miss..
01587       if(pass4 >= 2*minMaxWidth-1) Pass++;
01588       if(cfg(Debug) >= 6) oss << " (4)" << pass4
01589          << (pass4 >= 2*minMaxWidth-1 ? ">" : "<=") << 2*minMaxWidth-2;
01590       if(pass5 >= 2*minMaxWidth-1) Pass++;
01591       if(cfg(Debug) >= 6) oss << " (5)" << pass5
01592          << (pass5 >= 2*minMaxWidth-1 ? ">" : "<=") << 2*minMaxWidth-2;
01593 
01594       // CONDITION 7
01595       double ratio=(test-lim)/lim;
01596       if(cfg(Debug) >= 6) oss << " (6)" << ratio
01597          << (ratio > cfg(WLSlipSeparation) ? ">" : "<=") << cfg(WLSlipSeparation);
01598       if(ratio > cfg(WLSlipSeparation) ) Pass++;
01599 
01600       if(Pass == 3) {
01601          if(cfg(Debug) >= 6) oss << " possible WL slip";
01602          isSlip = true;
01603       }
01604 
01605       break;
01606    }  // end for(;;)
01607 
01608    if(cfg(Debug) >= 6) log << oss.str() << endl;
01609    return isSlip;
01610 }
01611 catch(Exception& e) { GPSTK_RETHROW(e); }
01612 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01613 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01614 }
01615 
01616 //------------------------------------------------------------------------------------
01617 //------------------------------------------------------------------------------------
01618 // estimate slips and adjust biases appropriately - ie fix slips - for both WL and GF
01619 // merge all data into one segment
01620 int GDCPass::fixAllSlips(string which) throw(Exception)
01621 {
01622 try {
01623    // find the largest segment and start there, always combine the largest with its
01624    // largest neighbor
01625    int i,nmax,ifirst;
01626    list<Segment>::iterator it, kt;
01627 
01628    // loop over all segments, erasing empty ones
01629    it = SegList.begin();
01630    while(it != SegList.end()) {
01631       if(it->npts == 0)
01632          it = SegList.erase(it);
01633       else
01634          it++;
01635    }
01636 
01637    if(SegList.empty())
01638       return NoData;
01639 
01640    // find the largest segment
01641    for(kt=SegList.end(),nmax=0,it=SegList.begin(); it != SegList.end(); it++) {
01642       if(it->npts > nmax) {
01643          nmax = it->npts;
01644          kt = it;
01645       }
01646    }
01647 
01648    // fix all the slips, starting with the largest segment
01649    // this will merge all segments into one
01650    GDCUniqueFix = 0;
01651    while(kt != SegList.end()) {
01652       fixOneSlip(kt,which);
01653    }
01654 
01655    // TD here to return should be a separate call...
01656 
01657    // now compute stats for the WL for the (single segment) whole pass
01658    kt = SegList.begin();
01659    if(which == string("WL")) {                                    // WL
01660       WLPassStats.Reset();
01661       for(i=kt->nbeg; i <= kt->nend; i++) {
01662          if(!(spdvector[i].flag & OK)) continue;
01663          WLPassStats.Add(spdvector[i].data[P1] - kt->bias1);
01664       }
01665       // NB Now you have a measure of range noise for the whole pass :
01666       // sigma(WLbias) ~ sigma(WLrange) = 0.71*sigma(range), so
01667       // range noise = WLPassStats.StdDev() * wlwl / 0.71;  // meters
01668       // 0.71 / wlwl = 0.83
01669 
01670       // TD mark the first slip 'fixed' - unmark it - or something
01671    }
01672    // change the biases - reset the GFP bias so that it matches the GFR
01673    else {                                                         // GF
01674       //dumpSegments("GFFbefRebias",2,true); //temp
01675       for(ifirst=-1,i=kt->nbeg; i <= kt->nend; i++) {
01676          if(!(spdvector[i].flag & OK)) continue;
01677          if(ifirst == -1) {
01678             ifirst = i;
01679             kt->bias2 = spdvector[ifirst].data[L2] + spdvector[ifirst].data[P2];
01680             kt->bias1 = spdvector[ifirst].data[P1];
01681          }
01682          // change the data - recompute GFR-GFP so it has one consistent bias
01683          spdvector[i].data[L1] = spdvector[i].data[L2] + spdvector[i].data[P2];
01684       }
01685    }
01686 
01687    if(cfg(Debug) >= 3) dumpSegments(which + string("F"),2,true);   // DSCWLF DSCGFF
01688 
01689    return ReturnOK;
01690 }
01691 catch(Exception& e) { GPSTK_RETHROW(e); }
01692 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01693 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01694 }
01695 
01696 //------------------------------------------------------------------------------------
01697 // called by fixAllSlips
01698 // assume there are no empty segments in the list
01699 void GDCPass::fixOneSlip(list<Segment>::iterator& kt, string which) throw(Exception)
01700 {
01701 try {
01702    if(kt->npts == 0) { kt++; return; }
01703 
01704    list<Segment>::iterator left,right,it;
01705 
01706    // kt points to the biggest segment
01707    // define left and right to be the two segments on each side of the slip to be
01708    // fixed; assume there are no empty segments in the list
01709    right = left = kt;
01710 
01711    // choose the next segment on the right of kt
01712    right++;
01713 
01714    // choose the next segment on the left of kt
01715    if(kt != SegList.begin())
01716       left--;
01717    else
01718       left = SegList.end();            // nothing on the left
01719 
01720    // no segment left of kt and no segment right of kt - nothing to do
01721    if(left == SegList.end() && right == SegList.end()) {
01722       kt++;
01723       return;
01724    }
01725 
01726    // Always define kt to == left, as it will be returned and right will be erased.
01727    if(left == SegList.end()) {         // no segment on left
01728       left = kt;
01729    }
01730    else if(right == SegList.end()      // no segment on right
01731       || left->npts >= right->npts) {  // or left is the bigger segment
01732       right = kt;
01733       kt = left;                       // fix between left and kt
01734    }
01735    else {                              // left and right exist, and right is bigger
01736       left = kt;                       // fix between kt and right
01737    }
01738 
01739    // fix the slip between left and right, making data in 'right' part of 'left'
01740    if(which == string("WL"))
01741       WLslipFix(left,right);
01742    else
01743       GFslipFix(left,right);
01744 
01745    left->npts += right->npts;
01746    left->nend = right->nend;
01747 
01748    // always delete right, otherwise on return kt(==left) will be invalid
01749    // (ignore return value = iterator to first element after the one erased)
01750    SegList.erase(right);
01751 
01752    return;
01753 }
01754 catch(Exception& e) { GPSTK_RETHROW(e); }
01755 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01756 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01757 }
01758 
01759 //------------------------------------------------------------------------------------
01760 // called by fixOneSlip
01761 void GDCPass::WLslipFix(list<Segment>::iterator& left,
01762                         list<Segment>::iterator& right)
01763 throw(Exception)
01764 {
01765 try {
01766    int i;
01767 
01768    GDCUniqueFix++;
01769 
01770    // full slip
01771    double dwl = right->bias1 + right->WLStats.Average()
01772          - (left->bias1 + left->WLStats.Average());
01773    long nwl = long(dwl + (dwl > 0 ? 0.5 : -0.5));
01774 
01775       // TD ? test gap size?
01776    //if(cfg(DT)*(right->nbeg - left->nend) > cfg(MaxGap)) break;
01777 
01778       // test that total variance is small
01779    //if(::sqrt(left->WLStats.Variance() + right->WLStats.Variance())
01780    //   / (left->WLStats.N() + right->WLStats.N()) < cfg(WLFixSigma)) {
01781    //   log << "Cannot fix WL slip (noisy) at " << right->nbeg
01782    //      << " " << time(right->nbeg).printf(outFormat)
01783    //      << endl;
01784    //   break;
01785    //}
01786 
01787       // TD ? test fractional part of offset fabs
01788    //if(fabs(dwl-nwl) > cfg(WLFixFrac)) break;
01789 
01790    if(cfg(Debug) >= 6) log << "Fix " << GDCUnique << " " << sat << " " << GDCUniqueFix
01791       << " WL " << time(right->nbeg).printf(outFormat)
01792       << " " << nwl           // put integer fix after time, all 'Fix'
01793       << " " << left->nseg << "-" << right->nseg
01794       << fixed << setprecision(2)
01795       << " right: " << right->bias1 << " + " << right->WLStats.Average()
01796       << " - left: " << left->bias1 << " + " << left->WLStats.Average()
01797       << " = " << dwl << " " << nwl
01798       << endl;
01799 
01800    // now do the fixing - change the data in the right segment to match left's
01801    for(i=right->nbeg; i<=right->nend; i++) {
01802       //if(!(spdvector[i].flag & OK)) continue;
01803       spdvector[i].data[P1] -= nwl;                                 // WLbias
01804       spdvector[i].data[L2] -= nwl * wl2;                           // GFP
01805       // add to WLStats
01806       //if(!(spdvector[i].flag & OK)) continue;
01807       //left->WLStats.Add(spdvector[i].data[P1] - left->bias1);
01808    }
01809 
01810    // fix the slips beyond the 'right' segment.
01811    // change the data in the GFP, and change both the data and the bias in the WL.
01812    // this way, WLStats is still valid, but if we change the GF bias, we will lose
01813    // that information before the GF slips get fixed.
01814    list<Segment>::iterator it = right;
01815    for(it++; it != SegList.end(); it++) {
01816       // Use real, not int, nwl b/c rounding error in a pass with many slips
01817       // can build up and produce errors.
01818       it->bias1 -= dwl;
01819       for(i=it->nbeg; i<=it->nend; i++) {
01820          //if(!(spdvector[i].flag & OK)) continue;                 // TD don't?
01821          spdvector[i].data[P1] -= nwl;                                 // WLbias
01822          spdvector[i].data[L2] -= nwl * wl2;                           // GFP
01823       }
01824    }
01825 
01826    // Add to slip list
01827    Slip newSlip(right->nbeg);
01828    newSlip.NWL = nwl;
01829    newSlip.msg = "WL";
01830    SlipList.push_back(newSlip);
01831 
01832    // mark it
01833    spdvector[right->nbeg].flag |= WLFIX;
01834 
01835    return;
01836 }
01837 catch(Exception& e) { GPSTK_RETHROW(e); }
01838 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01839 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01840 }
01841 
01842 //------------------------------------------------------------------------------------
01843 // fix one slip in the geometry-free phase
01844 // called by fixOneSlip
01845 void GDCPass::GFslipFix(list<Segment>::iterator& left,
01846                         list<Segment>::iterator& right) throw(Exception)
01847 {
01848 try {
01849       // use this number of data points on each side of slip
01850    const int Npts=int(cfg(GFFixNpts));
01851    int i,nb,ne,nl,nr,ilast;
01852    long n1,nadj;              // slip magnitude (cycles)
01853    double dn1,dnGFR;
01854    Stats<double> Lstats,Rstats;
01855 
01856    GDCUniqueFix++;
01857 
01858 //temp
01859 //log << "GFslipFix " << GDCUniqueFix << " biases L: " << left->bias2
01860 //    << " R: " << right->bias2 << endl;
01861    // find Npts points on each side of slip
01862    nb = left->nend;
01863    i = 1;
01864    nl = 0;
01865    ilast = -1;                               // ilast is last good point before slip
01866    while(nb > left->nbeg && i < Npts) {
01867       if(spdvector[nb].flag & OK) {
01868          if(ilast == -1) ilast = nb;
01869          i++; nl++;
01870          Lstats.Add(spdvector[nb].data[L1] - left->bias2);
01871 //log << "LDATA " << nb << " " << spdvector[nb].data[L1]-left->bias2 << endl;
01872       }
01873       nb--;
01874    }
01875    ne = right->nbeg;
01876    i = 1;
01877    nr = 0;
01878    while(ne < right->nend && i < Npts) {
01879       if(spdvector[ne].flag & OK) {
01880          i++; nr++;
01881          Rstats.Add(spdvector[ne].data[L1] - right->bias2);
01882 //log << "RDATA " << ne << " " << spdvector[ne].data[L1]-right->bias2 << endl;
01883       }
01884       ne++;
01885    }
01886 
01887    // first estimate of n1, without biases
01888    // need to use the GFR-GFP estimate here, and limit |nadj| to be well within
01889    // sigmas on the stats, b/c when ionosphere is very active, GFP and GFR will both
01890    // vary sharply, and fitting a polynomial to GFP is a BAD thing to do....
01891    // ultimately, GFR-GFP is accurate but noisy.
01892    // rms rof should tell you how much weight to put on rof
01893    // larger rof -> smaller npts and larger degree
01894    dn1 = spdvector[right->nbeg].data[L2] - right->bias2
01895          - (spdvector[ilast].data[L2] - left->bias2);
01896    // this screws up most fixes
01897    //dn1 = Rstats.Average() - right->bias2 - (Lstats.Average() - left->bias2);
01898    n1 = long(dn1 + (dn1 > 0 ? 0.5 : -0.5));
01899 
01900    // TD worry about too small pieces - nr or nl too small
01901 
01902    // estimate the slip using polynomial fits - this prints GFE data
01903    nadj = EstimateGFslipFix(left,right,nb,ne,n1);
01904 
01905    // adjust the adjustment if it is not consistent with Lstats vs Rstats
01906    // dn1+nadj                       - a. current best estimate
01907    // Rstats.Averge()-Lstats.Average - b. estimate from stats on GFR-GFP across slip
01908    // difference should be consistent with R/Lstats.StdDev
01909    // if not, replace nadj with b. - dn1
01910    dnGFR = Rstats.Average() - Lstats.Average();
01911 // temp add factor 10.*
01912    if(fabs(n1+nadj-dnGFR) > 10.*(Rstats.StdDev()+Lstats.StdDev())) {
01913       if(cfg(Debug) >= 6)
01914          log << "GFRadjust " << GDCUnique << " " << sat << " " << GDCUniqueFix
01915          << " GF " << time(right->nbeg).printf(outFormat)
01916          << fixed << setprecision(2)
01917          << " dbias(GFR): " << dnGFR
01918          << " n1+nadj: " << n1+nadj;
01919 
01920       nadj = long(dnGFR+(dnGFR > 0 ? 0.5 : -0.5)) - n1;
01921 
01922       if(cfg(Debug) >= 6)
01923          log << " new n1+nadj: " << n1+nadj << endl;
01924    }
01925 
01926    // output result
01927    if(cfg(Debug) >= 6) {
01928       log << "Fix " << GDCUnique << " " << sat << " " << GDCUniqueFix
01929       << " GF " << time(right->nbeg).printf(outFormat)
01930       << " " << nadj           // put integer fix after time, all 'Fix'
01931       << fixed << setprecision(2)
01932       << " dbias: " << right->bias2 - left->bias2
01933       << ", dn1: " << dn1 << ", n1: " << n1 << ", adj: " << nadj
01934       << " indexes " << nb << " " << ne << " " << nl << " " << nr
01935       << " segs " << left->nseg << " " << right->nseg
01936       << " GFR-GFP:L: "
01937       << Lstats.N() << " " << Lstats.Average() << " " << Lstats.StdDev()
01938       << "    R: "
01939       << Rstats.N() << " " << Rstats.Average() << " " << Rstats.StdDev()
01940       << " tests " << n1+nadj-dnGFR << " " << Rstats.StdDev()+Lstats.StdDev()
01941       << endl;
01942    }
01943 
01944    // full slip, including biases
01945    dn1 += right->bias2 - left->bias2;
01946    n1 = long(dn1 + (dn1 > 0 ? 0.5 : -0.5));
01947    n1 += nadj;
01948 
01949    // now do the fixing : 'change the data' within right segment
01950    // and through the end of the pass, to fix the slip
01951    for(i=right->nbeg; i<size(); i++) {
01952       //if(!(spdvector[i].flag & OK)) continue;                 // TD? don't?
01953       //spdvector[i].data[P1] -= nwl;                           // no change to WLbias
01954       spdvector[i].data[L2] -= n1;                              // GFP
01955       spdvector[i].data[L1] -= n1;                              // GFR+GFP
01956    }
01957 
01958    // 'change the bias'  for all segments in the future (although right to be deleted)
01959    list<Segment>::iterator kt;
01960    for(kt=right; kt != SegList.end(); kt++)
01961       kt->bias2 -= n1;
01962 
01963    // Add to slip list, but if one exists with same time tag, use it instead
01964    list<Slip>::iterator jt;
01965    for(jt=SlipList.begin(); jt != SlipList.end(); jt++)
01966       if(jt->index == right->nbeg) break;
01967 
01968    if(jt == SlipList.end()) {
01969       Slip newSlip(right->nbeg);
01970       newSlip.N1 = -n1;
01971       newSlip.msg = "GF only";
01972       SlipList.push_back(newSlip);
01973    }
01974    else {
01975       jt->N1 = -n1;
01976       jt->msg += string(" GF");
01977    }
01978 
01979    // mark it
01980    spdvector[right->nbeg].flag |= GFFIX;
01981 
01982    return;
01983 }
01984 catch(Exception& e) { GPSTK_RETHROW(e); }
01985 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
01986 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
01987 }
01988 
01989 //------------------------------------------------------------------------------------
01990 // called by GFslipFix
01991 // estimate GF slip using polynomial fit to data around it
01992 long GDCPass::EstimateGFslipFix(list<Segment>::iterator& left,
01993                                list<Segment>::iterator& right,
01994                                int nb, int ne, long n1)
01995 throw(Exception)
01996 {
01997 try {
01998    bool quit;
01999    int i,k,in[3];
02000    double rof,rmsrof[3];
02001    PolyFit<double> PF[3];
02002 
02003    // start at zero and limit |nadj| to ...TD
02004    long nadj = 0;
02005 
02006    // use a little indirect indexing array to avoid having to copy PolyFit objects....
02007    for(k=0; k<3; k++) {
02008       in[k]=k;
02009       PF[in[k]].Reset(int(cfg(GFFixDegree)));
02010    }
02011 
02012    while(1) {
02013       // compute 3 polynomial fits to this data, with slips of
02014       // (nadj-1, nadj and nadj+1) wavelengths added to left segment
02015       for(k=0; k<3; k++) {
02016          if(PF[in[k]].N() > 0) continue;
02017 
02018          // add all the data
02019          for(i=nb; i<=ne; i++) {
02020             if(!(spdvector[i].flag & OK)) continue;
02021             PF[in[k]].Add(
02022                // data
02023                spdvector[i].data[L2]
02024                // - (either               left bias - poss. slip : right bias)
02025                   - (i < right->nbeg ? left->bias2-n1-(nadj+k-1) : right->bias2),
02026                //  use a debiased count
02027                spdvector[i].ndt - spdvector[nb].ndt
02028             );
02029          }
02030 
02031          // TD check that it not singular
02032 
02033          // compute RMS residual of fit
02034          rmsrof[in[k]] = 0.0;
02035          for(i=nb; i<=ne; i++) {
02036             if(!(spdvector[i].flag & OK)) continue;
02037             rof =    // data minus fit
02038                spdvector[i].data[L2]
02039                   - (i < right->nbeg ? left->bias2-n1-(nadj+k-1) : right->bias2)
02040                - PF[in[k]].Evaluate(spdvector[i].ndt - spdvector[nb].ndt);
02041             rmsrof[in[k]] += rof*rof;
02042          }
02043          rmsrof[in[k]] = ::sqrt(rmsrof[in[k]]);
02044 
02045       }  // end loop over fits
02046 
02047       // the value of this is questionable, b/c with active ionosphere the real
02048       // GFP is NOT smooth
02049       for(quit=false,k=0; k<3; k++) if(rmsrof[in[k]] > cfg(GFFixMaxRMS)) {
02050          log << "Warning - large RMS ROF in GF slip fix at in,k = "
02051              << in[k] << " " << k << " " << rmsrof[in[k]] << " abort." << endl;
02052          quit = true;
02053       }
02054       if(quit) break;
02055 
02056       //if(cfg(Debug) >= 6) {
02057       //   log << "Fix GF slip RMSROF : adj: " << nadj;
02058       //   for(i=0; i<3; i++) log << " " << rmsrof[in[i]];
02059       //   // below log << endl;
02060       //}
02061 
02062       // three cases: (TD - exceptions?) :
02063       // rmsrof: 0 > 1 < 2   good
02064       //         0 > 1 > 2   shift 0,1,2 to 1,2,3
02065       //         0 < 1 < 2   shift 0,1,2 to -1,0,1
02066       //         0 < 1 > 2   local max! - ??
02067       if(rmsrof[in[0]] > rmsrof[in[1]]) {
02068          if(rmsrof[in[1]] < rmsrof[in[2]]) { // local min - done
02069             //if(cfg(Debug) >= 6) log << " done." << endl;
02070             break;
02071          }
02072          else {                              // shift 0,1,2 to 1,2,3
02073             k = in[0];
02074             in[0] = in[1];
02075             in[1] = in[2];
02076             in[2] = k;
02077             PF[in[2]].Reset();
02078             nadj += 1;
02079             //if(cfg(Debug) >= 6) log << " shift left" << endl;
02080          }
02081       }
02082       else {
02083          if(rmsrof[in[1]] < rmsrof[in[2]]) { // shift 0,1,2 to -1,0,1
02084             k = in[2];
02085             in[2] = in[1];
02086             in[1] = in[0];
02087             in[0] = k;
02088             PF[in[0]].Reset();
02089             nadj -= 1;
02090             //if(cfg(Debug) >= 6) log << " shift right" << endl;
02091          }
02092          else {                              // local max
02093             log << "Warning - local maximum in RMS residuals in EstimateGFslipFix"
02094                << endl;
02095             // TD do something
02096             break;
02097          }
02098       }
02099 
02100    }  // end while loop
02101 
02102    // dump the raw data with all the fits
02103    if(cfg(Debug) >= 4) for(i=nb; i<=ne; i++) {
02104       if(!(spdvector[i].flag & OK)) continue;
02105       log << "GFE " << GDCUnique << " " << sat
02106          << " " << GDCUniqueFix
02107          << " " << time(i).printf(outFormat)
02108          << " " << setw(2) << spdvector[i].flag << fixed << setprecision(3);
02109       for(k=0; k<3; k++) log << " " << spdvector[i].data[L2]
02110             - (i < right->nbeg ? left->bias2-n1-(nadj+k-1) : right->bias2)
02111          << " " << PF[in[k]].Evaluate(spdvector[i].ndt - spdvector[nb].ndt);
02112       log << " " << setw(3) << spdvector[i].ndt << endl;
02113    }
02114 
02115    return nadj;
02116 }
02117 catch(Exception& e) { GPSTK_RETHROW(e); }
02118 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02119 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02120 }
02121 
02122 //------------------------------------------------------------------------------------
02123 //------------------------------------------------------------------------------------
02124 // fit a polynomial to the GF range, and change the units of -gfr(P2) and gfp(L2)
02125 // to cycles of wl21 (=5.4cm)
02126 int GDCPass::prepareGFdata(void) throw(Exception)
02127 {
02128 try {
02129    bool first;
02130    int i,nbeg,nend;
02131    unsigned ndeg;
02132 
02133    // decide on the degree of fit
02134    nbeg = SegList.begin()->nbeg;
02135    nend = SegList.begin()->nend;
02136    ndeg = 2 + int(0.5 + (nend-nbeg+1)*cfg(DT)/3000.0);
02137    if(ndeg > 6) ndeg = 6;
02138    //if(ndeg > int(cfg(GFPolyMaxDegree))) ndeg = int(cfg(GFPolyMaxDegree));
02139    if(ndeg < 2) ndeg = 2;
02140 
02141    // global fit to the gfr
02142    GFPassFit.Reset(ndeg);
02143 
02144    for(first=true,i=nbeg; i <= nend; i++) {
02145       if(!(spdvector[i].flag & OK)) continue;
02146 
02147       // 'change the bias' (initial bias only) in the GFP by changing units, also
02148       // slip fixing in the WL may have changed the values of GFP
02149       if(first) {
02150 //temp uncomment next 3 lines then comment again
02151          //if(fabs(spdvector[i].data[L2] - SegList.begin()->bias2) > 10.*wl21) {
02152          //   SegList.begin()->bias2 = spdvector[i].data[L2];
02153          //}
02154          SegList.begin()->bias2 /= wl21;
02155          first = false;
02156       }
02157 
02158       // 'change the arrays'
02159       // change units on the GFP and the GFR
02160       spdvector[i].data[P2] /= wl21;                    // -gfr (cycles of wl21)
02161       spdvector[i].data[L2] /= wl21;                    // gfp (cycles of wl21)
02162 
02163       // compute polynomial fit
02164       GFPassFit.Add(spdvector[i].data[P2],spdvector[i].ndt);
02165 
02166       // 'change the data'
02167       // save in L1                          // gfp+gfr residual (cycles of wl21)
02168       // ?? spdvector[i].data[L1] =
02169       //    spdvector[i].data[L2] - spdvector[i].data[P2] - SegList.begin()->bias2;
02170 // temp add -bias2  then remove it again
02171       spdvector[i].data[L1] = spdvector[i].data[L2] - spdvector[i].data[P2];
02172                                  // - SegList.begin()->bias2;
02173    }
02174 
02175    if(GFPassFit.isSingular()) {
02176       log << "Polynomial fit to GF range is singular! .. abort." << endl;
02177       return Singular;
02178    }
02179 
02180    return ReturnOK;
02181 }
02182 catch(Exception& e) { GPSTK_RETHROW(e); }
02183 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02184 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02185 }
02186 
02187 //------------------------------------------------------------------------------------
02188 //------------------------------------------------------------------------------------
02189 // detect slips in the geometry-free phase
02190 int GDCPass::detectGFslips(void) throw(Exception)
02191 {
02192 try {
02193    int i,iret;
02194    //temp double bias;
02195    list<Segment>::iterator it;
02196 
02197    // places first difference of GF in A1 - 'change the arrays' A1
02198    if( (iret = detectObviousSlips("GF"))) return iret;
02199 
02200    GFPassStats.Reset();
02201    //temp bias = 0.0;
02202    for(it=SegList.begin(); it != SegList.end(); it++) {
02203       // use bias for debiasing below
02204       // TD what if this segment deleted?
02205       //temp if(it == SegList.begin()) bias += it->bias2;
02206 
02207       // compute stats on dGF/dt
02208       for(i=it->nbeg; i <= it->nend; i++) {
02209          if(!(spdvector[i].flag & OK)) continue;
02210 
02211          // compute first-diff stats in meters
02212          // skip the first point in a segment - it is an obvious GF slip
02213          if(i > it->nbeg) GFPassStats.Add(spdvector[i].data[A1]*wl21);
02214 
02215          // if a gross GF slip was found, must remove bias in L1=GF(R-P)
02216          // in all subsequent segments ; 'change the data' L1
02217          //temp if(it != SegList.begin()) spdvector[i].data[L1] += bias - it->bias2;
02218 
02219       }  // end loop over data in segment it
02220 
02221       // TD delete segments if sigma too high?
02222 
02223       // check number of good points
02224       if(it->npts < int(cfg(MinPts))) {
02225          deleteSegment(it,"insufficient data in segment");
02226          continue;
02227       }
02228 
02229       // fit polynomial to GFR in each segment
02230       // compute (1stD of) fit residual GFP-fit(GFR) -> A1 - 'change the arrays' A1
02231       // delete segment if polynomial is singular - probably due to too little data
02232       if( (iret = GFphaseResiduals(it))) {
02233          //return iret;
02234          deleteSegment(it,"polynomial fit to GF residual failed");
02235          continue;
02236       }
02237    }
02238 
02239    // 'change the arrays'
02240    // at this point:
02241    // L1 = GFP+GFR in cycles, by prepareGFdata()
02242    // L2 = GFP in cycles, by prepareGFdata()
02243    // P1 = wlbias
02244    // P2 = GFR in cycles, by prepareGFdata()
02245    // A1 = GFP-(local fit) OR its 1stD, by GFphaseResiduals()
02246    //      (was 1stD of GFP+GFR (in L1), by firstDifferences())
02247    // A2 = 1stD of GFP (in L2), by firstDifferences()
02248    if( (iret = detectGFsmallSlips())) return iret;
02249 
02250    // delete all segments that are too small
02251    for(it=SegList.begin(); it != SegList.end(); it++) {
02252       if(it->npts < int(cfg(MinPts)))
02253          deleteSegment(it,"insufficient data in segment");
02254    }
02255 
02256    if(cfg(Debug) >= 4) dumpSegments("GFD",2,true);
02257 
02258    return ReturnOK;
02259 }
02260 catch(Exception& e) { GPSTK_RETHROW(e); }
02261 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02262 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02263 }
02264 
02265 //------------------------------------------------------------------------------------
02266 // for each segment, fit a polynomial to the gfr, then compute and store the
02267 // residual of fit
02268 int GDCPass::GFphaseResiduals(list<Segment>::iterator& it) throw(Exception)
02269 {
02270 try {
02271    int i,j,ndeg,nprev;
02272    double fit,rbias,prev,tmp;
02273    Stats<double> rofStats;
02274 
02275    // decide on the degree of fit
02276    ndeg = 2 + int(0.5 + (it->nend-it->nbeg+1)*cfg(DT)/3000.0);
02277    //if(ndeg > int(cfg(GFPolyMaxDegree))) ndeg = int(cfg(GFPolyMaxDegree));
02278    if(ndeg > 6) ndeg = 6;
02279    if(ndeg < 2) ndeg = 2;
02280 
02281    it->PF.Reset(ndeg);     // for fit to GF range
02282 
02283    for(i=it->nbeg; i <= it->nend; i++) {
02284       if(!(spdvector[i].flag & OK)) continue;
02285       it->PF.Add(spdvector[i].data[P2],spdvector[i].ndt);
02286    }
02287 
02288    if(it->PF.isSingular()) {     // this should never happen
02289       log << "Polynomial fit to GF range is singular in segment " << it->nseg
02290          << "! .. abort." << endl;
02291       return Singular;
02292    }
02293 
02294    // now compute the residual of fit
02295    rbias = prev = 0.0;
02296    rofStats.Reset();
02297    for(i=it->nbeg; i <= it->nend; i++) {
02298       // skip bad data
02299       if(!(spdvector[i].flag & OK)) continue;
02300 
02301       // TD? Use whole pass for small segments?
02302       //fit = GFPassFit.Evaluate(spdvector[i].ndt);  // use fit to gfr for whole pass
02303       fit = it->PF.Evaluate(spdvector[i].ndt);
02304 
02305       // all (fit, resid, gfr and gfp) are in cycles of wl21 (5.4cm)
02306 
02307       // compute gfp-(fit to gfr), store in A1 - 'change the arrays' A1 and A2
02308       // OR let's try first difference of residual of fit
02309       //           residual =  phase                            - fit to range
02310       spdvector[i].data[A1] = spdvector[i].data[L2] - it->bias2 - fit;
02311       if(rbias == 0.0) {
02312          rbias = spdvector[i].data[A1];
02313          nprev = spdvector[i].ndt - 1;
02314       }
02315       spdvector[i].data[A1] -= rbias;                    // debias residual for plots
02316 
02317          // compute stats on residual of fit
02318       rofStats.Add(spdvector[i].data[A1]);
02319 
02320       if(1) { // 1stD of residual - remember A1 has just been debiased
02321          tmp = spdvector[i].data[A1];
02322          spdvector[i].data[A1] -= prev;       // diff with previous epoch's
02323          // 040809 should this be divided by delta n?
02324          // spdvector[i].data[A1] /= (spdvector[i].ndt - nprev);
02325          prev = tmp;          // store residual for next point
02326          nprev = spdvector[i].ndt;
02327       }
02328 
02329       // store fit in A2
02330       //spdvector[i].data[A2] = fit;                   // fit to gfr (cycles of wl21)
02331       // store raw residual GFP-GFR (cycles of wl21) in A2
02332       //spdvector[i].data[A2]
02333       //    = spdvector[i].data[L2] - it->bias2 - spdvector[i].data[P2];
02334    }
02335 
02336    // TD? need this? use this?
02337    //log << "GFDsum " << GDCUnique << " " << sat << " " << it->nseg << " " << ndeg
02338    //   << " " << it->nbeg << " " << it->npts << " " << it->nend
02339    //   << " " << rofStats.N() << fixed << setprecision(3)
02340    //   << " " << rofStats.Minimum()
02341    //   << " " << rofStats.Maximum()
02342    //   << " " << rofStats.Average()
02343    //   << " " << rofStats.StdDev() << endl;
02344 
02345    return ReturnOK;
02346 }
02347 catch(Exception& e) { GPSTK_RETHROW(e); }
02348 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02349 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02350 }
02351 
02352 //------------------------------------------------------------------------------------
02353 // detect small slips in the geometry-free phase
02354 // TD outliers at the beginning or end of the segment....
02355 int GDCPass::detectGFsmallSlips(void) throw(Exception)
02356 {
02357 try {
02358    const int width=int(cfg(GFSlipWidth));
02359    int i,j,iplus,inew,ifirst,nok;
02360    list<Segment>::iterator it;
02361    Stats<double> pastStats,futureStats;
02362 
02363    // loop over segments
02364    for(it=SegList.begin(); it != SegList.end(); it++) {
02365 
02366       if(it->npts < 2*width+1) continue;
02367 
02368       // Cartoon of the GF 'two-pane moving window'
02369       //          point of interest:|
02370       // windows:     'past window' | 'future window'
02371       // stats  :        pastStats  | futureStats  (5 pts in each window)
02372       // data   : ... x (x x x x x) x (x x x x x) x ...
02373       //                 |          |          |
02374       // indexes:        j          i        iplus
02375 
02376       deque<int> pastIndex, futureIndex;
02377       pastStats.Reset();
02378       futureStats.Reset();
02379       i = inew = ifirst = -1;
02380       nok = 0;                          // recount
02381 
02382       // loop over points in the segment
02383       for(iplus=it->nbeg; iplus<=it->nend+width; iplus++) {
02384 
02385          // ignore bad points
02386          if(iplus <= it->nend && !(spdvector[iplus].flag & OK)) continue;
02387          if(ifirst == -1) ifirst = iplus;
02388 
02389          // pop the new i from the future
02390          if(futureIndex.size() == width || iplus > it->nend) {
02391             inew = futureIndex.front();
02392             futureIndex.pop_front();
02393             futureStats.Subtract(spdvector[inew].data[A1]);
02394             nok++;
02395          }
02396 
02397          // put iplus into the future deque
02398          if(iplus <= it->nend) {
02399             futureIndex.push_back(iplus);
02400             futureStats.Add(spdvector[iplus].data[A1]);
02401          }
02402          else
02403             futureIndex.push_back(-1);
02404 
02405          // check for outliers
02406          // we now have:
02407          //                (  past   )     ( future  )
02408          // data   : ... x (x x x x x) x x (x x x x x) x ...
02409          //                            | |          |
02410          // indexes:                   i inew     iplus
02411          // outlier if: (i,inew) = opposite signs but ~= large magnitude
02412          // if found, mark i bad and replace A1(inew) = A1(inew)+A1(i)
02413          if(foundGFoutlier(i,inew,pastStats,futureStats)) {
02414             // check that i was not marked a slip in the last iteration
02415             // if so, let inew be the slip and i the outlier
02416             if(spdvector[i].flag & DETECT) {
02417                //log << "Warning - marking a slip point BAD in GF detect small "
02418                //   << GDCUnique << " " << sat
02419                //   << " " << time(i).printf(outFormat) << " " << i << endl;
02420                spdvector[inew].flag = spdvector[i].flag;
02421                it->nbeg = inew;
02422             }
02423             spdvector[i].flag = BAD;
02424             spdvector[inew].data[A1] += spdvector[i].data[A1];
02425             learn["points deleted: GF outlier"]++;
02426             i = inew;
02427             nok--;
02428          }
02429 
02430          // pop last from past
02431          if(pastIndex.size() == width) {
02432             j = pastIndex.front();
02433             pastIndex.pop_front();
02434             pastStats.Subtract(spdvector[j].data[A1]);
02435          }
02436 
02437          // move i into the past
02438          if(i > -1) {
02439             pastIndex.push_back(i);
02440             pastStats.Add(spdvector[i].data[A1]);
02441          }
02442 
02443          // return to original state
02444          i = inew;
02445 
02446          // test for slip .. foundGF...prints to log
02447          if(foundGFsmallSlip(i,it->nseg,it->nend,it->nbeg,
02448             pastIndex,futureIndex,pastStats,futureStats)) {
02449 
02450             // create a new segment
02451             it->npts = nok-1;
02452             it = createSegment(it,i,"GF slip small");
02453             nok = 1;
02454 
02455             // mark it
02456             spdvector[i].flag |= GFDETECT;
02457 
02458             // TD print the "possible GF slip" and timetag here - see WLS
02459          }
02460 
02461       }  // end loop over points in the pass
02462       it->npts = nok;
02463 
02464    }  // end loop over segments
02465 
02466    return ReturnOK;
02467 }
02468 catch(Exception& e) { GPSTK_RETHROW(e); }
02469 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02470 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02471 }
02472 
02473 //------------------------------------------------------------------------------------
02474 bool GDCPass::foundGFoutlier(int i, int inew,
02475    Stats<double>& pastSt, Stats<double>& futureSt)
02476    throw(Exception)
02477 {
02478 try {
02479    if(i < 0 || inew < 0) return false;
02480    bool ok;
02481    double pmag = spdvector[i].data[A1]; // -pastSt.Average();
02482    double fmag = spdvector[inew].data[A1]; // -futureSt.Average();
02483    double var = ::sqrt(pastSt.Variance() + futureSt.Variance());
02484 
02485    ostringstream oss;
02486    if(cfg(Debug) >= 6) oss << "GFoutlier " << GDCUnique
02487       << " " << sat << " " << setw(3) << inew
02488       << " " << time(inew).printf(outFormat)
02489       << fixed << setprecision(3)
02490       << " p,fave=" << fabs(pmag) << "," << fabs(fmag)
02491       << " snr=" << fabs(pmag)/var <<","<< fabs(fmag)/var;
02492 
02493    bool isOut=true;
02494    for(;;) {
02495 
02496       // 1. signs must be opposite
02497       if(pmag * fmag >= 0) isOut=false;
02498       if(cfg(Debug) >= 6) oss << " (1)" << (isOut?"ok":"no");
02499       if(!isOut) break;
02500 
02501       //if(fabs(pmag+fmag) > 0.15*fabs(pmag-fmag))         // approx equal magnitude
02502          //{ if(cfg(Debug) >= 6) oss << endl; return false; }
02503 
02504       // 2. magnitudes must be large compared to noise
02505       double noise=cfg(GFSlipOutlier)*var;
02506       if(fabs(pmag) < noise || fabs(fmag) < noise) isOut=false;
02507       if(cfg(Debug) >= 6) oss << " (2)" << fabs(pmag)/var << "or" << fabs(fmag)/var
02508          << (isOut?">=":"<") << cfg(GFSlipOutlier);
02509       if(!isOut) break;
02510 
02511       if(cfg(Debug) >= 6) oss << " possible GF outlier";
02512 
02513       break;
02514    }  // end for(;;)
02515 
02516    if(cfg(Debug) >= 6) log << oss.str() << endl;
02517 
02518    return isOut;
02519 }
02520 catch(Exception& e) { GPSTK_RETHROW(e); }
02521 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02522 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02523 }
02524 
02525 //------------------------------------------------------------------------------------
02526 // Better to find too many small ones than to miss them, since the fixing algorithm
02527 // will most likely refuse to act on the questionable ones.
02528 bool GDCPass::foundGFsmallSlip(int i,int nseg,int iend,int ibeg,
02529    deque<int>& pastIn,deque<int>& futureIn,
02530    Stats<double>& pastSt, Stats<double>& futureSt)
02531    throw(Exception)
02532 {
02533 try {
02534    if(i < 0) return false;
02535 
02536    int j,k;
02537    double mag,pmag,fmag,pvar,fvar;
02538 
02539    pmag = fmag = pvar = fvar = 0.0;
02540    // note when past.N == 1, this is first good point, which has 1stD==0
02541    // TD be very careful when N is small
02542    if(pastSt.N() > 0) pmag = spdvector[i].data[A1]-pastSt.Average();
02543    if(futureSt.N() > 0) fmag = spdvector[i].data[A1]-futureSt.Average();
02544    if(pastSt.N() > 1) pvar = pastSt.Variance();
02545    if(futureSt.N() > 1) fvar = futureSt.Variance();
02546    mag = (pmag + fmag) / 2.0;
02547 
02548    if(cfg(Debug) >= 6) log << "GFS " << GDCUnique
02549       << " " << sat << " " << nseg
02550       << " " << time(i).printf(outFormat)
02551       //<< " P( " << setw(3) << pastIn[0]            // don't print this...
02552       //<< " " << setw(3) << pastIn[1]
02553       //<< " " << setw(3) << pastIn[2]
02554       //<< " " << setw(3) << pastIn[3]
02555       //<< " " << setw(3) << pastIn[4]
02556       //<< ") " << setw(3) << i
02557       //<< " F( " << setw(3) << futureIn[0]
02558       //<< " " << setw(3) << futureIn[1]
02559       //<< " " << setw(3) << futureIn[2]
02560       //<< " " << setw(3) << futureIn[3]
02561       //<< " " << setw(3) << futureIn[4] << ")"       // ...to here
02562       << fixed << setprecision(3)
02563       << " " << setw(3) << pastSt.N()
02564       << " " << setw(7) << pastSt.Average()
02565       << " " << setw(7) << pastSt.StdDev()
02566       << " " << setw(3) << futureSt.N()
02567       << " " << setw(7) << futureSt.Average()
02568       << " " << setw(7) << futureSt.StdDev()
02569       << " " << setw(7) << mag
02570       << " " << setw(7) << ::sqrt(pvar+fvar)
02571       << " " << setw(9) << spdvector[i].data[A1]
02572       << " " << setw(7) << pmag
02573       << " " << setw(7) << pvar
02574       << " " << setw(7) << fmag
02575       << " " << setw(7) << fvar
02576       << " " << setw(3) << i
02577       << endl;
02578 
02579    //                    x                    -- mag
02580    //
02581    //    x   x   x   x                         - step
02582    //                       x    x   x   x   ---
02583    const double minMag=cfg(GFSlipSize);     // minimum slip magnitude
02584    const double STN=cfg(GFSlipStepToNoise); // step (past->future) to noise ratio
02585    const double MTS=cfg(GFSlipToStep);      // magnitude to step ratio
02586    const double MTN=cfg(GFSlipToNoise);     // magnitude to noise ratio
02587    const int Edge=int(cfg(GFSlipEdge));     // number of points before edge
02588    const double RangeCheckLimit = 2*cfg(WLSigma)/(0.83*wl21);
02589                                                 // 2 * range noise in units of wl21
02590    const double snr=fabs(pmag-fmag)/::sqrt(pvar+fvar);
02591 
02592    // 050109 if Debug=6, print only possible slips, if 7 print all
02593    bool isSlip=true;
02594    ostringstream oss;
02595 
02596    for(;;) {
02597       // NB last printed test is a failure unless it says possible GF slip
02598       if(cfg(Debug) >= 6) oss << "GFslip " << GDCUnique
02599          << " " << sat << " " << nseg
02600          << " " << setw(3) << i
02601          << " " << time(i).printf(outFormat) << fixed << setprecision(3)
02602          << " mag=" << mag << " snr=" << snr;      // no endl
02603 
02604       // 0. if WL slip here - ...?
02605 
02606       // 1. slip must be non-trivial
02607       if(fabs(mag) <= minMag) isSlip=false;
02608       if(cfg(Debug) >= 6) oss << " (1)" << fabs(mag) << (isSlip?">":"<=") << minMag;
02609       if(!isSlip) break;
02610 
02611       // 2. change in average is larger than noise
02612       if(snr <= STN) isSlip=false;
02613       if(cfg(Debug) >= 6) oss << " (2)" << snr << (isSlip?">":"<=") << STN;
02614       if(!isSlip) break;
02615 
02616       // 3. slip is large compared to change in average
02617       if(fabs(mag) <= MTS*fabs(pmag-fmag)) isSlip=false;
02618       if(cfg(Debug) >= 6) oss << " (3)" << fabs(mag/(pmag-fmag))
02619          << (isSlip?">":"<=") << MTS;
02620       if(!isSlip) break;
02621 
02622       // 4. magnitude is large compared to noise: a 3-sigma slip
02623       if(fabs(mag) <= MTN*::sqrt(pvar+fvar)) isSlip=false;
02624       if(cfg(Debug) >= 6) oss << " (4)" << fabs(mag)/::sqrt(pvar+fvar)
02625          << (isSlip?">":"<=") << MTN;
02626       if(!isSlip) break;
02627 
02628       // 5. if very close to edge, declare it an outlier
02629       if(pastSt.N() < Edge || futureSt.N() < Edge+1) isSlip=false;
02630       if(cfg(Debug) >= 6) oss << " (5)" << pastSt.N() << "," << futureSt.N()
02631          << (isSlip?">":"<") << Edge;
02632       if(!isSlip) break;
02633 
02634       // 5.5 TD? if slip is within a few epochs of WL slip - skip it
02635 
02636       // 6. large slips (compared to range noise): check the GFR-GFP for consistency
02637       if(fabs(mag) > RangeCheckLimit) {
02638          double magGFR,mtnGFR;
02639          Stats<double> pGFRmPh,fGFRmPh;
02640          for(j=0; j<pastIn.size(); j++) {
02641             if(pastIn[j] > -1) pGFRmPh.Add(spdvector[pastIn[j]].data[L1]);
02642             if(futureIn[j] > -1) fGFRmPh.Add(spdvector[futureIn[j]].data[L1]);
02643          }
02644          magGFR = spdvector[i].data[L1] - (pGFRmPh.Average()+fGFRmPh.Average())/2.0;
02645          mtnGFR = fabs(magGFR)/::sqrt(pGFRmPh.Variance()+fGFRmPh.Variance());
02646 
02647          if(cfg(Debug) >= 6)
02648             oss << "; GFR-GFP has mag: " << magGFR
02649                << ", |dmag|: " << fabs(mag-magGFR)
02650                << " and mag/noise " << mtnGFR;
02651 
02652          // TD test - mag must ~= magGFR if magGFR/noiseGFR >> 1
02653          // test - metz 54,56,57,58
02654          if(fabs(mag-magGFR) > fabs(magGFR)) isSlip=false;
02655          if(cfg(Debug) >= 6) oss << " (6a)" << fabs(mag-magGFR)
02656             << (isSlip?"<=":">") << fabs(magGFR);
02657          if(!isSlip) break;
02658 
02659          if(mtnGFR < 3) isSlip=false;
02660          if(cfg(Debug) >= 6) oss << " (6b)"
02661             << mtnGFR << "><3:can" << (isSlip?"":"not") << "_see_in_GFR";
02662          if(!isSlip) break;
02663       }
02664 
02665       // 7. small slips (compared to variations in dGF): extra careful
02666       // TD beware of small slips in the presence of noise >~ 1
02667       else { //if(fabs(mag) <= RangeCheckLimit)
02668          double magFD;
02669          Stats<double> fdStats;
02670          j = i-1; k=0;
02671          while(j >= ibeg && k < 15) {
02672             if(spdvector[j].flag & OK) { fdStats.Add(spdvector[j].data[A2]); k++; }
02673             j--;
02674          }
02675          j = i+1; k=0;
02676          while(j <= iend && k < 15) {
02677             if(spdvector[j].flag & OK) { fdStats.Add(spdvector[j].data[A2]); k++; }
02678             j++;
02679          }
02680          magFD = spdvector[i].data[A2] - fdStats.Average();
02681 
02682          if(cfg(Debug) >= 6)
02683             oss << " (7)1stD(GFP)mag=" << magFD
02684                << ",noise=" << fdStats.StdDev()
02685                << ",snr=" << fabs(magFD)/fdStats.StdDev()
02686                << ",maxima=" << fdStats.Minimum() << "," << fdStats.Maximum();
02687       }
02688 
02689       break;
02690    }  // end for(;;)
02691 
02692    if(isSlip) oss << " possible GF slip";
02693    if(cfg(Debug) >= 6) log << oss.str() << endl;
02694 
02695    return isSlip;
02696 }
02697 catch(Exception& e) { GPSTK_RETHROW(e); }
02698 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02699 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02700 }
02701 
02702 //------------------------------------------------------------------------------------
02703 //------------------------------------------------------------------------------------
02704 // check the consistency of WL slips where a GF slip, but not a WL slip, was detected.
02705 int GDCPass::WLconsistencyCheck(void) throw(Exception)
02706 {
02707 try {
02708    int i,k;
02709    const int N=2*int(cfg(WLWindowWidth));
02710    double mag,absmag,factor=wl2/wl21;
02711 
02712    // loop over the data and look for points with GFDETECT but not WLDETECT or WLFIX
02713    for(i=0; i<size(); i++) {
02714 
02715       if(!(spdvector[i].flag & OK)) continue;        // bad
02716       if(!(spdvector[i].flag & DETECT)) continue;    // no slips
02717       if(spdvector[i].flag & WLDETECT) continue;     // WL was detected
02718 
02719       // GF only slip - compute WL stats on both sides
02720       Stats<double> futureStats,pastStats;
02721       k = i;
02722       // fill future
02723       while(k < size() && futureStats.N() < N) {
02724          if(spdvector[k].flag & OK)                  // data is good
02725             futureStats.Add(spdvector[k].data[P1]);        // wlbias
02726          k++;
02727       }
02728       // fill past
02729       k = i-1;
02730       while(k >= 0 && pastStats.N() < N) {
02731          if(spdvector[k].flag & OK)                  // data is good
02732             pastStats.Add(spdvector[k].data[P1]);          // wlbias
02733          k--;
02734       }
02735 
02736       // is there a WL slip here?
02737       // 1. magnitude of slip > 0.75
02738       // 2. magnitude is > stddev on both sides
02739       // 3. N() > 10 on both sides TD??
02740       mag = futureStats.Average()-pastStats.Average();
02741       absmag = fabs(mag);
02742 
02743       if(absmag > cfg(WLSlipSize) &&      // 0.75 &&
02744          absmag > pastStats.StdDev() &&
02745          absmag > futureStats.StdDev()) {
02746 
02747          long nwl;
02748          nwl = long(mag + (mag > 0 ? 0.5 : -0.5));
02749 
02750          if(nwl == 0) continue;
02751 
02752          // now do the fixing - change the data to the future of the slip
02753          for(k=i; k<size(); k++) {
02754             //if(!(spdvector[i].flag & OK)) continue;
02755             spdvector[k].data[P1] -= nwl;                                 // WLbias
02756             spdvector[k].data[L2] -= nwl * factor;                        // GFP
02757          }
02758 
02759          // Add to slip list
02760          Slip newSlip(i);
02761          newSlip.NWL = nwl;
02762          newSlip.msg = "WL";
02763          SlipList.push_back(newSlip);
02764 
02765          // mark it
02766          spdvector[i].flag |= (WLDETECT + WLFIX);
02767 
02768          if(cfg(Debug) >= 7) log << "CHECK " << GDCUnique << " " << sat
02769             << " " << i
02770             << " " << time(i).printf(outFormat)
02771             << fixed << setprecision(3)
02772             << "  " << pastStats.N()
02773             //<< " " << pastStats.Average()
02774             << " " << pastStats.StdDev()
02775             << "  " << futureStats.N()
02776             //<< " " << futureStats.Average()
02777             << " " << futureStats.StdDev()
02778             << "  " << futureStats.Average() - pastStats.Average()
02779             << " " << nwl
02780             << endl;
02781 
02782       }
02783    }
02784 
02785    return ReturnOK;
02786 }
02787 catch(Exception& e) { GPSTK_RETHROW(e); }
02788 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
02789 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
02790 }
02791 
02792 //------------------------------------------------------------------------------------
02793 //------------------------------------------------------------------------------------
02794 // last call before returning:
02795 // generate editing commands for deleted (flagged) data,
02796 // use editing command (slips and deletes) to modify the original SatPass data
02797 // print ending summary, and also return it as a string
02798 string GDCPass::finish(int iret, SatPass& svp, vector<string>& editCmds)
02799    throw(Exception)
02800 {
02801 try {
02802    bool ok;
02803    int i,ifirst,ilast,npts;
02804    long N1,N2,prevN1,prevN2;
02805    double slipL1,slipL2,WLbias,GFbias;
02806    //SatPassData spd;
02807    list<Slip>::iterator jt;
02808    string retMessage;
02809 
02810    // ---------------------------------------------------------
02811    // sort the slips in time
02812    SlipList.sort();
02813 
02814    // ---------------------------------------------------------
02815    // merge *this GDCPass and the input SatPass...
02816    // use this->flag to generate edit commands for data marked bad
02817    // use the SlipList to fix slips
02818    // 'change the arrays' A1 and A2 - fill with range minus phase for output
02819    npts = 0;
02820    ilast = -1;                         // ilast is the index of the last good point
02821    ifirst = -1;                        // ifirst is the index of the first good point
02822    WLbias = GFbias = slipL1 = slipL2 = 0.0;
02823    prevN1 = prevN2 = 0L;
02824    jt = SlipList.begin();
02825    for(i=0; i<size(); i++) {
02826 
02827       // is this point bad?
02828       if(!(spdvector[i].flag & OK)) {       // data is bad
02829          ok = false;
02830          if(i == size() - 1) {    // but this is the last point
02831             i++;
02832             ok = true;
02833          }
02834       }
02835       else ok = true;                  // data is good
02836 
02837       if(ok) {
02838          if(ifirst == -1) ifirst = i;
02839 
02840          // generate edit commands: delete from ilast+1 to i-1
02841          if(i-ilast > 2 && cfg(OutputDeletes) != 0) {
02842             // delete 2, or a range of, points
02843             // -DS+<sat>,<time>
02844             ostringstream stst1;
02845             stst1 << "-DS";
02846             if(i-ilast > 3) stst1 << "+";
02847             stst1 << sat << ",";
02848             if(cfg(OutputGPSTime))
02849                stst1 << time(ilast+1).printf("%F,%.3g");
02850             else
02851                stst1 << time(ilast+1).printf("%Y,%m,%d,%H,%M,%f");
02852             if(i-ilast > 3) stst1 << " # begin delete of "
02853                   << asString(i+1-ilast) << " points";
02854             editCmds.push_back(stst1.str());
02855 
02856             // -DS-<sat>,<time>
02857             ostringstream stst2;
02858             stst2 << "-DS";
02859             if(i-ilast > 3) stst2 << "-";
02860             stst2 << sat << ",";
02861             if(cfg(OutputGPSTime))
02862                stst2 << time(i-1).printf("%F,%.3g");
02863             else
02864                stst2 << time(i-1).printf("%Y,%m,%d,%H,%M,%f");
02865             if(i-ilast > 3) stst2 << " # end delete of "
02866                << asString(i+1-ilast) << " points";
02867             editCmds.push_back(stst2.str());
02868          }
02869          else if(i-ilast > 1 && cfg(OutputDeletes) != 0) {
02870             // delete a single isolated point
02871             ostringstream stst;
02872             stst << "-DS" << sat << ",";
02873             if(cfg(OutputGPSTime))
02874                stst << time(i-1).printf("%F,%.3g");
02875             else
02876                stst << time(i-1).printf("%Y,%m,%d,%H,%M,%f");
02877             editCmds.push_back(stst.str());
02878          }
02879 
02880          ilast = i;
02881          npts++;
02882       }
02883 
02884       // keep track of net slip fix
02885       if(jt != SlipList.end() && i == jt->index) {          // there is a slip here
02886          // fix the slip by changing the bias added to phase
02887          N1 = jt->N1;
02888          N2 = jt->N1 - jt->NWL;
02889          slipL1 += double(N1);
02890          slipL2 += double(N2);
02891 
02892          // generate edit commands
02893          {
02894             ostringstream stst;
02895             stst << "-BD+" << sat << ",L1,";
02896             if(cfg(OutputGPSTime))
02897                stst << time(jt->index).printf("%F,%.3g");
02898             else
02899                stst << time(jt->index).printf("%Y,%m,%d,%H,%M,%f");
02900             stst << "," << N1-prevN1;
02901             if(!jt->msg.empty()) stst << " # " << jt->msg;
02902             //stst << " # WL: " << jt->NWL << " N1: " << jt->N1; //temp
02903             editCmds.push_back(stst.str());
02904          }
02905          {
02906             ostringstream stst;
02907             stst << "-BD+" << sat << ",L2,";
02908             if(cfg(OutputGPSTime))
02909                stst << time(jt->index).printf("%F,%.3g");
02910             else
02911                stst << time(jt->index).printf("%Y,%m,%d,%H,%M,%f");
02912             stst << "," << N2-prevN2;
02913             if(!jt->msg.empty()) stst << " # " << jt->msg;
02914             editCmds.push_back(stst.str());
02915          }
02916 
02917          prevN1 = N1;
02918          prevN2 = N2;
02919          jt++;
02920       }
02921 
02922       if(i >= size()) break;
02923 
02924       // 'change the data' for the last time
02925       spdvector[i].data[L1] = svp.data(i,DCobstypes[L1]) - slipL1;
02926       spdvector[i].data[L2] = svp.data(i,DCobstypes[L2]) - slipL2;
02927       spdvector[i].data[P1] = svp.data(i,DCobstypes[P1]);
02928       spdvector[i].data[P2] = svp.data(i,DCobstypes[P2]);
02929 
02930       // compute range minus phase for output
02931       // do the same at the beginning ("BEG")
02932 
02933       // compute WL and GFP
02934          // narrow lane range (m)
02935       double wlr = wl1r * spdvector[i].data[P1] + wl2r * spdvector[i].data[P2];
02936          // wide lane phase (m)
02937       double wlp = wl1p * spdvector[i].data[L1] + wl2p * spdvector[i].data[L2];
02938          // geo-free range (m)
02939       double gfr = gf1r * spdvector[i].data[P1] + gf2r * spdvector[i].data[P2];
02940          // geo-free phase (m)
02941       double gfp = gf1p * spdvector[i].data[L1] + gf2p * spdvector[i].data[L2];
02942       if(i == ifirst) {
02943          WLbias = (wlp-wlr)/wlwl;
02944          GFbias = gfp;
02945       }
02946       spdvector[i].data[A1] = (wlp-wlr)/wlwl - WLbias; // wide lane bias (cyc)
02947       spdvector[i].data[A2] = gfp - GFbias;            // geo-free phase (m)
02948       //spdvector[i].data[A2] = gfr - gfp;             // geo-free range - phase (m)
02949 
02950    } // end loop over all data
02951 
02952    // first fix the segment for dump - TD? is this necessary?
02953    if(SegList.begin() != SegList.end()) {
02954       SegList.begin()->bias1 = SegList.begin()->bias2 = 0;     // not necessary..
02955       SegList.begin()->nbeg = 0;
02956       SegList.begin()->nend = size()-1;
02957       SegList.begin()->npts = npts;
02958    }
02959    // dump the corrected data
02960    if(cfg(Debug) >= 2) dumpSegments("AFT",2,true);
02961 
02962    // dump the edit commands to log
02963    for(i=0; i<editCmds.size(); i++)
02964       if(cfg(Debug) >= 2)
02965          log << "EditCmd: " << GDCUnique << " " << editCmds[i] << endl;
02966 
02967    // ---------------------------------------------------------
02968    // copy corrected data into original SatPass, without disturbing other obs types
02969    for(i=0; i<size(); i++) {
02970       svp.data(i,DCobstypes[L1]) = spdvector[i].data[L1];
02971       svp.data(i,DCobstypes[L2]) = spdvector[i].data[L2];
02972       svp.data(i,DCobstypes[P1]) = spdvector[i].data[P1];
02973       svp.data(i,DCobstypes[P2]) = spdvector[i].data[P2];
02974 
02975       // change the flag for use by SatPass
02976       //const unsigned short SatPass::OK  = 1; good data
02977       //const unsigned short SatPass::BAD = 0; used by caller and DC to mark bad data
02978       //const unsigned short SatPass::LL1 = 2; discontinuity on L1 only
02979       //const unsigned short SatPass::LL2 = 4; discontinuity on L2 only
02980       //const unsigned short SatPass::LL3 = 6; discontinuity on L1 and L2
02981       //const unsigned short GDCPass::DETECT   =   6;  // = WLDETECT | GFDETECT
02982       //const unsigned short GDCPass::FIX      =  24;  // = WLFIX | GFFIX
02983       if(spdvector[i].flag & OK) {
02984 //??     if(((spdvector[i].flag & DETECT)!=0 && (spdvector[i].flag & FIX)==0)
02985          if(((spdvector[i].flag & DETECT)==0 && (spdvector[i].flag & FIX)!=0)
02986             || i == ifirst)
02987             spdvector[i].flag = LL3 + OK;
02988          else
02989             spdvector[i].flag = OK;
02990       }
02991       else
02992          spdvector[i].flag = BAD;
02993 
02994       svp.LLI(i,DCobstypes[L1]) = (spdvector[i].flag & LL1) ? 1 : 0;
02995       svp.LLI(i,DCobstypes[L2]) = (spdvector[i].flag & LL2) ? 1 : 0;
02996       svp.setFlag(i,spdvector[i].flag);
02997    }
02998 
02999    // ---------------------------------------------------------
03000    // print stuff at the end
03001    if(cfg(Debug) >= 1) retMessage = dumpSegments("GDC",1);  // dump also prints it
03002 
03003    ostringstream oss;
03004    // print WL & GF stats for whole pass
03005    if(cfg(Debug) > 0 && WLPassStats.N() > 2) {
03006       oss << "GDC " << GDCUnique << " " << sat
03007          << " " << fixed << setprecision(3) << WLPassStats.StdDev()
03008          << " WL sigma in cycles"
03009          << " N=" << WLPassStats.N()
03010          << " Min=" << WLPassStats.Minimum()
03011          << " Max=" << WLPassStats.Maximum()
03012          << " Ave=" << WLPassStats.Average();
03013       if(WLPassStats.StdDev() > cfg(WLSigma))
03014          oss << " Warning - WL sigma > input (" << cfg(WLSigma) << ")";
03015       oss << endl;
03016    }
03017 
03018    if(cfg(Debug) > 0 && GFPassStats.N() > 2) {
03019       oss << "GDC " << GDCUnique << " " << sat
03020          << " " << fixed << setprecision(3) << GFPassStats.StdDev()
03021          << " sigma GF variation in meters per DT"
03022          << " N=" << GFPassStats.N()
03023          << " Min=" << GFPassStats.Minimum()
03024          << " Max=" << GFPassStats.Maximum()
03025          << " Ave=" << GFPassStats.Average()
03026          << endl;
03027       oss << "GDC " << GDCUnique << " " << sat
03028          << " " << fixed << setprecision(3)
03029          << (fabs(GFPassStats.Minimum()) > fabs(GFPassStats.Maximum()) ?
03030             fabs(GFPassStats.Minimum()) : fabs(GFPassStats.Maximum()))
03031          << " maximum GF variation in meters per DT"
03032          << " N=" << GFPassStats.N()
03033          << " Ave=" << GFPassStats.Average()
03034          << " Std=" << GFPassStats.StdDev()
03035          << endl;
03036    }
03037 
03038    // print 'learn' summary
03039    if(cfg(Debug) > 0) {
03040       map<string,int>::const_iterator kt;
03041       for(kt=learn.begin(); kt != learn.end(); kt++)
03042          oss << "GDC " << GDCUnique << " " << sat
03043             << " " << setw(3) << kt->second << " " << kt->first << endl;
03044 
03045       int n = int((lastTime-firstTime)/cfg(DT)) + 1;
03046       double percent = 100.0*double(ngood)/n;
03047       if(cfg(Debug) > 0) oss << "GDC# " << setw(2) << GDCUnique << ", SAT " << sat
03048          << ", Pts: " << setw(4) << n << " total " << setw(4) << ngood
03049          << " good " << fixed << setprecision(1) << setw(5) << percent << "%"
03050          << ", start " << firstTime.printf(outFormat)
03051          << endl;
03052    }
03053 
03054    if(iret) {
03055       oss << "GDC " << setw(3) << GDCUnique << " " << sat
03056          << " " << firstTime.printf(outFormat)
03057          << " is returning with error code: "
03058          << (iret == NoData ? "insufficient data" :
03059             (iret == Singular ? "singularity" :
03060             (iret == FatalProblem ? "fatal problem" : "unknown problem")
03061             //(iret == PrematureEnd ? "premature end" : "unknown problem")
03062             ))
03063          << endl;
03064    }
03065 
03066    retMessage += oss.str();
03067    if(oss.str().size() > 0) log << oss.str();
03068 
03069    if(cfg(Debug) >= 2) log << "======== End GPSTK Discontinuity Corrector "
03070                            << GDCUnique << " ================================================" << endl;
03071 
03072    return retMessage;
03073 }
03074 catch(Exception& e) { GPSTK_RETHROW(e); }
03075 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
03076 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
03077 }
03078 
03079 //------------------------------------------------------------------------------------
03080 // create, delete and dump Segments
03081 //------------------------------------------------------------------------------------
03082 // create a new segment from the given one, starting at index ibeg,
03083 // and insert it after the given iterator.
03084 // Return an iterator pointing to the new segment. String msg is for debug output
03085 list<Segment>::iterator GDCPass::createSegment(list<Segment>::iterator sit,
03086                                                int ibeg, string msg) throw(Exception)
03087 {
03088 try {
03089    Segment s;
03090    s = *sit;
03091    s.nbeg = ibeg;
03092    s.nend = sit->nend;
03093    sit->nend = ibeg-1;
03094 
03095    // 'trim' beg and end indexes
03096    while(s.nend > s.nbeg && !(spdvector[s.nend].flag & OK)) s.nend--;
03097    while(sit->nend > sit->nbeg && !(spdvector[sit->nend].flag & OK)) sit->nend--;
03098 
03099    // get the segment number right
03100    s.nseg++;
03101    list<Segment>::iterator skt=sit;
03102    for(skt++; skt != SegList.end(); skt++) skt->nseg++;
03103 
03104    if(cfg(Debug) >= 6)
03105       log << "SEG " << GDCUnique << " " << sat
03106          << " " << msg
03107          << " " << time(ibeg).printf(outFormat)
03108          << " " << s.nbeg << " - " << s.nend
03109          << " biases " << fixed << setprecision(3) << s.bias1 << " " << s.bias2
03110          << endl;
03111 
03112    learn["breaks found: " + msg]++;
03113 
03114    return SegList.insert(++sit,s); // insert puts s before ++sit
03115 }
03116 catch(Exception& e) { GPSTK_RETHROW(e); }
03117 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
03118 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
03119 }
03120 
03121 //------------------------------------------------------------------------------------
03122 // dump a list of the segments
03123 // level=0 one line summary (number of segments)
03124 // level=1 one per line list of segments
03125 // level=2 dump all data, including (if extra) temporary arrays
03126 // return level 1 output as string
03127 string GDCPass::dumpSegments(string label, int level, bool extra) throw(Exception)
03128 {
03129 try {
03130    int i,ifirst,ilast;
03131    list<Segment>::iterator it;
03132    string msg;
03133    ostringstream oss;
03134 
03135       // summary of SegList
03136    oss << label << " " << GDCUnique
03137       << " list of Segments (" << SegList.size() << "):"
03138       << endl;
03139 
03140    if(level < 1) { msg = oss.str(); log << msg; return msg; }
03141 
03142       // one line per segment
03143    ilast = -1;                               // last good point
03144    for(it=SegList.begin(); it != SegList.end(); it++) {
03145       //if(it->npts > 0) {
03146       //   biaswl = spdvector[it->nbeg].data[P1];
03147       //   biasgf = spdvector[it->nbeg].data[L2];
03148       //}
03149       //else biaswl = biasgf = 0.0;
03150 
03151       i = (it->nend - it->nbeg + 1);         // total number of points
03152 
03153       oss << label << " " << GDCUnique << " " << sat
03154          << " #" << setw(2) << it->nseg << ": "
03155          << setw(4) << it->npts << "/" << setw(4) << i << " pts, # "
03156          << setw(4) << it->nbeg << "-" << setw(4) << it->nend
03157          << " (" << time(it->nbeg).printf(outFormat)
03158          << " - " << time(it->nend).printf(outFormat)
03159          << ")";
03160 
03161       if(it->npts > 0) {
03162          oss << fixed << setprecision(3)
03163             << " bias(wl)=" << setw(13) << it->bias1 //biaswl
03164             << " bias(gf)=" << setw(13) << it->bias2; //biasgf;
03165          if(ilast > -1) {
03166             ifirst = it->nbeg;
03167             while(ifirst <= it->nend && !(spdvector[ifirst].flag & OK)) ifirst++;
03168             i = spdvector[ifirst].ndt - spdvector[ilast].ndt;
03169             oss << " Gap " << setprecision(1) << setw(5)
03170                << cfg(DT)*i << " s = " << i << " pts.";
03171          }
03172          ilast = it->nend;
03173          while(ilast >= it->nbeg && !(spdvector[ilast].flag & OK)) ilast--;
03174       }
03175 
03176       oss << endl;
03177    }
03178 
03179    msg = oss.str();
03180    log << msg;
03181    if(level < 2) return msg;
03182 
03183       // dump the data
03184    for(it=SegList.begin(); it != SegList.end(); it++) {
03185       for(i=it->nbeg; i<=it->nend; i++) {
03186          //if(!(spdvector[i].flag & OK)) continue;  //dfplot ignores bad data
03187 
03188          log << "DSC" << label << " " << GDCUnique << " " << sat << " " << it->nseg
03189             << " " << time(i).printf(outFormat)
03190             << " " << setw(3) << spdvector[i].flag
03191             << fixed << setprecision(3)
03192             << " " << setw(13) << spdvector[i].data[L1] - it->bias2 //biasgf  //temp
03193             << " " << setw(13) << spdvector[i].data[L2] - it->bias2 //biasgf
03194             << " " << setw(13) << spdvector[i].data[P1] - it->bias1 //biaswl
03195             << " " << setw(13) << spdvector[i].data[P2];
03196          if(extra) log
03197             << " " << setw(13) << spdvector[i].data[A1]
03198             << " " << setw(13) << spdvector[i].data[A2];
03199          log << " " << setw(4) << i;          // TD? make this spdvector[i].ndt?
03200          if(i == it->nbeg) log
03201             << " " << setw(13) << it->bias1 //biaswl
03202             << " " << setw(13) << it->bias2; //biasgf;
03203          log << endl;
03204       }
03205    }
03206    return msg;
03207 }
03208 catch(Exception& e) { GPSTK_RETHROW(e); }
03209 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
03210 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
03211 }
03212 
03213 //------------------------------------------------------------------------------------
03214 void GDCPass::deleteSegment(list<Segment>::iterator& it, string msg) throw(Exception)
03215 {
03216 try {
03217    int i;
03218 
03219    if(cfg(Debug) >= 6) log << "Delete segment "
03220       << GDCUnique << " " << sat << " " << it->nseg
03221       << " pts " << it->npts
03222       << " indexes " << it->nbeg << " - " << it->nend
03223       << " start " << firstTime.printf(outFormat)
03224       << " : " << msg
03225       << endl;
03226 
03227    it->npts = 0;
03228    for(i=it->nbeg; i<=it->nend; i++) if(spdvector[i].flag & OK) {
03229       // count these : learn
03230       learn["points deleted: " + msg]++;
03231       spdvector[i].flag = BAD;
03232    }
03233 
03234    learn["segments deleted: " + msg]++;
03235 }
03236 catch(Exception& e) { GPSTK_RETHROW(e); }
03237 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
03238 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
03239 }
03240 
03241 //------------------------------------------------------------------------------------
03242 //------------------------------------------------------------------------------------

Generated on Wed Feb 8 03:30:58 2012 for GPS ToolKit Software Library by  doxygen 1.3.9.1