How to calculate position of GPS Satellite with transmitted ephemeris data?

How can I calculate the position of a satellite using the ephemeris data transmitted by each satellite?

I will then use positions generated from each GPS satellite's self-reported ephemeris to calculate the distance to a certain location on earth as discussed in my related question Skyfiled | How to satellite-to-site range given site coordinates and satellite ephemeris?

I've searched and read through a lot of material on GPS satellites but could not find anything that specifically explaina how to use the satellite ephemeris encoded in the data stream into a position in a standard coordinate system.

• Andrew Holmes has built a homebrew GPS receiver -- hardware and software: extremetech.com/extreme/… -- I don't know if he's posted his software or not, and it seems his website has disappeared. Presumably he found an answer to all your questions. Now aholme.co.uk – Bit Chaser Oct 28 '19 at 23:53
• Thx, for your answer, he did post the source code which is brilliant work but does not contain the calculations of the satellite position although he does mention that it can be done. – cy8berpunk Oct 29 '19 at 16:50
• Somewhere between ephemeris.cpp and solve.cpp that's calculated in his code. I suggest finding the equations on wiki, and hope you can find a concise description in documentation. I'll bet calculate current position from TLE is answered somewhere on this site. – Bit Chaser Oct 29 '19 at 17:03
• oh, I over read it, thank you very much! – cy8berpunk Oct 29 '19 at 17:18

You can find a C++ method contained in the source code of Andrew Holme homemade GPS receiver project. The method is called GetXYZ and is the EPHEM (ephemeris) namespace and looks as the following:

void EPHEM::GetXYZ(double *x, double *y, double *z, double t) { // Get satellite position at time t

// Time from ephemeris reference epoch
double t_k = TimeFromEpoch(t, t_oe);

// Eccentric Anomaly
double E_k = EccentricAnomaly(t_k);

// True Anomaly
double v_k = atan2(
sqrt(1-e*e) * sin(E_k),
cos(E_k) - e);

// Argument of Latitude
double AOL = v_k + omega;

// Second Harmonic Perturbations
double du_k = C_us*sin(2*AOL) + C_uc*cos(2*AOL);    // Argument of Latitude Correction
double dr_k = C_rs*sin(2*AOL) + C_rc*cos(2*AOL);    // Radius Correction
double di_k = C_is*sin(2*AOL) + C_ic*cos(2*AOL);    // Inclination Correction

// Corrected Argument of Latitude; Radius & Inclination
double u_k = AOL + du_k;
double r_k = A*(1-e*cos(E_k)) + dr_k;
double i_k = i_0 + di_k + IDOT*t_k;

// Positions in orbital plane
double x_kp = r_k*cos(u_k);
double y_kp = r_k*sin(u_k);

// Corrected longitude of ascending node
double OMEGA_k = OMEGA_0 + (OMEGA_dot-OMEGA_E)*t_k - OMEGA_E*t_oe;

// Earth-fixed coordinates
*x = x_kp*cos(OMEGA_k) - y_kp*cos(i_k)*sin(OMEGA_k);
*y = x_kp*sin(OMEGA_k) + y_kp*cos(i_k)*cos(OMEGA_k);
*z = y_kp*sin(i_k);


}

Thx to @bitchaser!

The required equations and much more are contained in the GPS Interface Control Document (IS-GPS-200L Table 20-IV. Broadcast Navigation User Equations), available on the GPS Website:

https://www.gps.gov/technical/icwg/

This is the source of the equations given above by cy8berpunk