rough code for coverage maps

experimental
bert hubert 2019-10-12 14:57:44 +02:00
parent 34ff892bc2
commit 80a062a69f
5 changed files with 283 additions and 157 deletions

View File

@ -22,7 +22,7 @@ clean:
rm -f ext/fmt-5.2.1/src/format.o
navparse: navparse.o ext/fmt-5.2.1/src/format.o $(H2OPP) $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o
navparse: navparse.o ext/fmt-5.2.1/src/format.o $(H2OPP) $(SIMPLESOCKETS) minicurl.o ubx.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o navmon.o coverage.o
$(CXX) -std=gnu++17 $^ -o $@ -pthread -L/usr/local/lib -L/usr/local/opt/openssl/lib/ -lh2o-evloop -lssl -lcrypto -lz -lcurl -lprotobuf $(WSLAY)
navdump: navdump.o ext/fmt-5.2.1/src/format.o bits.o navmon.pb.o gps.o ephemeris.o beidou.o glonass.o navmon.o $(patsubst %.cc,%.o,$(wildcard ext/sgp4/libsgp4/*.cc)) tle.o sp3.o

82
coverage.cc 100644
View File

@ -0,0 +1,82 @@
#include <map>
#include "galileo.hh"
#include "minivec.hh"
#include "navmon.hh"
#include "navparse.hh"
#include "ephemeris.hh"
#include "fmt/format.h"
#include "fmt/printf.h"
#include <fstream>
using namespace std;
extern GetterSetter<map<int, GalileoMessage::Almanac>> g_galileoalmakeeper;
extern GetterSetter<svstats_t> g_statskeeper;
// a vector of pairs of latidude,vector<longitude,numsats>
typedef vector<pair<double,vector<pair<double, int> > > > covmap_t;
covmap_t emitCoverage()
{
covmap_t ret;
ofstream cmap("covmap.csv");
cmap<<"latitude longitude count5 count10 count20"<<endl;
map<int, Point> sats;
auto galileoalma = g_galileoalmakeeper.get();
auto svstats = g_statskeeper.get();
auto pseudoTow = (time(0) - 820368000) % (7*86400);
cout<<"pseudoTow "<<pseudoTow<<endl;
for(const auto &g : galileoalma) {
Point sat;
getCoordinates(pseudoTow, g.second, &sat);
if(g.first < 0)
continue;
if(svstats[{2,g.first,1}].completeIOD() && svstats[{2,g.first,1}].liveIOD().sisa == 255) {
cout<<g.first<<" NAPA!"<<endl;
continue;
}
sats[g.first]=sat;
}
cout<<"Have "<<sats.size()<<" SVs active"<<endl;
double R = 6371000;
for(int latitude = 90 ; latitude > -90; latitude-=1) { // north-south
double phi = M_PI* latitude / 180;
double longsteps = 1 + 360.0 * cos(phi);
double step = 360.0 / longsteps;
cout<<"lat "<<latitude<<" Step size: "<<step<<endl;
vector<pair<double, int>> latvect;
for(double longitude = -180; longitude < 180; longitude += step) { // east - west
cout<<" long "<<longitude<<endl;
Point p;
// phi = latitude, lambda = longitude
double lambda = M_PI* longitude / 180;
p.x = R * cos(phi) * cos(lambda);
p.y = R * cos(phi) * sin(lambda);
p.z = R * sin(phi);
if(longitude == -180) {
auto longlat = getLongLat(p.x, p.y, p.z);
// cout<<fmt::sprintf("%3.0f ", 180.0*longlat.second/M_PI);
}
int numsats5=0, numsats10=0, numsats20=0;
for(const auto& s : sats) {
// double getElevationDeg(const Point& sat, const Point& our);
double elev = getElevationDeg(s.second, p);
if(elev > 5.0)
numsats5++;
if(elev > 10.0)
numsats10++;
if(elev > 20.0)
numsats20++;
}
latvect.push_back(make_pair(longitude, numsats10));
cmap << longitude <<" " <<latitude <<" " << numsats5 << " " <<numsats10<<" "<<numsats20<<endl;
}
ret.push_back(make_pair(latitude, latvect));
cout<<endl;
}
return ret;
}

View File

@ -4,6 +4,7 @@
#include <time.h>
#include <string>
#include <tuple>
#include <mutex>
struct EofException{};
size_t readn2(int fd, void* buffer, size_t len);
@ -19,3 +20,28 @@ struct SatID
return std::tie(gnss, sv, sigid) < std::tie(rhs.gnss, rhs.sv, rhs.sigid);
}
};
template<typename T>
class GetterSetter
{
public:
void set(const T& t)
{
std::lock_guard<std::mutex> mut(d_lock);
d_t = t;
}
T get()
{
T ret;
{
std::lock_guard<std::mutex> mut(d_lock);
ret = d_t;
}
return ret;
}
private:
T d_t;
std::mutex d_lock;
};

View File

@ -27,6 +27,7 @@
#include <optional>
#include "navmon.hh"
#include <Tle.h>
#include "navparse.hh"
using namespace std;
struct ObserverPosition
@ -36,31 +37,6 @@ struct ObserverPosition
};
std::map<int, ObserverPosition> g_srcpos;
template<typename T>
class GetterSetter
{
public:
void set(const T& t)
{
std::lock_guard<std::mutex> mut(d_lock);
d_t = t;
}
T get()
{
T ret;
{
std::lock_guard<std::mutex> mut(d_lock);
ret = d_t;
}
return ret;
}
private:
T d_t;
std::mutex d_lock;
};
map<int, BeidouAlmanacEntry> g_beidoualma;
map<int, pair<GlonassMessage, GlonassMessage>> g_glonassalma;
@ -130,64 +106,6 @@ string humanUra(uint8_t ura)
}
struct SVIOD
{
std::bitset<32> words;
int gnssid;
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};
// 60 seconds
uint16_t t0c; // clock epoch, stored UNSCALED, since it is not in the same place as GPS
// 2^-34 2^-46
int32_t af0{0} , af1{0};
// 2^-59
int8_t af2{0};
uint8_t sisa;
uint32_t wn{0}, tow{0};
bool complete() const
{
if(gnssid==2)
return words[1] && words[2] && words[3] && words[4];
else
return words[2] && words[3];
}
void addGalileoWord(std::basic_string_view<uint8_t> page);
double getMu() const
{
if(gnssid == 2) return 3.986004418 * pow(10.0, 14.0);
if(gnssid == 0) return 3.986005 * pow(10.0, 14.0);
throw std::runtime_error("getMu() called for unsupported gnssid "+to_string(gnssid));
} // 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; }
uint32_t getT0c() const { return 60*t0c; }
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
};
void SVIOD::addGalileoWord(std::basic_string_view<uint8_t> page)
{
@ -231,78 +149,7 @@ void SVIOD::addGalileoWord(std::basic_string_view<uint8_t> page)
}
struct SVPerRecv
{
int el{-1}, azi{-1}, db{-1};
time_t deltaHzTime{-1};
double deltaHz{-1};
double prres{-1};
time_t t; // last seen
};
/* Most of thes fields are raw, EXCEPT:
t0t = seconds, raw fields are 3600 seconds (galileo), 4096 seconds (GPS)
*/
struct SVStat
{
uint8_t e5bhs{0}, e1bhs{0};
uint8_t gpshealth{0};
uint16_t ai0{0};
int16_t ai1{0}, ai2{0};
bool sf1{0}, sf2{0}, sf3{0}, sf4{0}, sf5{0};
int BGDE1E5a{0}, BGDE1E5b{0};
bool e5bdvs{false}, e1bdvs{false};
uint16_t wn{0}; // we put the "unrolled" week number here!
uint32_t tow{0}; // "last seen"
//
// 2^-30 2^-50 1 8-bit week
int32_t a0{0}, a1{0}, t0t{0}, wn0t{0};
int32_t a0g{0}, a1g{0}, t0g{0}, wn0g{0};
int8_t dtLS{0}, dtLSF{0};
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, 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;
BeidouMessage lastBeidouMessage1, lastBeidouMessage2;
TLERepo::Match tleMatch;
double lastTLELookupX{0};
// new galileo
GalileoMessage galmsg;
map<int, GalileoMessage> galmsgTyped;
// Glonass
GlonassMessage glonassMessage;
pair<uint32_t, GlonassMessage> glonassAlmaEven;
map<uint64_t, SVPerRecv> perrecv;
double latestDisco{-1}, latestDiscoAge{-1}, timeDisco{-1000};
map<int, SVIOD> iods;
void addGalileoWord(std::basic_string_view<uint8_t> page);
bool completeIOD() const;
uint16_t getIOD() const;
SVIOD liveIOD() const;
SVIOD& getEph(int i) { return iods[i]; } // XXXX gps adaptor
pair<int,SVIOD> prevIOD{-1, SVIOD()};
void clearPrev()
{
prevIOD.first = -1;
}
void checkCompleteAndClean(int iod);
};
bool SVStat::completeIOD() const
{
@ -409,7 +256,6 @@ void SVStat::addGalileoWord(std::basic_string_view<uint8_t> page)
}
typedef std::map<SatID, SVStat> svstats_t;
svstats_t g_svstats;
GetterSetter<svstats_t> g_statskeeper;
@ -657,7 +503,7 @@ std::optional<double> getHzCorrection(time_t now, int src, unsigned int gnssid,
std::string humanBhs(int bhs)
{
static vector<string> options{"ok", "out of service", "will be out of service", "test"};
if(bhs >= options.size()) {
if(bhs >= (int)options.size()) {
cerr<<"Asked for humanBHS "<<bhs<<endl;
return "??";
}
@ -1121,6 +967,31 @@ try
return ret;
}
);
h2s.addHandler("/galcov.json", [](auto handler, auto req) {
auto cov = emitCoverage();
auto ret = nlohmann::json::array();
// ret =
// [ [90, [[-180, 3], [-179, 3], ... [180,3] ]]
// [89, [[-180, 4], [-179, 4], ... [180,2] ]]
// ]
for(const auto& latvect : cov) {
auto jslatvect = nlohmann::json::array();
auto jslongvect = nlohmann::json::array();
for(const auto& longpair : latvect.second) {
auto jsdatum = nlohmann::json::array();
jsdatum.push_back((int)longpair.first);
jsdatum.push_back(longpair.second);
jslongvect.push_back(jsdatum);
}
jslatvect.push_back(latvect.first);
jslatvect.push_back(jslongvect);
ret.push_back(jslatvect);
}
return ret;
});
h2s.addHandler("/svs.json", [](auto handler, auto req) {
auto svstats = g_statskeeper.get();

147
navparse.hh 100644
View File

@ -0,0 +1,147 @@
#include "navmon.hh"
#include "galileo.hh"
#include "gps.hh"
#include "beidou.hh"
#include "glonass.hh"
#include <map>
#include "tle.hh"
using namespace std; // XXX
struct SVPerRecv
{
int el{-1}, azi{-1}, db{-1};
time_t deltaHzTime{-1};
double deltaHz{-1};
double prres{-1};
time_t t; // last seen
};
struct SVIOD
{
std::bitset<32> words;
int gnssid;
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};
// 60 seconds
uint16_t t0c; // clock epoch, stored UNSCALED, since it is not in the same place as GPS
// 2^-34 2^-46
int32_t af0{0} , af1{0};
// 2^-59
int8_t af2{0};
uint8_t sisa;
uint32_t wn{0}, tow{0};
bool complete() const
{
if(gnssid==2)
return words[1] && words[2] && words[3] && words[4];
else
return words[2] && words[3];
}
void addGalileoWord(std::basic_string_view<uint8_t> page);
double getMu() const
{
if(gnssid == 2) return 3.986004418 * pow(10.0, 14.0);
if(gnssid == 0) return 3.986005 * pow(10.0, 14.0);
throw std::runtime_error("getMu() called for unsupported gnssid "+to_string(gnssid));
} // 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; }
uint32_t getT0c() const { return 60*t0c; }
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
};
/* Most of thes fields are raw, EXCEPT:
t0t = seconds, raw fields are 3600 seconds (galileo), 4096 seconds (GPS)
*/
struct SVStat
{
uint8_t e5bhs{0}, e1bhs{0};
uint8_t gpshealth{0};
uint16_t ai0{0};
int16_t ai1{0}, ai2{0};
bool sf1{0}, sf2{0}, sf3{0}, sf4{0}, sf5{0};
int BGDE1E5a{0}, BGDE1E5b{0};
bool e5bdvs{false}, e1bdvs{false};
uint16_t wn{0}; // we put the "unrolled" week number here!
uint32_t tow{0}; // "last seen"
//
// 2^-30 2^-50 1 8-bit week
int32_t a0{0}, a1{0}, t0t{0}, wn0t{0};
int32_t a0g{0}, a1g{0}, t0g{0}, wn0g{0};
int8_t dtLS{0}, dtLSF{0};
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, 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;
BeidouMessage lastBeidouMessage1, lastBeidouMessage2;
TLERepo::Match tleMatch;
double lastTLELookupX{0};
// new galileo
GalileoMessage galmsg;
map<int, GalileoMessage> galmsgTyped;
// Glonass
GlonassMessage glonassMessage;
pair<uint32_t, GlonassMessage> glonassAlmaEven;
map<uint64_t, SVPerRecv> perrecv;
double latestDisco{-1}, latestDiscoAge{-1}, timeDisco{-1000};
map<int, SVIOD> iods;
void addGalileoWord(std::basic_string_view<uint8_t> page);
bool completeIOD() const;
uint16_t getIOD() const;
SVIOD liveIOD() const;
SVIOD& getEph(int i) { return iods[i]; } // XXXX gps adaptor
pair<int,SVIOD> prevIOD{-1, SVIOD()};
void clearPrev()
{
prevIOD.first = -1;
}
void checkCompleteAndClean(int iod);
};
typedef std::map<SatID, SVStat> svstats_t;
typedef vector<pair<double,vector<pair<double, int> > > > covmap_t;
covmap_t emitCoverage();