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

Add support for GLONASS ephemeris integration. 

* glonass.hh (getCoordinates): Specialization of getCoordinates from ephemeris.hh.
* glonass.cc (ff, rk4step, getCoordinates): New functions.
* navparse.cc (almanac.json and cov.json): Compute Glonass ECEF coords using the above.
* html/geo/coverage.{html, js}: Add Glonass.
pull/82/head
Serge Belyshev 2020-01-12 21:32:07 +03:00
parent 1c50a76b58
commit 9c6897fcec
5 changed files with 98 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
glonass.cc
View File

@ -2,6 +2,11 @@
#include <math.h>
#include <string.h>
static const double ae = 6378.136; // km // IERS: 6378136.6
static const double mu = 398.6004418E3; // km3/s2 // IERS: 3.986004418
static const double J2 = 1082625.75E-9; // IERS: 1.0826359
static const double oe = 7.2921151467E-5; // rad/s // IERS: 7.292115
// this strips out spare bits + parity, and leaves 10 clean 24 bit words
std::basic_string<uint8_t> getGlonassMessage(std::basic_string_view<uint8_t> payload)
{
@ -47,3 +52,69 @@ uint32_t getGlonassT0e(time_t referencetime, int Tb)
tm.tm_sec = 0;
return timegm(&tm)-3*3600; // and back to UTC
}
// y[0] .. y[2] --> X, Y, Z
// y[3] .. y[5] --> Vx, Vy, Vz
static void ff (const double A[3], const double y[6], double f[6])
{
double X = y[0];
double Y = y[1];
double Z = y[2];
double Vx = y[3];
double Vy = y[4];
double Vz = y[5];
double R = sqrt (X*X + Y*Y + Z*Z);
double C0 = pow(R, -3);
double C1 = 1.5*J2*ae*ae*pow(R, -5);
double C2 = 5*pow(Z/R, 2);
double CXY = C0 + C1 * (1 - C2);
double CZ = C0 + C1 * (3 - C2);
f[0] = Vx;
f[1] = Vy;
f[2] = Vz;
f[3] = (- mu*CXY + oe*oe)*X + 2*oe*Vy + A[0];
f[4] = (- mu*CXY + oe*oe)*Y - 2*oe*Vx + A[1];
f[5] = - mu*CZ*Z + A[2];
}
static void rk4step (const double A[3], double y[6], double h)
{
double k1[6], k2[6], k3[6], k4[6], z[6];
ff (A, y, k1);
for (int j = 0; j < 6; j ++)
z[j] = y[j] + 0.5*h*k1[j];
ff (A, z, k2);
for (int j = 0; j < 6; j ++)
z[j] = y[j] + 0.5*h*k2[j];
ff (A, z, k3);
for (int j = 0; j < 6; j ++)
z[j] = y[j] + h*k3[j];
ff (A, z, k4);
for (int j = 0; j < 6; j ++)
y[j] = y[j] + h * (k1[j] + 2*(k2[j] + k3[j]) + k4[j]) / 6;
}
double getCoordinates(double tow, const GlonassMessage& eph, Point* p)
{
double y0[6] = {ldexp(eph.x, -11), ldexp(eph.y, -11), ldexp(eph.z, -11),
ldexp(eph.dx, -20), ldexp(eph.dy, -20), ldexp(eph.dz, -20)};
double A[3] = {ldexp(eph.ddx, -30), ldexp(eph.ddy, -30), ldexp(eph.ddz, -30)};
// These all fake
uint32_t glotime = eph.getGloTime();
uint32_t gloT0e = getGlonassT0e(glotime + 820368000, eph.Tb);
uint32_t ephtow = (gloT0e - 820368000) % (7*86400);
double DT = tow - ephtow;
int n = abs (DT / 100); // integrate in roughly 100 second steps
double h = DT / n;
for (int j = 0; j < n; j ++)
rk4step (A, y0, h);
*p = Point (1E3*y0[0], 1E3*y0[1], 1E3*y0 [2]);
return 0;
}

2
glonass.hh
View File

@ -4,6 +4,7 @@
#include "bits.hh"
#include <iostream>
#include <math.h>
#include "minivec.hh"
std::basic_string<uint8_t> getGlonassessage(std::basic_string_view<uint8_t> payload);
struct GlonassMessage
@ -195,3 +196,4 @@ struct GlonassMessage
};
uint32_t getGlonassT0e(time_t referencetime, int Tb);
double getCoordinates(double tow, const GlonassMessage& eph, Point* p);

1
html/geo/coverage.html
View File

@ -26,6 +26,7 @@
<input type="checkbox" id="GalE1" onclick="do_timer();" checked> <label for="GalE1">Galileo E1</label> &nbsp;&nbsp;
<input type="checkbox" id="GPSL1CA" onclick="do_timer();"> <label for="GPSL1CA">GPS L1C/A</label> &nbsp;&nbsp;
<input type="checkbox" id="Beidou" onclick="do_timer();"> <label for="Beidou">Beidou</label> &nbsp;&nbsp;
<input type="checkbox" id="Glonass" onclick="do_timer();"> <label for="GloL1OF">Glonass</label> &nbsp;&nbsp;
</p>
</div>
</div>

7
html/geo/coverage.js
View File

@ -1,7 +1,7 @@
'use strict';
//
var wantGPS=0, wantGalileo=0, wantBeidou=0;
var wantGPS=0, wantGalileo=0, wantBeidou=0, wantGlonass=0;
//
//
@ -315,14 +315,17 @@ function do_timer()
wantGPS = 0;
wantGalileo = 0;
wantBeidou = 0;
wantGlonass = 0;
if(d3.select("#GPSL1CA").property("checked"))
wantGPS=1;
if(d3.select("#GalE1").property("checked"))
wantGalileo=1;
if(d3.select("#Beidou").property("checked"))
wantBeidou=1;
if(d3.select("#Glonass").property("checked"))
wantGlonass=1;
var galcovurl="../cov.json?gps="+wantGPS+"&galileo="+wantGalileo+"&beidou="+wantBeidou;
var galcovurl="../cov.json?gps="+wantGPS+"&galileo="+wantGalileo+"&beidou="+wantBeidou+"&glonass="+wantGlonass;
d3.queue(1)
.defer(d3.json, fileAlmanac)

23
navparse.cc
View File

@ -798,9 +798,7 @@ try
getCoordinates(latestTow(sv.first.gnss, svstats), sv.second.oldBeidouMessage, &sat);
}
if(sv.first.gnss == 6) {
sat.x = sv.second.glonassMessage.x;
sat.y = sv.second.glonassMessage.y;
sat.z = sv.second.glonassMessage.z;
getCoordinates(latestTow(6, svstats), sv.second.glonassMessage, &sat);
}
if(sat.x) {
Point our = g_srcpos[iter->first].pos;
@ -990,7 +988,7 @@ try
// cout<<"pseudoTow "<<pseudoTow<<endl;
string_view path = convert(req->path);
bool doGalileo{true}, doGPS{false}, doBeidou{false};
bool doGalileo{true}, doGPS{false}, doBeidou{false}, doGlonass{false};
auto pos = path.find("gps=");
if(pos != string::npos) {
doGPS = (path[pos+4]=='1');
@ -1003,6 +1001,10 @@ try
if(pos != string::npos) {
doBeidou = (path[pos+7]=='1');
}
pos = path.find("glonass=");
if(pos != string::npos) {
doGlonass = (path[pos+8]=='1');
}
if(doGalileo)
for(const auto &g : galileoalma) {
@ -1060,6 +1062,19 @@ try
sats.push_back(sat);
}
// fake almanac from svstats -- to avoid ejecting Glonasses to another galaxy due
// to excessive extrapolation, they freeze where they were last seen.
if(doGlonass)
for(const auto &s : svstats) {
if(s.first.gnss == 6 && s.second.wn > 0) {
Point sat;
getCoordinates(s.second.tow, s.second.glonassMessage, &sat);
if(svstats[s.first].glonassMessage.Bn & 7) {
continue;
}
sats.push_back(sat);
}
}
auto cov = emitCoverage(sats);