doubledifference.cpp

gpstk1.5的有关内容 对于刚刚接触gps有一定的帮助 很有用的啊

               sddata.L1.push_back(it1->second.L1[i] - it2->second.L1[j]);
               sddata.L2.push_back(it1->second.L2[i] - it2->second.L2[j]);
               sddata.P1.push_back(it1->second.P1[i] - it2->second.P1[j]);
               sddata.P2.push_back(it1->second.P2[i] - it2->second.P2[j]);
               sddata.ER.push_back(it1->second.ER[i] - it2->second.ER[j]);
               sddata.elev.push_back(it1->second.elev[i]);

            }  // end if elevation is ok

               // next epoch
            i++;
            j++;
         }
            // i is behind j in time(count)
         else if(it1->second.count[i] < it2->second.count[j])
            i++;
            // i is ahead of j in time(count)
         else
            j++;

      }  // end while

         // save it in the map
      SDmap[sdid] = sddata;

   }  // end loop over satellites at first site

}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
// Assume SDmap is all for the same baseline
int ComputeDoubleDifferences(map<SDid,RawData>& SDmap)
{
try {
   bool frst,ok;
   int i,j,k,indx,count,ddsign;
   long nn1,nn2;
   double ddL1,ddL2,ddER,ddP1,ddP2,dd,db1,db2;
   DayTime tt,ttnext;   // ttnext is the time of the next reference satellite switch
   //SDid sid,ref;        // SDid of the current satellite and reference satellite
   map<SDid,int> Inext; // index in count (all) buffers which is to be processed next
   map<SDid,RawData>::const_iterator it;

   if(SDmap.size() == 0) return 0;

      // initialize the 'next index' map to zero
      // find the smallest (earliest) count
   frst = true;
   for(it=SDmap.begin(); it != SDmap.end(); it++) {
      int jj=0;
      Inext[it->first] = jj;  
      if(frst || it->second.count[0] < count) {
         count = it->second.count[0];
         frst = false;
      }
   }
   
      // ref will be the SDid for the reference satellite
   SDid ref = SDmap.begin()->first;        // ref.sat is TBD by timetable

      // loop over epochs in the SDs
   ttnext = DayTime::BEGINNING_OF_TIME;
   while(1) {
         // time at this count
      tt = FirstEpoch + count * CI.DataInterval;

         // get the reference satellite at this time
      if(tt > ttnext) {
         ttnext = tt;
         if(QueryTimeTable(ref, ttnext)) {         // error - timetable failed
            oflog << "DD: Error - failed to find reference from timetable at "
               << tt.printf("%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << " count "
               << count << " for baseline " << ref.site1 << "-" << ref.site2 << endl;
            return 1;
         }
         if(CI.Verbose) oflog << "DD: reference is set to " << ref << " at "
            << tt.printf("%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g")
            << " count " << count << endl;
      }

         // does reference satellite have data at this count?
      if(SDmap[ref].count[Inext[ref]] != count) {
         oflog << "Error - failed to find reference data " << ref << " at "
            << tt.printf("%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g") << endl;
            // TD return here, or just skip the epoch?
            // question is do we allow 'holes' in ref sat's data?
         return 1;
      }

         // compute DDs
      for(it=SDmap.begin(); it != SDmap.end(); it++) {

            // sid is the SDid for the current satellite
         SDid sid = it->first;
         indx = Inext[sid];

            // end of buffer has been reached
         if(Inext[sid] >= SDmap[sid].count.size()) continue;

         if(sid == ref) continue;                     // ignore the reference

            // no data for this satellite at this count
         if(SDmap[sid].count[indx] != count) continue;

            // compute DD phases and DD nominal range
         ddL1 = wl1 * (SDmap[sid].L1[indx] - SDmap[ref].L1[Inext[ref]]);
         ddL2 = wl2 * (SDmap[sid].L2[indx] - SDmap[ref].L2[Inext[ref]]);
         ddP1 = SDmap[sid].P1[indx] - SDmap[ref].P1[Inext[ref]];
         ddP2 = SDmap[sid].P2[indx] - SDmap[ref].P2[Inext[ref]];
         ddER = (SDmap[sid].ER[indx] - SDmap[ref].ER[Inext[ref]]);

            // get the appropriate DDData from the map, or create a new one
         map<DDid,DDData>::iterator jt;
         DDid ddid((ref.ssite == 1 ? ref.site1 : ref.site2),
                   (ref.ssite == 1 ? ref.site2 : ref.site1),sid.sat,ref.sat);
         if(DDDataMap.find(ddid) == DDDataMap.end()) {
               // create a new DDData
            DDData tddb;
            dd = (-ddL1+ddER)/wl1;
            nn1 = int(dd + (dd > 0 ? 0.5 : -0.5));
            tddb.L1bias = wl1 * nn1;
            dd = (-ddL2+ddER)/wl2;
            nn2 = int(dd + (dd > 0 ? 0.5 : -0.5));
            tddb.L2bias = wl2 * nn2;
            oflog << " Phase bias (initial) on " << ddid
               << " at " << setw(3) << count << " "
               << tt.printf("%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g")
               << " L1: " << setw(14) << nn1 << " L2: " << setw(14) << nn2 << endl;
            //tddb.lastresetcount = count;
            tddb.resets.push_back(tddb.count.size());    // always one at beginning
            tddb.prevL1 = (ddL1-ddER)+tddb.L1bias;
            tddb.prevL2 = (ddL2-ddER)+tddb.L2bias;
            DDDataMap[ddid] = tddb;
         }
               
            // get the current DDData structure, and relative sign
         jt = DDDataMap.find(ddid); // never fail...
         ddsign = DDid::compare(ddid,jt->first);
         DDData& ddb=jt->second;
         ok = true;                 // if ok, buffer this DDData = ddb

            // reset bias?
         db1 = (ddsign*(ddL1-ddER) + ddb.L1bias - ddb.prevL1)/wl1;
         db2 = (ddsign*(ddL2-ddER) + ddb.L2bias - ddb.prevL2)/wl2;
         if((CI.Frequency != 2 && fabs(db1) > CI.PhaseBiasReset) ||
            (CI.Frequency != 1 && fabs(db2) > CI.PhaseBiasReset)) {
            long ndb1 = long(db1 + (db1 > 0 ? 0.5 : -0.5));
            long ndb2 = long(db2 + (db2 > 0 ? 0.5 : -0.5));
            oflog << " Phase bias (reset  ) on " << ddid
               << " at " << setw(3) << count << " "
               << tt.printf("%Y/%02m/%02d %2H:%02M:%6.3f=%F/%10.3g")
               << " L1: " << setw(14) << ndb1 << " L2: " << setw(14) << ndb2 << endl;
            ddb.L1bias -= wl1 * ndb1;
            ddb.L2bias -= wl2 * ndb2;
            //ddb.lastresetcount = count;
            ddb.resets.push_back(ddb.count.size());
         }

            // remove the bias from the data
         ddL1 += ddsign * ddb.L1bias;
         ddL2 += ddsign * ddb.L2bias;

            // save for next time
         ddb.prevL1 = ddsign*(ddL1-ddER);
         ddb.prevL2 = ddsign*(ddL2-ddER);

            // buffer the debiased DDs
         ddb.DDL1.push_back(ddsign*ddL1);
         ddb.DDL2.push_back(ddsign*ddL2);
         ddb.DDP1.push_back(ddsign*ddP1);
         ddb.DDP2.push_back(ddsign*ddP2);
         ddb.DDER.push_back(ddsign*ddER);
         ddb.count.push_back(count);

            // increment Inext
         Inext[sid]++;

      }  // end loop over SD's in SDmap

      Inext[ref]++;
      
         // find the next count
         // quit when all Inext are at end
      frst = true;
      ok = false;
      for(it=SDmap.begin(); it != SDmap.end(); it++) {
         if(Inext[it->first] < SDmap[it->first].count.size()) {
            if(frst || it->second.count[Inext[it->first]] < count) {
               count = it->second.count[Inext[it->first]];
               frst = false;
            }
            ok = true;
         }
      }

      if(!ok) break;

   }  // end while loop over epochs

   return 0;
}
catch(Exception& e) { GPSTK_RETHROW(e); }
catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); }
catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------