Galileo Notes: Ephemeris, Anomalies, I/NAV, SFRBX

I’ve recently been spending quite some time on the EU’s Galileo Navigation Satellite System. In this post you’ll find some rough notes on things I’ve found out, both about how to receive data & how to understand the Galileo ephemeris parameters.

Feedback is very welcome on bert@hubertnet.nl or @bert_hu_bert! I’m no Galileo specialist, so I am sure this page could be better. Please also let me know if this page has been useful to you.

In learning about Galileo I benefited greatly from reading the source of RTKLIB which provided key insights on both Galileo and the u-blox 8 chipset.

Receiver

I am using the Navilock NL-8012U Multi GNSS Receiver. This has the fabulous u-blox 8 chipset which is very well documented.

I think almost any device with a u-blox 8 chipset should be fine for your needs.

U-blox ships a rather rough tool called u-center to configure the device. This tool is Windows-only, but runs under Wine if you setup Wine to emulate Windows Vista. The USB serial port showed up for me as COM33, and then things worked. Control-F9 is where you send ‘messages’ to configure the device.

The chip speaks both the NMEA protocol and its own binary UBX format. To get Galileo data, you must specifically enable Galileo, because out of the box it will stick to GPS and GLONASS.

In addition, before you’ll see Galileo specific data in NMEA, you’ll have to enable NMEA 4.1.

To get raw Galileo (and GPS, Glonass, BeiDou and other) frames, enable the UBX-RXM-SFRBX messages via u-center.

U-blox UBX-RXM-SFRBX messages to I/Nav

This.. is non-trivial. Galileo transmits navigation frames and subframes and words and pages, and parts of these are sent in reverse order. Wordswapping is also going on. Other receivers will do this differently, but it will always be very fiddly work.

The u-blox procedure is first to reverse the words in the SFRBX message:

std::basic_string<uint8_t> payload;
for(unsigned int i = 0 ; i < (msg.size() - 8) / 4; ++i)
	for(int j=1; j <= 4; ++j)
		payload.append(1, msg[8 + (i+1) * 4 -j]);

So this replaces ‘12345678’ by ‘43218765’. Ok.

Then from ‘payload’ we need to lift the actual I/NAV contents. These are hidden in two parts, and are not on whole byte boundaries.

std::basic_string<uint8_t> inav;
unsigned int i,j;
for (i=0,j=2; i<14; i++, j+=8)
	inav.append(1, (unsigned char)getbitu(payload.c_str()   ,j,8));
for (i=0,j=2; i< 2; i++, j+=8)
	inav.append(1, (unsigned char)getbitu(payload.c_str()+16,j,8));

In inav, we now have a ready to parse I/NAV message, as described in 4.3.5 of the Galileo - Open Service - Signal In Space Interface Control Document (OS SIS ICD V1.3).

Incidentally, here are getbitu and getbits from RTKLIB:

unsigned int getbitu(const unsigned char *buff, int pos, int len)
{
  unsigned int bits=0;
  int i;
  for (i=pos;i<pos+len;i++) bits=(bits<<1)+((buff[i/8]>>(7-i%8))&1u);
  return bits;
}

int getbits(const unsigned char *buff, int pos, int len)
{
    unsigned int bits=getbitu(buff,pos,len);
    if (len<=0||32<=len||!(bits&(1u<<(len-1)))) return (int)bits;
    return (int)(bits|(~0u<<len)); /* extend sign */
}

Ephemeris

Once you have extracted the I/NAV data from the complex UBX-RXM-SFRBX messages, you can get all kinds of raw information about the Galileo satellites.

The hardest part to understand however is the coolest part: the detailed ephemerides that tell you where the satellite vehicles are.

All details are in the SIS ICD document, but some things are not obvious, at least to this layperson.

The ephemeris parameters described in table 57 are what you need to know. I wasted half a day over three things:

Semi-circles?

First, $M_0$, $\Delta{}n$, $\Omega_0$, $i_0$, $\omega$, $\dot \Omega$, $\dot i$ are specified in ‘semi-circle’ units. A footnote tells us that a semi-circle is $\pi$ radians, so it seems like we are done. But that is not that they mean - all these parameters need to be multiplied by $\pi$ before you can do any calculations!

Eccentric anomaly iteration

Second, calculating the eccentric anomaly from table 58. This helpfully says this may be done by iteration, and this is true. To save you from a wonderful trip through Numerical Recipes to find an interesting way to get to E, it turns out the algorithm is trivial:

double M = m0 + n * tk;
double E = M;
for(int k =0 ; k < 10; ++k) {
    E = M + e * sin(E);
}

You could be smarter about it and not just do 10 iterations, but perhaps terminate when old E and new E are identical. The Wikipedia pages on the eccentric anomaly and on Kepler’s equation are a fun read, where we find that the algorithm above stems from 1621.

True anomaly

The true anomaly, $\nu$ can be calculated in various ways. The way found in table 58 of the SIS is a lot of typing and (at least for me) gets the sign wrong from time to time.

On the Wikipedia true anomaly page, we find several other formulas to calculate $\nu$, the easiest of which is:

$$ \nu = 2\ \arctan \bigg(\sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}\bigg) $$

This version however also sometimes gets the wrong sign, or is off by $\pi$. The reason for this is likely in how $\tan^{-1}$ operates, where it will only return values between $-\frac{1}{2}\pi$ and $\frac{1}{2}\pi$, even though this result may not be ‘physical’, much like $x^2 = 9$ has -3 and 3 as solutions, even though -3 might not be possible given the physical context.

The Galileo orbital and technical parameters page helpfully lists yet a third way to calculate $\nu$:

$$ \nu = M + e\ 2\ \sin(M) + e^2\frac{5}{4}\sin(2M) - e^3\big(\frac{1}{4}\sin(M)-\frac{13}{12}\sin(3M)\big) + \cdots $$

This one has no inverse operations & consistently gets the right results for me. Note that this function is an approximation that works really well for circular orbits, where $e$ is almost zero, which is true for most Galileo orbits. There are two Galileo satellites in elliptical orbits where $e = 0.162$, for which the ESA supplied formula is not quite precise enough.

On the Wikipedia page for this equation of the center, we find that the following terms are:

$$ + \frac{e^4}{96}(103\sin(4M) - 44\sin(2M)) $$ $$ + \frac{e^5}{960}(1097\sin(5M) - 645\sin(3M) + 50\sin(M)) $$ $$ + \frac{e^6}{960}(1223\sin(6M) - 902\sin(4M) + 85\sin(2M))+ \cdots $$

The Galileo eccentric satellites do transmit ephemerides, but these are marked as ’testing’, causing receivers to not use these two vehicles for position determination. I think this is wise because not all receivers might be equipped with the math to deal with really eccentric orbits.

Extracting the parameters from the I/NAV messages

To save you a lot of error prone typing, here is what works for me:

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};
  
uint8_t sisa;

void SVIOD::addWord(std::basic_string_view<uint8_t> page)
{
  uint8_t wtype = getbitu(&page[0], 0, 6);
  if(wtype == 1) {
    t0e = getbitu(&page[0],   16, 14);
    m0 = getbits(&page[0],    30, 32);
    e = getbitu(&page[0],     62, 32);
    sqrtA = getbitu(&page[0], 94, 32);
  }
  else if(wtype == 2) {
    omega0 = getbits(&page[0], 16, 32);
    i0 = getbits(&page[0],     48, 32);
    omega = getbits(&page[0],  80, 32);
    idot = getbits(&page[0],   112, 14);
  }
  else if(wtype == 3) {
    omegadot = getbits(&page[0], 16, 24);
    deltan = getbits(&page[0],   40, 16);
    cuc = getbits(&page[0],      56, 16);
    cus = getbits(&page[0],      72, 16);
    crc = getbits(&page[0],      88, 16);
    crs = getbits(&page[0],     104, 16);
    sisa = getbitu(&page[0],    120, 8);
  }
  else if(wtype == 4) {
    cic = getbits(&page[0], 22, 16);
    cis = getbits(&page[0], 38, 16);
  }
}

Bringing it all together

I can’t yet vouch for all this code being correct. Do not use it to target your new ICBM in any case! But it does generate plausible results. This code assumes you’ve copied all the parameters from the I/NAV frames, and that you’ve done so correctly, including paying attention to the sign-bits for the fields that are two’s complement encoded:

void getCoordinates(int wn, int tow, const SVIOD& iod, 
		    double* x, double* y, double* z)
{
  constexpr double mu = 3.986004418 * pow(10.0, 14.0);
  constexpr double omegaE = 7.2921151467 * pow(10.0, -5);

  // iod contains the raw I/NAV integers
  
  double sqrtA = 1.0*iod.sqrtA / (1ULL<<19);
  double deltan = M_PI * 1.0*iod.deltan / (1ULL<<43);
  double t0e = 60.0*iod.t0e;
  double m0 = M_PI * 1.0*iod.m0 / (1ULL<<31);
  double e = 1.0*iod.e / (1ULL<<33);
  double omega = M_PI * 1.0*iod.omega / (1ULL<<31);

  double cuc = 1.0*iod.cuc / (1ULL<<29);
  double cus = 1.0*iod.cus / (1ULL<<29);

  double crc = 1.0*iod.crc / (1ULL<<5);
  double crs = 1.0*iod.crs / (1ULL<<5);

  double cic = 1.0*iod.cic / (1ULL<<29);
  double cis = 1.0*iod.cis / (1ULL<<29);

  double idot = M_PI * 1.0*iod.idot / (1ULL<<43);
  double i0 = M_PI * 1.0*iod.i0 / (1ULL << 31);

  double Omegadot = M_PI * 1.0*iod.omegadot / (1ULL << 43);
  double Omega0 = M_PI * 1.0*iod.omega0 / (1ULL << 31);
  
  // NO IOD BEYOND THIS POINT!
  
  double A = pow(sqrtA, 2.0);
  double A3 = pow(sqrtA, 6.0);
  double n0 = sqrt(mu/A3);
  double tk = tow - t0e; // in seconds, ignores WN!

  double n = n0 + deltan;
  double M = m0 + n * tk;
  double E = M;
  for(int k =0 ; k < 10; ++k) {
    E = M + e * sin(E);
  }

  double nu =  M + 
    e* 2*sin(M) +
    e*e * (5.0/4.0) * sin(2*M) -
    e*e*e * (0.25*sin(M) - (13.0/12.0)*sin(3*M)) +
    e*e*e*e * (103*sin(4*M) - 44*sin(2*M)) / 96.0 +
    e*e*e*e*e * (1097*sin(5*M) - 645*sin(3*M) + 50 *sin(M))/960.0 +
    e*e*e*e*e*e * (1223*sin(6*M) - 902*sin(4*M) + 85 *sin(2*M))/960.0;
  
  double psi = nu + omega;

  double deltau = cus * sin(2*psi) + cuc * cos(2*psi);
  double deltar = crs * sin(2*psi) + crc * cos(2*psi);
  double deltai = cis * sin(2*psi) + cic * cos(2*psi);

  double u = psi + deltau;
  double r = A * (1- e * cos(E)) + deltar;

  double xprime = r*cos(u), yprime = r*sin(u);
  double Omega = Omega0 + (Omegadot - omegaE)*tk - omegaE * t0e;
  double i = i0 + deltai + idot * tk;

  *x = xprime * cos(Omega) - yprime * cos(i) * sin(Omega);
  *y = xprime * sin(Omega) + yprime * cos(i) * cos(Omega);
  *z = yprime * sin(i);
}

Note that this code will not work across week-number boundaries! This is the line that starts with double tk.

The output of this function corresponds to the algorithm described in table 58 of the Galileo SIS ICD document, and gets you the coordinates in the Earth Centric Earth Facing Galileo Terrestrial Reference Frame, or ECEF GTRF.

X, Y, Z, latitude, longitude

Now that we have the x, y and z coordinates, an easy sanity check is to calculate the elevation of a Galileo satellite over the horizon, not in the least because this number is also easily compared to what popular websites tell us. Many receivers will also inform us of their elevation calculation.

Now, if you are skilled with trigonometry, this is a trivial calculation, not worth writing a whole web page about. It also appears that if you are not skilled with vector calculations or trigonometry, it is very tempting to write up a page that misinforms or confuses readers on how to calculate the elevation of a satellite.

To save you the world of pain I went through as someone very bad at 3D math, here is how to do it.

For starters, the elevation of course depends on where the viewer is. A better name for elevation perhaps is “viewing angle”.

On this helpful page on the GSSC “Navipedia”, we find how the Galileo ECEF coordinates work. The X-axis points straight up if you are standing on the zero meridian on the equator.

The Y-axis does the same if you stand at 90 east on the equator.

Finally, the Z-axis is your distance ‘up’ from the equator, so it goes from the center of the earth to the north pole.

On the Navipedia we can find the math for a flattened ellipsoidal earth, here I’ll show the simpler version where we assume the earth is round, and we are at sea level:

$$ x = R \cos\varphi \cos \lambda $$ $$ y = R \cos\varphi \sin \lambda $$ $$ z = R \sin\varphi $$

$\varphi$ is the latitude, $\lambda$ is the longitude, $R$ is the radius of the earth (6371KM for our purposes).

Nootdorp, The Netherlands is at 52.0440N, 4.3909E (latitude, longitude), which gives us $\varphi=0.90834, \lambda=0.7664$ (in radians). If we fill out the rest, this gives us (3907, 300, 5023) as coordinates, all in kilometers. As a sanity check, you can try this page which offers an online converter.

Determining the elevation

Now that we know our own ECEF coordinates, we can take the ECEF coordinates of a Galileo satellite, and determine the elevation or viewing angle.

The key insight is that what we are looking for is the angle between the vector from the center of the earth to us, and the vector pointing from us to the satellite. If the satellite is exactly overhead from us, this angle is 0 degrees, as these vectors are then parallel.

If the satellite is exactly on our horizon, the angle between the our ‘center of earth’ vector and the vector to the satellite is exactly 90 degrees.

Or in other words, our viewing angle is 90 minus the angle between our vectors.

For any two vectors A and B, the following holds, where $\omega$ is the angle between the vectors:

$$ \cos \omega = \frac{\vec A\cdot{}\vec B}{\lVert \vec A \rVert \lVert \vec B\rVert} $$

This is the inner product of A and B, divided by the product of the lengths of A and B.

A denotes our own ECEF $(x,y,z)$ coordinates, B is $\Delta{}x, \Delta{}y, \Delta{}z$ or $(satx - x, saty - y, satz - z)$, the vector pointing from us to the satellite.

Filling this out gets us:

$$ \cos \omega = \frac{ x\cdot{}\Delta{}x + y\cdot{}\Delta{}y + z\cdot{}\Delta{}z}{\sqrt{x^2+y^2+z^2} \cdot \sqrt{\Delta{}x^2 + \Delta{}y^2 +\Delta{}z^2}} $$

To get the actual viewing angle, take $\cos^{-1}$ of this, convert $\omega$ to degrees, and subtract from 90.

Summarising

I realise that parts of the above are trivial for seasoned professionals, but I hope this will help other people get started. I sincerely hope this page will be helpful.

Feedback is more than welcome on bert@hubertnet.nl or @bert_hu_bert.