spacecruft
/
galmon
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

add coverage map with hdop/pdop/vdop

experimental
bert hubert 2019-10-23 10:33:13 +02:00
parent 003164e516
commit 4534472f5e
3 changed files with 136 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
coverage.cc
View File

@ -40,40 +40,47 @@ xDOP getDOP(Point& us, vector<Point> sats)
// cout<<G<<endl; // cout<<G<<endl;
MatrixXd Q = (G.transpose() * G).inverse(); MatrixXd Q = (G.transpose() * G).inverse();
auto [lambda, phi] = getLongLat(us.x, us.y, us.z);
// https://gssc.esa.int/navipedia/index.php/Positioning_Error
Eigen::Matrix3d Renu;
Renu <<
(-sin(lambda)) , (-sin(phi)*cos(lambda)) , (cos(phi)*cos(lambda)),
(cos(lambda)) , (-sin(phi)*sin(lambda)) , (cos(phi)*sin(lambda)),
(0.0) , (cos(phi)) , (sin(phi));
Eigen::Matrix3d Qxyz;
for(int x=0; x<3; ++x) // feels like there should be a better way for this, but not sure
for(int y=0; y<3; ++y)
Qxyz(x,y) = Q(x,y);
Eigen::Matrix3d Qenu = Renu.transpose()*Qxyz*Renu;
// if(Qenu(0,0) < 0 || Qenu(1,1) < 0 || Qenu(2,2) < 0)
// cout << "Original: \n"<<Qxyz<<"\nRotated: \n"<<Qenu<<endl;
ret.pdop = sqrt(Q(0,0) + Q(1,1) + Q(2,2)); // already squared ret.pdop = sqrt(Q(0,0) + Q(1,1) + Q(2,2)); // already squared
// ret.pdop = sqrt(Qenu(0,0) + Qenu(1,1) + Qenu(2,2)); // already squared
ret.tdop = sqrt(Q(3,3)); ret.tdop = sqrt(Q(3,3));
ret.gdop = sqrt(ret.pdop*ret.pdop + ret.tdop*ret.tdop); ret.gdop = sqrt(ret.pdop*ret.pdop + ret.tdop*ret.tdop);
if(Qenu(0,0) >=0 && Qenu(1,1) >=0)
ret.hdop = sqrt(Qenu(0,0) + Qenu(1,1));
if(Qenu(2,2)>=0)
ret.vdop = sqrt(Qenu(2,2));
return ret; return ret;
}; };
covmap_t emitCoverage() covmap_t emitCoverage(const vector<Point>& sats)
{ {
covmap_t ret; covmap_t ret;
ofstream cmap("covmap.csv"); ofstream cmap("covmap.csv");
cmap<<"latitude longitude count5 count10 count20"<<endl; 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,(uint32_t)g.first,1}].completeIOD() && svstats[{2,(uint32_t)g.first,1}].liveIOD().sisa == 255) {
// cout<<g.first<<" NAPA!"<<endl;
continue;
}
sats[g.first]=sat;
}
double R = 6371000; double R = 6371000;
for(double latitude = 90 ; latitude > -90; latitude-=0.5) { // north-south for(double latitude = 90 ; latitude > -90; latitude-=2) { // north-south
double phi = M_PI* latitude / 180; double phi = M_PI* latitude / 180;
double longsteps = 1 + 360.0 * cos(phi); double longsteps = 1 + 360.0 * cos(phi);
double step = 180.0 / longsteps; double step = 4*180.0 / longsteps;
vector<tuple<double, int, int, int, double, double, double>> latvect; vector<tuple<double, int, int, int, double, double, double, double, double, double,double, double, double>> latvect;
for(double longitude = -180; longitude < 180; longitude += step) { // east - west for(double longitude = -180; longitude < 180; longitude += step) { // east - west
Point p; Point p;
// phi = latitude, lambda = longitude // phi = latitude, lambda = longitude
@ -93,17 +100,17 @@ covmap_t emitCoverage()
vector<Point> satposs5, satposs10, satposs20; vector<Point> satposs5, satposs10, satposs20;
for(const auto& s : sats) { for(const auto& s : sats) {
// double getElevationDeg(const Point& sat, const Point& our); // double getElevationDeg(const Point& sat, const Point& our);
double elev = getElevationDeg(s.second, p); double elev = getElevationDeg(s, p);
if(elev > 5.0) { if(elev > 5.0) {
satposs5.push_back(s.second); satposs5.push_back(s);
numsats5++; numsats5++;
} }
if(elev > 10.0) { if(elev > 10.0) {
satposs10.push_back(s.second); satposs10.push_back(s);
numsats10++; numsats10++;
} }
if(elev > 20.0) { if(elev > 20.0) {
satposs20.push_back(s.second); satposs20.push_back(s);
numsats20++; numsats20++;
} }
} }
@ -111,7 +118,14 @@ covmap_t emitCoverage()
numsats5, numsats10, numsats20, numsats5, numsats10, numsats20,
getDOP(p, satposs5).pdop, getDOP(p, satposs5).pdop,
getDOP(p, satposs10).pdop, getDOP(p, satposs10).pdop,
getDOP(p, satposs20).pdop getDOP(p, satposs20).pdop,
getDOP(p, satposs5).hdop,
getDOP(p, satposs10).hdop,
getDOP(p, satposs20).hdop,
getDOP(p, satposs5).vdop,
getDOP(p, satposs10).vdop,
getDOP(p, satposs20).vdop
)); ));
// cmap << longitude <<" " <<latitude <<" " << numsats5 << " " <<numsats10<<" "<<numsats20<<endl; // cmap << longitude <<" " <<latitude <<" " << numsats5 << " " <<numsats10<<" "<<numsats20<<endl;
} }

94
navparse.cc
View File

@ -1016,8 +1016,88 @@ try
} }
); );
h2s.addHandler("/galcov.json", [](auto handler, auto req) { h2s.addHandler("/cov.json", [](auto handler, auto req) {
auto cov = emitCoverage(); vector<Point> sats;
auto galileoalma = g_galileoalmakeeper.get();
auto gpsalma = g_gpsalmakeeper.get();
auto beidoualma = g_beidoualmakeeper.get();
auto svstats = g_statskeeper.get();
// cout<<"pseudoTow "<<pseudoTow<<endl;
string_view path = convert(req->path);
bool doGalileo{true}, doGPS{false}, doBeidou{false};
auto pos = path.find("gps=");
if(pos != string::npos) {
doGPS = (path[pos+4]=='1');
}
pos = path.find("galileo=");
if(pos != string::npos) {
doGalileo = (path[pos+8]=='1');
}
pos = path.find("beidou=");
if(pos != string::npos) {
doBeidou = (path[pos+7]=='1');
}
if(doGalileo)
for(const auto &g : galileoalma) {
Point sat;
getCoordinates(latestTow(2, svstats), g.second, &sat);
if(g.first < 0)
continue;
SatID id{2,(uint32_t)g.first,1};
if(svstats[id].completeIOD() && svstats[id].liveIOD().sisa == 255) {
continue;
}
if(svstats[id].e1bhs || svstats[id].e1bdvs)
continue;
sats.push_back(sat);
}
if(doGPS)
for(const auto &g : gpsalma) {
Point sat;
getCoordinates(latestTow(0, svstats), g.second, &sat);
if(g.first < 0)
continue;
SatID id{0,(uint32_t)g.first,0};
if(svstats[id].completeIOD() && svstats[id].ura == 16) {
// cout<<"Skipping G"<<id.sv<<" because of URA"<<endl;
continue;
}
if(svstats[id].gpshealth) {
// cout<<"Skipping G"<<id.sv<<" because of health"<<endl;
continue;
}
sats.push_back(sat);
}
if(doBeidou)
for(const auto &g : beidoualma) {
Point sat;
getCoordinates(latestTow(3, svstats), g.second.alma, &sat);
if(g.first < 0)
continue;
SatID id{3,(uint32_t)g.first,0};
/*
if(svstats[id].completeIOD() && svstats[id].ura == 16) {
cout<<"Skipping G"<<id.sv<<" because of URA"<<endl;
continue;
}
*/
if(svstats[id].gpshealth) {
// cout<<"Skipping C"<<id.sv<<" because of health"<<endl;
continue;
}
sats.push_back(sat);
}
auto cov = emitCoverage(sats);
auto ret = nlohmann::json::array(); auto ret = nlohmann::json::array();
// ret = // ret =
@ -1037,6 +1117,16 @@ try
jsdatum.push_back((int)(10*get<4>(longpair))); jsdatum.push_back((int)(10*get<4>(longpair)));
jsdatum.push_back((int)(10*get<5>(longpair))); jsdatum.push_back((int)(10*get<5>(longpair)));
jsdatum.push_back((int)(10*get<6>(longpair))); jsdatum.push_back((int)(10*get<6>(longpair)));
jsdatum.push_back((int)(10*get<7>(longpair)));
jsdatum.push_back((int)(10*get<8>(longpair)));
jsdatum.push_back((int)(10*get<9>(longpair)));
jsdatum.push_back((int)(10*get<10>(longpair)));
jsdatum.push_back((int)(10*get<11>(longpair)));
jsdatum.push_back((int)(10*get<12>(longpair)));
jslongvect.push_back(jsdatum); jslongvect.push_back(jsdatum);
} }
jslatvect.push_back(latvect.first); jslatvect.push_back(latvect.first);

9
navparse.hh
View File

@ -143,17 +143,16 @@ struct SVStat
typedef std::map<SatID, SVStat> svstats_t; typedef std::map<SatID, SVStat> svstats_t;
// a vector of pairs of latidude,vector<longitude,numsats> // a vector of pairs of latidude,vector<longitude,numsats>
typedef vector<pair<double,vector<tuple<double, int, int, int, double, double, double> > > > covmap_t; typedef vector<pair<double,vector<tuple<double, int, int, int, double, double, double, double, double, double, double, double, double> > > > covmap_t;
covmap_t emitCoverage();
covmap_t emitCoverage(const vector<Point>& sats);
struct xDOP struct xDOP
{ {
double gdop{-1}; double gdop{-1};
double pdop{-1}; double pdop{-1};
double tdop{-1}; double tdop{-1};
// double hdop; double hdop{-1};
// double vdop; double vdop{-1};
}; };
xDOP getDOP(Point& us, vector<Point> sats); xDOP getDOP(Point& us, vector<Point> sats);