diff --git a/beidou.hh b/beidou.hh index 930cf53..4b877af 100644 --- a/beidou.hh +++ b/beidou.hh @@ -222,9 +222,33 @@ struct BeidouMessage return alma.pageno; } + // 2^-30 2^-50 int a0gps, a1gps, a0gal, a1gal, a0glo, a1glo, a0utc, a1utc; int8_t deltaTLS; + // in Beidou the offset is a0utc + SOW * a1utc + std::pair getUTCOffset(int tow) const + { + // 2^-30 2^-50 + // a0utc a1utc + double cur = a0utc + ldexp(1.0*tow*a1utc, -20); + double trend = ldexp(a1utc, -20); + + // now in units of 2^-30 seconds, which are ~1.1 nanoseconds each + + double factor = ldexp(1000000000, -30); + return {factor * cur, factor * trend}; + } + + // in Beidou the offset is a0GPS + SOW * a1GPS + std::pair getGPSOffset(int tow) const + { + double cur = a0gps/10.0 + tow*a1gps/10.0; + double trend = a1gps/1.0; + return {cur, trend}; + } + + int parse5(std::basic_string_view cond) { alma.pageno = bbitu(44, 7); diff --git a/ephemeris.hh b/ephemeris.hh index 39bef97..64b1f19 100644 --- a/ephemeris.hh +++ b/ephemeris.hh @@ -164,6 +164,15 @@ struct DopplerData time_t t; }; +template +void getSpeed(int wn, double tow, const T& eph, Vector* v) +{ + Point a, b; + getCoordinates(wn, tow-0.5, eph, &a); + getCoordinates(wn, tow+0.5, eph, &b); + *v = Vector(a, b); +} + template DopplerData doDoppler(int wn, int tow, const Point& us, const T& eph, double freq) { diff --git a/galileo.hh b/galileo.hh index 803a99e..5f78690 100644 --- a/galileo.hh +++ b/galileo.hh @@ -68,8 +68,9 @@ struct GalileoMessage bool disturb1{false}, disturb2{false}, disturb3{false}, disturb4{false}, disturb5{false}; // - // 2^-30 2^-50 1 8-bit week - int32_t a0{0}, a1{0}, t0t{0}, wn0t{0}; + // 2^-30 2^-50 3600 8-bit week + int32_t a0{0}, a1{0}, t0t{0}, wn0t{0}; + // 2^-35 2^-51 3600 8-bit week int32_t a0g{0}, a1g{0}, t0g{0}, wn0g{0}; int8_t dtLS{0}, dtLSF{0}; uint16_t wnLSF{0}; @@ -162,22 +163,68 @@ struct GalileoMessage sisa = getbitu(&page[0], 120, 8); } - int getT0c() + int getT0c() const { return t0c * 60; } + int getT0t() const + { + return t0t * 3600; + } + int getT0g() const + { + return t0g * 3600; + } - std::pair getAtomicOffset(int tow) + std::pair getAtomicOffset(int tow) const { int delta = ephAge(tow, getT0c()); - double cur = af0 + ldexp(delta*af1, -12) + ldexp(delta*delta*af2, -25); - double trend = ldexp(af1, -12) + ldexp(2*delta*af2, -25); + double cur = af0 + ldexp(1.0*delta*af1, -12) + ldexp(1.0*delta*delta*af2, -25); + double trend = ldexp(af1, -12) + ldexp(2.0*delta*af2, -25); // now in units of 2^-34 seconds, which are ~0.058 nanoseconds each double factor = ldexp(1000000000, -34); return {factor * cur, factor * trend}; } + std::pair getUTCOffset(int tow, int wn) const + { + int dw = (int)(uint8_t)wn - (int)(uint8_t) wn0t; + int delta = dw*7*86400 + tow - getT0t(); // NOT ephemeris age tricks + + // 2^-30 2^-50 3600 + // a0 a1 t0t + double cur = a0 + ldexp(1.0*delta*a1, -20); + double trend = ldexp(a1, -20); + + // now in units of 2^-30 seconds, which are ~1.1 nanoseconds each + + double factor = ldexp(1000000000, -30); + // std::cout<<"dw: "< getGPSOffset(int tow, int wn) const + { + int dw = (int)(uint8_t)wn - (int)(uint8_t) wn0g; + int delta = dw*7*86400 + tow - getT0g(); // NOT ephemeris age tricks + + // 2^-35 2^-51 3600 + // a0g a1g t0g + double cur = a0g + ldexp(1.0*delta*a1g, -16); + double trend = ldexp(1.0*a1g, -16); + + // now in units of 2^-35 seconds + + double factor = ldexp(1000000000, -35); // turn into nanoseconds + // std::cout<<"gps dw: "< page) diff --git a/glonass.cc b/glonass.cc index 51aa9a7..183fdfd 100644 --- a/glonass.cc +++ b/glonass.cc @@ -15,6 +15,7 @@ std::basic_string getGlonassMessage(std::basic_string_view pay } +// this does NOT turn it into unix time!! uint32_t GlonassMessage::getGloTime() const { struct tm tm; @@ -33,3 +34,16 @@ uint32_t GlonassMessage::getGloTime() const t += 3600 * (hour) + 60 * minute + seconds; return t - 820368000; // this starts GLONASS time at 31st of december 1995, 00:00 UTC } + +// the 'referencetime' must reflect the time when the frame with Tb was received +uint32_t getGlonassT0e(time_t referencetime, int Tb) +{ + time_t now = referencetime + 3*3600; // this is so we get the Moscow day + struct tm tm; + memset(&tm, 0, sizeof(tm)); + gmtime_r(&now, &tm); + tm.tm_hour = (Tb/4.0); + tm.tm_min = (Tb % 4)*15; + tm.tm_sec = 0; + return timegm(&tm)-3*3600; // and back to UTC +} diff --git a/glonass.hh b/glonass.hh index 32e7aae..2be8525 100644 --- a/glonass.hh +++ b/glonass.hh @@ -82,7 +82,6 @@ struct GlonassMessage y=getbitsglonass(&gstr[0], 85-35, 27); // 2^-11, in kilometers dy=getbitsglonass(&gstr[0], 85-64, 24); // 2^-20, in kilometers ddy=getbitsglonass(&gstr[0], 85-40, 5); // 2^-30, in kilometers - } int32_t z{0}, dz, ddz; @@ -195,3 +194,4 @@ struct GlonassMessage } }; +uint32_t getGlonassT0e(time_t referencetime, int Tb); diff --git a/gps.hh b/gps.hh index 6b30d79..c707c06 100644 --- a/gps.hh +++ b/gps.hh @@ -7,6 +7,74 @@ #include std::basic_string getCondensedGPSMessage(std::basic_string_view payload); +struct GPSAlmanac +{ + int dataid{-1}; + int sv; + uint32_t t0a{0}; + uint32_t e{0}, sqrtA{0}; + int32_t M0, Omega0, deltai, omega, omegadot; + int health; + int af0, af1; + + double getMu() const + { + return 3.986005 * pow(10.0, 14.0); + } // m^3/s^2 + // same for galileo & gps + double getOmegaE() const { return 7.2921151467 * pow(10.0, -5.0);} // rad/s + + + double getE() const + { + return ldexp(e, -21); + } + double getT0e() const + { + return ldexp(t0a, 12); + } + + double getI0() const + { + return M_PI*0.3 + ldexp(M_PI*deltai, -19); + } + + double getOmegadot() const + { + return ldexp(M_PI * omegadot, -38); + } + + double getSqrtA() const + { + return ldexp(sqrtA, -11); + } + + double getOmega0() const + { + return ldexp(M_PI * Omega0, -23); + } + double getOmega() const + { + return ldexp(M_PI * omega, -23); + } + + double getM0() const + { + return ldexp(M_PI * M0, -23); + } + double getIdot() const { return 0; } // radians/s + double getCic() const { return 0; } // radians + double getCis() const { return 0; } // radians + double getCuc() const { return 0; } // radians + double getCus() const { return 0; } // radians + double getCrc() const { return 0; } // meters + double getCrs() const { return 0; } // meters + double getDeltan()const { return 0; } //radians/s + + +}; + + struct GPSState { struct SVIOD @@ -16,6 +84,7 @@ struct GPSState uint32_t t0e; uint32_t e, sqrtA; int32_t m0, omega0, i0, omega, idot, omegadot, deltan; + int16_t cuc{0}, cus{0}, crc{0}, crs{0}, cic{0}, cis{0}; // 16 seconds @@ -26,8 +95,36 @@ struct GPSState // ??? int8_t af2; uint32_t wn{0}, tow{0}; + + double getMu() const + { + return 3.986005 * pow(10.0, 14.0); + } // m^3/s^2 + // same for galileo & gps + double getOmegaE() const { return 7.2921151467 * pow(10.0, -5.0);} // rad/s + + uint32_t getT0e() const { return t0e; } + double getSqrtA() const { return ldexp(sqrtA, -19); } + double getE() const { return ldexp(e, -33); } + double getCuc() const { return ldexp(cuc, -29); } // radians + double getCus() const { return ldexp(cus, -29); } // radians + double getCrc() const { return ldexp(crc, -5); } // meters + double getCrs() const { return ldexp(crs, -5); } // meters + double getM0() const { return ldexp(m0 * M_PI, -31); } // radians + double getDeltan()const { return ldexp(deltan *M_PI, -43); } //radians/s + double getI0() const { return ldexp(i0 * M_PI, -31); } // radians + double getCic() const { return ldexp(cic, -29); } // radians + double getCis() const { return ldexp(cis, -29); } // radians + double getOmegadot() const { return ldexp(omegadot * M_PI, -43); } // radians/s + double getOmega0() const { return ldexp(omega0 * M_PI, -31); } // radians + double getIdot() const { return ldexp(idot * M_PI, -43); } // radians/s + double getOmega() const { return ldexp(omega * M_PI, -31); } // radians + + }; + GPSAlmanac gpsalma; + uint8_t gpshealth{0}; uint16_t ai0{0}; int16_t ai1{0}, ai2{0}; @@ -46,10 +143,26 @@ struct GPSState uint8_t dn; // leap second day number // 1 2^-31 2^-43 2^-55 16 second int ura, af0, af1, af2, t0c; // GPS parameters that should not be here XXX + + int gpsiod{-1}; std::map iods; SVIOD& getEph(int i) { return iods[i]; } // XXXX gps adaptor - void checkCompleteAndClean(int iod){} + void checkCompleteAndClean(int iod) + { + if(iods[iod].words[2] && iods[iod].words[3]) { + if(iods.size() > 1) { + auto tmp = iods[iod]; + iods.clear(); + iods[iod] = tmp; + } + } + } + bool isComplete(int iod) + { + return iods[iod].words[2] && iods[iod].words[3]; + } + }; template @@ -59,7 +172,7 @@ int getT0c(const T& eph) } template -std::pair getAtomicOffset(int tow, const T& eph) +std::pair getGPSAtomicOffset(int tow, const T& eph) { int delta = ephAge(tow, getT0c(eph)); double cur = eph.af0 + ldexp(delta*eph.af1, -12) + ldexp(delta*delta*eph.af2, -24); @@ -70,6 +183,25 @@ std::pair getAtomicOffset(int tow, const T& eph) return {factor * cur, factor * trend}; } +template +std::pair getGPSUTCOffset(int tow, int wn, const T& eph) +{ + // 2^-30 2^-50 3600 + // a0 a1 t0t + + int dw = (int)(uint8_t)wn - (int)(uint8_t) eph.wn0t; + int delta = dw*7*86400 + tow - eph.t0t; // this is pre-scaled for GPS.. + + double cur = eph.a0 + ldexp(1.0*delta*eph.a1, -20); + double trend = ldexp(eph.a1, -20); + + // now in units of 2^-30 seconds, which are ~1.1 nanoseconds each + + double factor = ldexp(1000000000, -30); + return {factor * cur, factor * trend}; +} + + // expects input as 24 bit read to to use messages, returns frame number template @@ -109,8 +241,8 @@ int parseGPSMessage(std::basic_string_view cond, T& out, uint8_t* pagep // out.af0 < cond, T& out, uint8_t* pagep eph.cuc = getbits(&cond[0], 5*24, 16); // 2^-29 RADIANS eph.cus = getbits(&cond[0], 7*24, 16); // 2^-29 RADIANS - out.checkCompleteAndClean(iod); + out.checkCompleteAndClean(out.gpsiod); } else if(frame == 3) { - int iod = getbitu(&cond[0], 9*24, 8); - auto& eph = out.getEph(iod); + out.gpsiod = getbitu(&cond[0], 9*24, 8); + auto& eph = out.getEph(out.gpsiod); eph.words[3]=1; eph.cic = getbits(&cond[0], 2*24, 16); // 2^-29 RADIANS eph.omega0 = getbits(&cond[0], 2*24 + 16, 32); // 2^-31 semi-circles @@ -146,14 +278,14 @@ int parseGPSMessage(std::basic_string_view cond, T& out, uint8_t* pagep eph.omegadot = getbits(&cond[0], 8*24, 24); // 2^-43, semi-circles/s eph.idot = getbits(&cond[0], 9*24+8, 14); // 2^-43, semi-cirlces/s - out.checkCompleteAndClean(iod); + out.checkCompleteAndClean(out.gpsiod); } else if(frame == 4) { // this is a carousel frame int page = getbitu(&cond[0], 2*24 + 2, 6); if(pageptr) *pageptr=0; // cerr<<"Frame 4, page "< cond, T& out, uint8_t* pagep // page 25 -> 63 // 2-10 -> 25 -> 32 ?? } - else if(frame == 5) { // this is a caroussel frame + + if(frame == 5 || frame==4) { // this is a caroussel frame + out.gpsalma.dataid = getbitu(&cond[0], 2*24, 2); + out.gpsalma.sv = getbitu(&cond[0], 2*24+2, 6); + + if(pageptr) + *pageptr= out.gpsalma.sv; + + + out.gpsalma.e = getbitu(&cond[0], 2*24 + 8, 16); + out.gpsalma.t0a = getbitu(&cond[0], 3*24, 8); + out.gpsalma.deltai = getbits(&cond[0], 3*24 +8 , 16); + out.gpsalma.omegadot = getbits(&cond[0], 4*24, 16); + out.gpsalma.health = getbitu(&cond[0], 4*24 +16, 8); + out.gpsalma.sqrtA = getbitu(&cond[0], 5*24, 24); + out.gpsalma.Omega0 = getbits(&cond[0], 6*24, 24); + out.gpsalma.omega = getbits(&cond[0], 7*24, 24); + out.gpsalma.M0 = getbits(&cond[0], 8*24, 24); + out.gpsalma.af0 = (getbits(&cond[0], 9*24, 8) << 3) + getbits(&cond[0], 9*24 +19, 3); + out.gpsalma.af1 = getbits(&cond[0], 9*24 + 8, 11); // cerr<<"Frame 5, SV: "<'; + ret.value = ''; // ret.value=""; ret.value += " "+row.sv; } @@ -74,9 +79,23 @@ function maketable(str, arr) if(row["best-tle-dist"] != null) ret.value = row["best-tle-dist"].toFixed(1) + " km"; } + else if((column == "alma-dist")) { + if(row["alma-dist"] != null) + ret.value = row["alma-dist"].toFixed(1) + " km"; + } + else if(column == "norad") { ret.value = row["best-tle-norad"]; } + else if(column == "delta-utc" && row["delta-utc"] != null) { + ret.value = row["delta-utc"]+'a0: '+row["a0"]+'
a1: '+row["a1"]+'
wn0t: ' + row["wn0t"]+'
t0t: '+row["t0t"]+'
'; + ret.Class = 'CellWithComment'; + } + else if(column == "delta-gps" && row["delta-gps"] != null) { + ret.value = row["delta-gps"]+'a0g: '+row["a0g"]+'
a1g: '+row["a1g"]+'
wn0g: ' +row["wn0g"]+'
t0g: '+row["t0g"]+'
'; + ret.Class = 'CellWithComment'; + } + else if(column == "last-seen-s") { var b = moment.duration(row[column], 's'); ret.value = b.humanize(); @@ -121,7 +140,7 @@ function maketable(str, arr) return ret; }) }) - .enter().append("td").html(function(d) { return d.value; }).attr("align", "right").style("background-color", function(d) { + .enter().append("td").attr("class", function(d) { return d.Class; }).html(function(d) { return d.value; }).attr("align", "right").style("background-color", function(d) { return d.color; }); @@ -193,6 +212,6 @@ function update() }); } -repeat=update(); +update(); diff --git a/html/index.html b/html/index.html index 37f2ea5..aaf7390 100644 --- a/html/index.html +++ b/html/index.html @@ -16,9 +16,33 @@ font-family: monospace; } tr:nth-child(even) {background: #CCC} tr:nth-child(odd) {background: #FFF} + +.CellWithComment{ + position:relative; +} + +.CellComment{ + display:none; + position:absolute; + z-index:100; + border:1px; + background-color:white; + border-style:solid; + border-width:1px; + border-color:red; + padding:3px; + color:red; + top:20px; + left:20px; +} + +.CellWithComment:hover span.CellComment{ + display:block; +} + - Last update:
+ Last update: . More information about this Galileo/GPS/BeiDou/Glonass open source monitor can be found here.

diff --git a/navdump.cc b/navdump.cc index b44d695..bed8234 100644 --- a/navdump.cc +++ b/navdump.cc @@ -26,17 +26,6 @@ #include using namespace std; -static std::string humanTime(time_t t) -{ - struct tm tm={0}; - gmtime_r(&t, &tm); - - char buffer[80]; - strftime(buffer, sizeof(buffer), "%a, %d %b %Y %T %z", &tm); - return buffer; -} - - string beidouHealth(int in) { string ret; @@ -61,6 +50,12 @@ string beidouHealth(int in) return ret; } +double utcFromGPS(int wn, double tow) +{ + return (315964800 + wn * 7*86400 + tow - 18); +} + + int main(int argc, char** argv) try { @@ -71,9 +66,17 @@ try tles.parseFile("gps-ops.txt"); tles.parseFile("beidou.txt"); + bool skipGPS{false}; + bool skipBeidou{false}; + bool skipGalileo{false}; + bool skipGlonass{false}; + ofstream almanac("almanac.txt"); ofstream iodstream("iodstream.csv"); iodstream << "timestamp gnssid sv iodnav t0e age" << endl; + + ofstream csv("delta.csv"); + csv <<"timestamp gnssid sv tow tle_distance alma_distance utc_dist x y z vx vy vz rad inclination e iod"< inav((uint8_t*)nmm.gi().contents().c_str(), nmm.gi().contents().size()); static map gms; static map, GalileoMessage> gmwtypes; @@ -181,29 +186,93 @@ try cout<((uint8_t*)nmm.gpsi().contents().c_str(), nmm.gpsi().contents().size())); struct GPSState gs; + static map eph; + static map almas; uint8_t page; + static int gpswn; int frame=parseGPSMessage(cond, gs, &page); - cout<<"GPS "< oldgs1s; - + gpswn = gs.wn; cout << "gpshealth = "<<(int)gs.gpshealth<<", wn "<second.t0c != gs.t0c) { - auto oldOffset = getAtomicOffset(gs.tow, iter->second); - auto newOffset = getAtomicOffset(gs.tow, gs); + auto oldOffset = getGPSAtomicOffset(gs.tow, iter->second); + auto newOffset = getGPSAtomicOffset(gs.tow, gs); cout<<" Timejump: "<second) )<<" seconds, old t0c "<second.t0c; } oldgs1s[sv] = gs; } else if(frame == 2) { - cout << "t0e = "<second.t0e << " " <second.t0e); + parseGPSMessage(cond, eph[sv], &page); + cout << "t0e = "<second.t0e << " " <second.t0e) << " iod "<= 25 && gs.gpsalma.sv <= 32) || gs.gpsalma.sv==57 ) { // see table 20-V of the GPS ICD + cout << " data-id "<((uint8_t*)nmm.bid1().contents().c_str(), nmm.bid1().contents().size())); uint8_t pageno; @@ -282,6 +351,9 @@ try } else if(nmm.type() == NavMonMessage::GlonassInavType) { + if(skipGlonass) + continue; + static map gms; auto& gm = gms[nmm.gloi().gnsssv()]; @@ -292,10 +364,12 @@ try cout << ", hour "<<(int)gm.hour <<" minute " <<(int)gm.minute <<" seconds "<<(int)gm.seconds; // start of period is 1st of January 1996 + (n4-1)*4, 03:00 UTC time_t glotime = gm.getGloTime(); - cout<<" 'wn' " << glotime / (7*86400)<<" 'tow' "<< (glotime % (7*86400)); + cout<<" 'wn' " << glotime / (7*86400)<<" 'tow' "<< (glotime % (7*86400)) << " x " < "; cout << "("< #include +#include +#include + struct EofException{}; size_t readn2(int fd, void* buffer, size_t len); +std::string humanTime(time_t t); diff --git a/navparse.cc b/navparse.cc index ee2a5ca..a5f28de 100644 --- a/navparse.cc +++ b/navparse.cc @@ -10,6 +10,7 @@ #include #include #include +#include #include "ext/powerblog/h2o-pp.hh" #include "minicurl.hh" #include @@ -30,9 +31,42 @@ using namespace std; std::map g_srcpos; +template +class GetterSetter +{ +public: + void set(const T& t) + { + std::lock_guard mut(d_lock); + d_t = t; + } + + T get() + { + T ret; + { + std::lock_guard mut(d_lock); + ret = d_t; + } + return ret; + } +private: + T d_t; + std::mutex d_lock; +}; + + + map g_beidoualma; map> g_glonassalma; map g_galileoalma; +map g_gpsalma; + +GetterSetter> g_beidoualmakeeper; +GetterSetter>> g_glonassalmakeeper; +GetterSetter> g_galileoalmakeeper; +GetterSetter> g_gpsalmakeeper; + TLERepo g_tles; struct GNSSReceiver { @@ -222,8 +256,10 @@ struct SVStat uint16_t wnLSF{0}; uint8_t dn; // leap second day number // 1 2^-31 2^-43 2^-55 16 second - int ura, af0, af1, af2, t0c; // GPS parameters that should not be here XXX + int ura, af0, af1, af2, t0c, gpsiod; // GPS parameters that should not be here XXX (or perhaps they should) + GPSAlmanac gpsalma; + // beidou: int t0eMSB{-1}, t0eLSB{-1}, aode{-1}, aodc{-1}; BeidouMessage beidouMessage, oldBeidouMessage; @@ -326,11 +362,13 @@ void SVStat::addGalileoWord(std::basic_string_view page) sf5 = getbitu(&page[0], 46, 1); BGDE1E5a = getbits(&page[0], 47, 10); BGDE1E5b = getbits(&page[0], 57, 10); - + e5bhs = getbitu(&page[0], 67, 2); e1bhs = getbitu(&page[0], 69, 2); e5bdvs = getbitu(&page[0], 71, 1); e1bdvs = getbitu(&page[0], 72, 1); + + wn = getbitu(&page[0], 73, 12); if(tow != getbitu(&page[0], 85, 20)) { @@ -361,41 +399,35 @@ void SVStat::addGalileoWord(std::basic_string_view page) } } -template -void getSpeed(int wn, double tow, const T& eph, Vector* v) -{ - Point a, b; - getCoordinates(wn, tow-0.5, eph, &a); - getCoordinates(wn, tow+0.5, eph, &b); - *v = Vector(a, b); -} +typedef std::map, SVStat> svstats_t; +svstats_t g_svstats; -std::map, SVStat> g_svstats; +GetterSetter, SVStat>> g_statskeeper; -int latestWN(int gnssid) +int latestWN(int gnssid, const svstats_t& stats) { map> ages; - for(const auto& s: g_svstats) + for(const auto& s: stats) if(s.first.first == gnssid) ages[7*s.second.wn*86400 + s.second.tow]= s.first; if(ages.empty()) throw runtime_error("Asked for latest WN for "+to_string(gnssid)+": we don't know it yet"); - return g_svstats[ages.rbegin()->second].wn; + + return stats.find(ages.rbegin()->second)->second.wn; } -int latestTow(int gnssid) +int latestTow(int gnssid, const svstats_t& stats) { map> ages; - for(const auto& s: g_svstats) + for(const auto& s: stats) if(s.first.first == gnssid) ages[7*s.second.wn*86400 + s.second.tow]= s.first; if(ages.empty()) throw runtime_error("Asked for latest TOW for "+to_string(gnssid)+": we don't know it yet"); - return g_svstats[ages.rbegin()->second].tow; + return stats.find(ages.rbegin()->second)->second.tow; } - uint64_t nanoTime(int gnssid, int wn, int tow) { int offset; @@ -493,9 +525,10 @@ struct InfluxPusher /* The GST start epoch is defined as 13 seconds before midnight between 21st August and 22nd August 1999, i.e. GST was equal to 13 seconds at 22nd August 1999 00:00:00 UTC. */ -std::string humanTime(int wn, int tow) +std::string humanTime(int gnssid, int wn, int tow) { - time_t t = utcFromGST(wn, tow); + time_t t = nanoTime(gnssid, wn, tow)/1000000000; + struct tm tm; gmtime_r(&t, &tm); @@ -504,12 +537,12 @@ std::string humanTime(int wn, int tow) return buffer; } -std::optional getHzCorrection(time_t now) +std::optional getHzCorrection(time_t now, const svstats_t svstats) { int galcount{0}, gpscount{0}, allcount{0}; double galtot{0}, gpstot{0}, alltot{0}; - for(const auto& s: g_svstats) { + for(const auto& s: svstats) { if(now - s.second.deltaHz.first < 60) { alltot+=s.second.deltaHz.second; allcount++; @@ -565,9 +598,11 @@ double getElevationDeg(const Point& p, int sourceid) return 90.0 - deg; } - - - +std::string humanBhs(int bhs) +{ + static vector options{"ok", "out of service", "will be out of service", "test"}; + return options.at(bhs); +} int main(int argc, char** argv) try { @@ -577,6 +612,7 @@ try g_tles.parseFile("glo-ops.txt"); g_tles.parseFile("gps-ops.txt"); g_tles.parseFile("beidou.txt"); + signal(SIGPIPE, SIG_IGN); InfluxPusher idb(argc > 3 ? argv[3] : "galileo"); @@ -587,15 +623,17 @@ try h2s.addHandler("/global", [](auto handler, auto req) { nlohmann::json ret = nlohmann::json::object(); + auto svstats = g_statskeeper.get(); ret["leap-seconds"] = g_dtLS; try { - ret["last-seen"]=utcFromGST(latestWN(2), latestTow(2)); + ret["last-seen"]=utcFromGST(latestWN(2, svstats), latestTow(2, svstats)); } catch(...) {} map utcstats, gpsgststats, gpsutcstats; - for(const auto& s: g_svstats) { + + for(const auto& s: svstats) { if(!s.second.wn) // this will suck in 20 years continue; @@ -604,8 +642,8 @@ try if(s.first.first == 0) { // GPS int sv = s.first.second; - int dw = (uint8_t)g_svstats[{0,sv}].wn - g_svstats[{0,sv}].wn0t; - int age = dw * 7 * 86400 + g_svstats[{0,sv}].tow - g_svstats[{0,sv}].t0t; // t0t is PRESCALED + int dw = (uint8_t)svstats[{0,sv}].wn - svstats[{0,sv}].wn0t; + int age = dw * 7 * 86400 + svstats[{0,sv}].tow - svstats[{0,sv}].t0t; // t0t is PRESCALED gpsutcstats[age]=s.first.second; continue; @@ -626,9 +664,9 @@ try } else { int sv = utcstats.begin()->second; // freshest SV - long shift = g_svstats[{2,sv}].a0 * (1LL<<20) + g_svstats[{2,sv}].a1 * utcstats.begin()->first; // in 2^-50 seconds units + long shift = svstats[{2,sv}].a0 * (1LL<<20) + svstats[{2,sv}].a1 * utcstats.begin()->first; // in 2^-50 seconds units ret["utc-offset-ns"] = 1.073741824*ldexp(1.0*shift, -20); - ret["leap-second-planned"] = (g_svstats[{2,sv}].dtLSF != g_svstats[{2,sv}].dtLS); + ret["leap-second-planned"] = (svstats[{2,sv}].dtLSF != svstats[{2,sv}].dtLS); } if(gpsgststats.empty()) { @@ -636,7 +674,7 @@ try } else { int sv = gpsgststats.begin()->second; // freshest SV - long shift = g_svstats[{2,sv}].a0g * (1L<<16) + g_svstats[{2,sv}].a1g * gpsgststats.begin()->first; // in 2^-51 seconds units + long shift = svstats[{2,sv}].a0g * (1L<<16) + svstats[{2,sv}].a1g * gpsgststats.begin()->first; // in 2^-51 seconds units ret["gps-offset-ns"] = 1.073741824*ldexp(shift, -21); } @@ -646,7 +684,7 @@ try } else { int sv = gpsutcstats.begin()->second; // freshest SV - long shift = g_svstats[{0,sv}].a0 * (1LL<<20) + g_svstats[{0,sv}].a1 * gpsutcstats.begin()->first; // In 2^-50 seconds units + long shift = svstats[{0,sv}].a0 * (1LL<<20) + svstats[{0,sv}].a1 * gpsutcstats.begin()->first; // In 2^-50 seconds units ret["gps-utc-offset-ns"] = 1.073741824*ldexp(shift, -20); } @@ -657,13 +695,17 @@ try }); h2s.addHandler("/almanac", [](auto handler, auto req) { + auto beidoualma = g_beidoualmakeeper.get(); + auto svstats = g_statskeeper.get(); nlohmann::json ret = nlohmann::json::object(); - for(const auto& ae : g_beidoualma) { + for(const auto& ae : beidoualma) { + nlohmann::json item = nlohmann::json::object(); + item["gnssid"]=3; if(ae.second.alma.getT0e() > 7*86400) continue; Point sat; - getCoordinates(0, latestTow(3), ae.second.alma, &sat); + getCoordinates(0, latestTow(3, svstats), ae.second.alma, &sat); item["eph-ecefX"]= sat.x/1000; item["eph-ecefY"]= sat.y/1000; item["eph-ecefZ"]= sat.z/1000; @@ -672,9 +714,10 @@ try item["eph-longitude"] = 180*longlat.first/M_PI; item["eph-latitude"]= 180*longlat.second/M_PI; item["t0e"] = ae.second.alma.getT0e(); - item["t"]= ephAge(ae.second.alma.getT0e(), latestTow(3))/86400.0; - if(ephAge(ae.second.alma.getT0e(), latestTow(3)) < 0) { - auto match = g_tles.getBestMatch(nanoTime(3, latestWN(3), latestTow(3))/1000000000.0, + item["t"]= ephAge(ae.second.alma.getT0e(), latestTow(3, svstats))/86400.0; + item["inclination"] = 180 * ae.second.alma.getI0() /M_PI; + if(ephAge(ae.second.alma.getT0e(), latestTow(3, svstats)) < 0) { + auto match = g_tles.getBestMatch(nanoTime(3, latestWN(3, svstats), latestTow(3, svstats))/1000000000.0, sat.x, sat.y, sat.z); if(match.distance < 200000) { @@ -700,11 +743,12 @@ try ret[fmt::sprintf("C%02d", ae.first)] = item; } - for(const auto& ae : g_glonassalma) { + auto glonassalma = g_glonassalmakeeper.get(); + for(const auto& ae : glonassalma) { nlohmann::json item = nlohmann::json::object(); // ae.second.first -> even ae.second.sceond -> odd - + item["gnssid"]=6; item["e"] = ae.second.first.getE(); item["inclination"] = 180 * ae.second.first.getI0() /M_PI; item["health"] = ae.second.first.CnA; @@ -715,18 +759,20 @@ try ret[fmt::sprintf("R%02d", ae.first)] = item; } - for(const auto& ae : g_galileoalma) { + auto galileoalma = g_galileoalmakeeper.get(); + for(const auto& ae : galileoalma) { nlohmann::json item = nlohmann::json::object(); + item["gnssid"]=2; item["e"] = ae.second.getE(); item["e1bhs"] = ae.second.e1bhs; item["e5bhs"] = ae.second.e5bhs; item["t0e"] = ae.second.getT0e(); - item["t"]= ephAge(ae.second.getT0e(), latestTow(2))/86400.0; - item["eph-age"] = ephAge(latestTow(2), ae.second.getT0e()); + item["t"]= ephAge(ae.second.getT0e(), latestTow(2, svstats))/86400.0; + item["eph-age"] = ephAge(latestTow(2, svstats), ae.second.getT0e()); item["i0"] = 180.0 * ae.second.getI0()/ M_PI; - + item["inclination"] = 180 * ae.second.getI0() /M_PI; Point sat; - getCoordinates(0, latestTow(2), ae.second, &sat); + getCoordinates(0, latestTow(2, svstats), ae.second, &sat); item["eph-ecefX"]= sat.x/1000; item["eph-ecefY"]= sat.y/1000; item["eph-ecefZ"]= sat.z/1000; @@ -735,7 +781,7 @@ try item["eph-longitude"] = 180*longlat.first/M_PI; item["eph-latitude"]= 180*longlat.second/M_PI; - auto match = g_tles.getBestMatch(nanoTime(2, latestWN(2), latestTow(2))/1000000000.0, + auto match = g_tles.getBestMatch(nanoTime(2, latestWN(2, svstats), latestTow(2, svstats))/1000000000.0, sat.x, sat.y, sat.z); if(match.distance < 200000) { @@ -756,20 +802,81 @@ try item["tle-longitude"] = 180*match.longitude/M_PI; item["tle-altitude"] = match.altitude; } - - ret[fmt::sprintf("E%02d", ae.first)] = item; } + + auto gpsalma = g_gpsalmakeeper.get(); + for(const auto& ae : gpsalma) { + nlohmann::json item = nlohmann::json::object(); + item["gnssid"]=0; + item["e"] = ae.second.getE(); + item["health"] = ae.second.health; + item["t0e"] = ae.second.getT0e(); + item["t"]= ephAge(ae.second.getT0e(), latestTow(2, svstats))/86400.0; + item["eph-age"] = ephAge(latestTow(2, svstats), ae.second.getT0e()); + item["i0"] = 180.0 * ae.second.getI0()/ M_PI; + item["inclination"] = 180 * ae.second.getI0() /M_PI; + Point sat; + getCoordinates(0, latestTow(0, svstats), ae.second, &sat); + item["eph-ecefX"]= sat.x/1000; + item["eph-ecefY"]= sat.y/1000; + item["eph-ecefZ"]= sat.z/1000; + + auto longlat = getLongLat(sat.x, sat.y, sat.z); + item["eph-longitude"] = 180*longlat.first/M_PI; + item["eph-latitude"]= 180*longlat.second/M_PI; + + auto match = g_tles.getBestMatch(nanoTime(0, latestWN(0, svstats), latestTow(0, svstats))/1000000000.0, + sat.x, sat.y, sat.z); + + if(match.distance < 200000) { + item["best-tle"] = match.name; + item["best-tle-norad"] = match.norad; + item["best-tle-int-desig"] = match.internat; + item["best-tle-dist"] = match.distance/1000.0; + + item["tle-ecefX"] = match.ecefX/1000; + item["tle-ecefY"] = match.ecefY/1000; + item["tle-ecefZ"] = match.ecefZ/1000; + + item["tle-eciX"] = match.eciX/1000; + item["tle-eciY"] = match.eciY/1000; + item["tle-eciZ"] = match.eciZ/1000; + + item["tle-latitude"] = 180*match.latitude/M_PI; + item["tle-longitude"] = 180*match.longitude/M_PI; + item["tle-altitude"] = match.altitude; + } + ret[fmt::sprintf("G%02d", ae.first)] = item; + } return ret; }); + + h2s.addHandler("/observers", [](auto handler, auto req) { + nlohmann::json ret = nlohmann::json::array(); + for(const auto& src : g_srcpos) { + nlohmann::json obj; + obj["id"] = src.first; + auto longlat = getLongLat(src.second.x, src.second.y, src.second.z); + longlat.first *= 180.0/M_PI; + longlat.second *= 180.0/M_PI; + longlat.first = ((int)(10*longlat.first))/10.0; + longlat.second = ((int)(10*longlat.second))/10.0; + obj["longitude"] = longlat.first; + obj["latitude"] = longlat.second; + ret.push_back(obj); + } + return ret; + }); h2s.addHandler("/svs", [](auto handler, auto req) { + auto svstats = g_statskeeper.get(); nlohmann::json ret = nlohmann::json::object(); - auto hzCorrection = getHzCorrection(time(0)); + auto hzCorrection = getHzCorrection(time(0), svstats); - for(const auto& s: g_svstats) { + for(const auto& s: svstats) { nlohmann::json item = nlohmann::json::object(); if(!s.second.tow) // I know, I know, will suck briefly continue; @@ -778,7 +885,6 @@ try item["svid"] = s.first.second; // perhaps check oldBeidouMessage for sow >=0 as 'completeIOD'? - if(s.first.first == 3) { // beidou item["sisa"]=humanUra(s.second.ura); if(s.second.t0eMSB >= 0 && s.second.t0eLSB >=0) @@ -794,7 +900,14 @@ try item["best-tle-dist"] = s.second.tleMatch.distance /1000.0; item["best-tle-norad"] = s.second.tleMatch.norad; item["best-tle-int-desig"] = s.second.tleMatch.internat; - + } + Point p; + getCoordinates(0, s.second.tow, s.second.oldBeidouMessage, &p); + auto beidoualma = g_beidoualmakeeper.get(); + if(auto iter = beidoualma.find(s.first.second); iter != beidoualma.end()) { + Point almapos; + getCoordinates(0, s.second.tow, iter->second.alma, &almapos); + item["alma-dist"] = Vector(almapos, p).length()/1000.0; } } else if(s.first.first == 6) { // glonass @@ -802,6 +915,14 @@ try item["sisa"] = humanFt(s.second.glonassMessage.FT); item["aode"] = s.second.aode; item["iod"] = s.second.glonassMessage.Tb; + + time_t glonow = nanoTime(6, s.second.wn, s.second.tow)/1000000000.0; + // the 820368000 stuff is to rebase to 'day 1' so the % works + auto pseudoTow = (getGlonassT0e(glonow, s.second.glonassMessage.Tb) - 820368000) % (7*86400); + + // cout<=0) { item["best-tle"] = s.second.tleMatch.name; item["best-tle-dist"] = s.second.tleMatch.distance /1000.0; @@ -850,10 +971,57 @@ try if(hzCorrection) item["delta_hz_corr"] = s.second.deltaHz.second - *hzCorrection; } + if(s.first.first == 0) { + auto gpsalma = g_gpsalmakeeper.get(); + if(auto iter = gpsalma.find(s.first.second); iter != gpsalma.end()) { + Point almapos; + getCoordinates(0, s.second.tow, iter->second, &almapos); + item["alma-dist"] = Vector(almapos, p).length()/1000.0; + } + } + if(s.first.first == 2) { + auto galileoalma = g_galileoalmakeeper.get(); + if(auto iter = galileoalma.find(s.first.second); iter != galileoalma.end()) { + Point almapos; + getCoordinates(0, s.second.tow, iter->second, &almapos); + item["alma-dist"] = Vector(almapos, p).length()/1000.0; + } + } + + } item["a0"]=s.second.a0; item["a1"]=s.second.a1; + if(s.first.first == 0) { // GPS + auto deltaUTC = getGPSUTCOffset(s.second.tow, s.second.wn, s.second); + item["delta-utc"] = fmt::sprintf("%.1f %+.1f/d", deltaUTC.first, deltaUTC.second * 86400); + item["t0t"] = s.second.t0t; + item["wn0t"] = s.second.wn0t; + } + if(s.first.first == 2) { + auto deltaUTC = s.second.galmsg.getUTCOffset(s.second.tow, s.second.wn); + item["delta-utc"] = fmt::sprintf("%.1f %+.1f/d", deltaUTC.first, deltaUTC.second * 86400); + auto deltaGPS = s.second.galmsg.getGPSOffset(s.second.tow, s.second.wn); + item["delta-gps"] = fmt::sprintf("%.1f %+.1f/d", deltaGPS.first, deltaGPS.second * 86400); + item["t0t"] = s.second.galmsg.t0t; + item["wn0t"] = s.second.galmsg.wn0t; + } + if(s.first.first == 3) { + auto deltaUTC = s.second.oldBeidouMessage.getUTCOffset(s.second.oldBeidouMessage.sow); + item["delta-utc"] = fmt::sprintf("%.1f %+.1f/d", deltaUTC.first, deltaUTC.second * 86400); + + auto deltaGPS = s.second.oldBeidouMessage.getGPSOffset(s.second.oldBeidouMessage.sow); + item["delta-gps"] = fmt::sprintf("%.1f %+.1f/d", deltaGPS.first, deltaGPS.second * 86400); + item["t0g"] =0; + item["wn0g"] = 0; + + item["t0t"] = 0; + item["wn0t"] = 0; + } + + + item["dtLS"]=s.second.dtLS; if(s.first.first == 3) { // beidou @@ -870,10 +1038,12 @@ try if(s.first.first == 2) { // galileo item["a0g"]=s.second.a0g; item["a1g"]=s.second.a1g; - vector options{"ok", "out of service", "will be out of service", "test"}; + item["t0g"]=s.second.t0g; + item["wn0g"]=s.second.wn0g; + item["health"] = - options[s.second.e1bhs] +"/" + - options[s.second.e5bhs] +"/" + + humanBhs(s.second.e1bhs) +"/" + + humanBhs(s.second.e5bhs) +"/" + (s.second.e1bdvs ? "NG" : "val") +"/"+ (s.second.e5bdvs ? "NG" : "val"); item["e5bdvs"]=s.second.e5bdvs; @@ -902,7 +1072,7 @@ try Point sat; if((s.first.first == 0 || s.first.first == 2) && s.second.completeIOD()) - getCoordinates(0, latestTow(s.first.first), s.second.liveIOD(), & sat); + getCoordinates(0, latestTow(s.first.first, svstats), s.second.liveIOD(), & sat); if(sat.x) { det["elev"] = getElevationDeg(sat, pr.first); @@ -937,18 +1107,24 @@ try const char *address = argc > 1 ? argv[1] : "127.0.0.1:29599"; std::thread ws([&h2s, address]() { auto actx = h2s.addContext(); - h2s.addListener(ComboAddress(address), actx); - cout<<"Listening on "<< address < ["<< humanBhs(gm.e5bhs)<<", "<< humanBhs(gm.e1bhs)<<", "<< (int)gm.e5bdvs <<", " << (int)gm.e1bdvs<<"], lastseen "<=1 && wtype <= 4) { // ephemeris uint16_t iod = getbitu(&inav[0], 6, 10); @@ -1109,7 +1294,7 @@ try int ephage = ephAge(ent.second.tow, ent.second.prevIOD.second.t0e); if(ent.second.liveIOD().sisa != ent.second.prevIOD.second.sisa) { - cout< id{nmm.rfd().gnssid(), nmm.rfd().gnsssv()}; if(id.first == 3 && g_svstats[id].oldBeidouMessage.sow >= 0 && g_svstats[id].oldBeidouMessage.sqrtA != 0) { auto res = doDoppler(nmm.rfd().rcvwn(), nmm.rfd().rcvtow(), g_srcpos[nmm.sourceid()], g_svstats[id].oldBeidouMessage, 1561.098 * 1000000); @@ -1172,18 +1356,11 @@ try } time_t t = utcFromGPS(nmm.rfd().rcvwn(), nmm.rfd().rcvtow()); - if(0) { - dopplercsv << std::fixed << t <<" " << nmm.rfd().gnssid() <<" " < 10) { g_svstats[id].deltaHz = {t, nmm.rfd().doppler() - res.preddop}; idb.addValue(id, "delta_hz", nmm.rfd().doppler() - res.preddop); - auto corr = getHzCorrection(t); + auto corr = getHzCorrection(t, g_svstats); if(corr) { idb.addValue(id, "delta_hz_cor", nmm.rfd().doppler() - res.preddop - (1561.098/1575.42) * (*corr)); } @@ -1192,20 +1369,15 @@ try else if(g_svstats[id].completeIOD()) { auto res = doDoppler(nmm.rfd().rcvwn(), nmm.rfd().rcvtow(), g_srcpos[nmm.sourceid()], g_svstats[id].liveIOD(), 1575.42 * 1000000); time_t t = utcFromGPS(nmm.rfd().rcvwn(), nmm.rfd().rcvtow()); - if(0) - dopplercsv << std::fixed << t <<" " << nmm.rfd().gnssid() <<" " < 10) { g_svstats[id].deltaHz = {t, nmm.rfd().doppler() - res.preddop}; idb.addValue(id, "delta_hz", nmm.rfd().doppler() - res.preddop); - auto corr = getHzCorrection(t); + auto corr = getHzCorrection(t, g_svstats); if(corr) { idb.addValue(id, "delta_hz_cor", nmm.rfd().doppler() - res.preddop - *corr); } } - - // cout<<"Had doppler for "< 1) { // XX find better way to check + cout< [" << humanUra(svstat.ura)<<"] "<<(int)svstat.ura<<", lastseen "< [" << (int)svstat.gpshealth <<"], lastseen "< ["<< (int)svstat.gpshealth<<"] , lastseen "<=25 && page<=32)) { + g_gpsalma[svstat.gpsalma.sv] = svstat.gpsalma; + } g_svstats[id].perrecv[nmm.sourceid()].t = nmm.localutcseconds(); g_svstats[id].tow = nmm.gpsi().gnsstow(); g_svstats[id].wn = nmm.gpsi().gnsswn(); @@ -1266,7 +1452,7 @@ try uint8_t pageno; auto cond = getCondensedBeidouMessage(std::basic_string((uint8_t*)nmm.bid1().contents().c_str(), nmm.bid1().contents().size())); auto& bm = svstat.beidouMessage; - + auto oldbm = bm; int fraid=bm.parse(cond, &pageno); svstat.tow = nmm.bid1().gnsstow(); svstat.wn = nmm.bid1().gnsswn(); @@ -1278,11 +1464,17 @@ try svstat.af2 = bm.a2; svstat.aode = bm.aode; svstat.aodc = bm.aodc; + + if(oldbm.sath1 != bm.sath1) { + cout<=0 && svstat.lastBeidouMessage1.t0c != bm.t0c) { auto oldOffset = svstat.lastBeidouMessage1.getAtomicOffset(bm.sow); auto newOffset = bm.getAtomicOffset(bm.sow); @@ -1311,9 +1503,10 @@ try if(svstat.oldBeidouMessage.getT0e() != svstat.beidouMessage.getT0e()) { getCoordinates(svstat.wn, svstat.tow, svstat.oldBeidouMessage, &oldpoint); Vector jump(oldpoint ,newpoint); - cout< (%f, %f, %f), jump: %f, seconds: %f\n", + /* cout< (%f, %f, %f), jump: %f, seconds: %f\n", id.second, oldpoint.x, oldpoint.y, oldpoint.z, newpoint.x, newpoint.y, newpoint.z, jump.length(), (double)bm.getT0e() - svstat.oldBeidouMessage.getT0e()); + */ double hours = (bm.getT0e() - svstat.oldBeidouMessage.getT0e())/3600; if(hours < 4) { svstat.latestDisco = jump.length(); @@ -1350,11 +1543,6 @@ try g_dtLSBeidou = bm.deltaTLS; // cout<<"Beidou leap seconds: "<((uint8_t*)nmm.bid2().contents().c_str(), nmm.bid2().contents().size())); @@ -1381,9 +1569,14 @@ try if(strno == 1 && gm.n4 != 0 && gm.NT !=0) { uint32_t glotime = gm.getGloTime(); // this starts GLONASS time at 31st of december 1995, 00:00 UTC svstat.wn = glotime / (7*86400); - svstat.tow = glotime % (7*86400); + svstat.tow = glotime % (7*86400); + // cout<<"Glonass now: "<