Code for downloading, parsing and propagating GNSS Satellites. Uses acfutils and opengpws.

.Space.cpp

#include <string>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <ctime>
#include <iostream>
#include <curl/curl.h>
#include <XPLMProcessing.h>

#include "opengpws/xplane_api.h"
#include "Space.hpp"

extern const egpws_intf_t* egpws_intf;

static bool FileExists(const char* file) {
	std::ifstream infile(file);
	return infile.good();
}

static float FlightLoop(float inElapsedSinceLastCall, float inElapsedTimeSinceLastFlightLoop, int inCounter, void* inRefcon) {
	UNUSED(inElapsedSinceLastCall);
	UNUSED(inElapsedTimeSinceLastFlightLoop);
	UNUSED(inCounter);
	SpaceLoop* space = (SpaceLoop*)inRefcon;

	struct tm t;
	time_t unix_time;
	int days = XPLMGetDatai(space->days_dref);
	t.tm_year = space->year - 1900;
	lacf_yday_to_mon_mday(days, &t.tm_mon, &t.tm_mday);
	t.tm_hour = XPLMGetDatai(space->zhrs_dref);
	t.tm_min = XPLMGetDatai(space->zmins_dref);
	t.tm_sec = XPLMGetDatai(space->zsecs_dref);
	t.tm_isdst = 0;
	unix_time = mktime(&t);
	space->gps_time = space->UnixToGPS(unix_time);

	vect3_t ecef;
#if XPLM400
	XPLMDataRef ecef_x = XPLMFindDataRef("sim/flightmodel2/position/ecef_x");
	XPLMDataRef ecef_y = XPLMFindDataRef("sim/flightmodel2/position/ecef_y");
	XPLMDataRef ecef_z = XPLMFindDataRef("sim/flightmodel2/position/ecef_z");
	ecef = { XPLMGetDatad(ecef_x), XPLMGetDatad(ecef_y), XPLMGetDatad(ecef_z) };
#else
	XPLMDataRef geo_lat = XPLMFindDataRef("sim/flightmodel/position/latitude");
	XPLMDataRef geo_lon = XPLMFindDataRef("sim/flightmodel/position/longitude");
	XPLMDataRef geo_elev = XPLMFindDataRef("sim/flightmodel/position/elevation");
	ecef = geo2ecef_mtr({ XPLMGetDatad(geo_lat), XPLMGetDatad(geo_lon), XPLMGetDatad(geo_elev) }, &wgs84);
#endif

	mutex_enter(&space->prop_mute);
	space->d_time = space->gps_time - space->GNSS.week_epoch;
	space->ecef_pos = ecef;
	mutex_exit(&space->prop_mute);

	return 0.25;
}
static size_t WriteData(void* ptr, size_t size, size_t nmemb, void* stream) {
	std::string data((const char*)ptr, (size_t)size * nmemb);
	*((std::stringstream*)stream) << data;
	return size * nmemb;
}

static bool_t WorkerFunction(void* usrInf) {
	SpaceLoop* space = (SpaceLoop*)usrInf;

	mutex_enter(&space->prop_mute);
	vect3_t ecef = space->ecef_pos;
	time_t delta_time = space->d_time;
	mutex_exit(&space->prop_mute);

	/* Satellite Look angles */
	for (size_t i = 0; i < space->GNSS.satellites.size(); i++) {
		space->CalcSatPos(&space->GNSS.satellites[i], delta_time);
		space->SatLookAngle(&space->GNSS.satellites[i], ecef);
	}
	for (int i = 0; i < 3; i++) {
		space->SatLookAngle(&space->WAAS.satellites[i], ecef);
	}
	for (int i = 0; i < 2; i++) {
		space->SatLookAngle(&space->EGNOS.satellites[i], ecef);
	}
	for (int i = 0; i < 1; i++) {
		space->SatLookAngle(&space->MSAS.satellites[i], ecef);
	}

	/* Satellite Terrain Mask Angles */
	geo_pos3_t obs = ecef2geo(ecef, &wgs84);
	std::vector<Satellite*>list;

	for (size_t i = 0; i < space->GNSS.satellites.size(); i++) {
		list.push_back(&space->GNSS.satellites[i]);
	}
	for (int i = 0; i < 3; i++) {
		list.push_back(&space->WAAS.satellites[i]);
	}
	for (int i = 0; i < 2; i++) {
		list.push_back(&space->EGNOS.satellites[i]);
	}
	for (int i = 0; i < 1; i++) {
		list.push_back(&space->MSAS.satellites[i]);
	}

	for (size_t i = 0; i < list.size(); i++) {
		space->SatMaskAngle(list[i], obs);
	}

	return B_TRUE;
}

/* Public */
SpaceLoop::SpaceLoop(const char* xp, const char* plug) {
	xp_dir = xp;
	plug_dir = plug;

	DownloadAlmanac();
	FindAlmanac();

	mutex_init(&prop_mute);
	zhrs_dref = XPLMFindDataRef("sim/cockpit2/clock_timer/zulu_time_hours");
	zmins_dref = XPLMFindDataRef("sim/cockpit2/clock_timer/zulu_time_minutes");
	zsecs_dref = XPLMFindDataRef("sim/cockpit2/clock_timer/zulu_time_seconds");
	days_dref = XPLMFindDataRef("sim/time/local_date_days");

	XPLMCreateFlightLoop_t params = {
		.structSize = sizeof(XPLMCreateFlightLoop_t),
		.phase = xplm_FlightLoop_Phase_BeforeFlightModel,
		.callbackFunc = FlightLoop,
		.refcon = this
	};
	flight_id = XPLMCreateFlightLoop(&params);
	XPLMScheduleFlightLoop(flight_id, 0.25, false);
	worker_init(&worker, WorkerFunction, 250000, this, "GPS Space Loop Worker");

	Log("SpaceLoop Inited.");
}
SpaceLoop::~SpaceLoop() {
	worker_fini(&worker);
	XPLMDestroyFlightLoop(flight_id);

	zhrs_dref = NULL;
	zmins_dref = NULL;
	zsecs_dref = NULL;
	days_dref = NULL;

	mutex_destroy(&prop_mute);

	Log("SpaceLoop Destroyed.");
}

void SpaceLoop::CalcSatPos(Satellite* sat, time_t delta_time) {
	float GM = 3.986005e14;
	float We = 7.2921151467e-5;

	float A = pow(sat->semi_maj, 2);

	// mean motion
	float n0 = sqrt(GM / pow(A, 3));

	// Mean anomaly
	float Mk = sat->anom + (n0 * delta_time);

	// solve for eccentric anomaly
	float Ek = SolveE(Mk, sat->ecc);
	// True anomaly
	float Vk = 2 * atan(sqrt((1 + sat->ecc) / (1 - sat->ecc)) * tan(Ek / 2));

	// argument of latitude
	float PHIk = Vk + sat->perigee;

	// corrected radius
	float Rk = A * (1 - (sat->ecc * cos(Ek)));

	// corected inclination
	float Inck = 0.94247779608 + sat->incl;

	// position in orbital plane
	float Xprime = Rk * cos(PHIk);
	float Yprime = Rk * sin(PHIk);

	// corrected longitude of ascending node
	float Ok = sat->lng + ((sat->ascen - We) * delta_time) - (We * GNSS.applic_secs);

	// earth fixed geoncentric satellite coordinate
	float x = (Xprime * cos(Ok)) - (Yprime * cos(Inck) * sin(Ok));
	float y = (Xprime * sin(Ok)) + (Yprime * cos(Inck) * cos(Ok));
	float z = Yprime * sin(Inck);
	sat->ecef = { (double)x , (double)y, (double)z };
	sat->llh = ecef2geo(sat->ecef, &wgs84);
}
void SpaceLoop::SatLookAngle(Satellite* sat, vect3_t obs) {
	vect3_t path = vect3_sub(sat->ecef, obs);

	/* Incorporate aircraft orientation into this too */
	vect3_t z = { 0, 0, 1 };
	vect3_t up = vect3_unit(obs, NULL);

	vect3_t east = vect3_unit(vect3_xprod(z, up), NULL);
	vect3_t north = vect3_unit(vect3_xprod(up, east), NULL);

	double e = vect3_dotprod(east, path);
	double n = vect3_dotprod(north, path);
	double u = vect3_dotprod(up, path);

	double r = sqrt(e * e + n * n);
	double az = atan2(e, n);
	double elev = atan2(u, r);

	sat->range = sqrt(pow(path.x, 2) + pow(path.y, 2) + pow(path.z, 2));
	sat->El = elev;
	sat->Az = az;
}
void SpaceLoop::SatMaskAngle(Satellite* sat, geo_pos3_t obs) {
	if (!egpws_intf) {
		sat->terr_mask_angle = 0;
		return;
	}

	geo_pos2_t rec = { obs.lat, obs.lon };
	fpp_t fpp = fpp_init(rec, 0, EARTH_MSL, NULL, B_TRUE);

	enum {
		MAX_PTS = 400,
		SPACING = 100,		        /* meters */
		MIN_DIST = 1000,		    /* meters */
		MAX_DIST = 185200,		    /* meters */
		WATER_OCEAN_MIN = 40000,	/* meters */
		WATER_OCEAN_MAX = 185200	/* meters */
	};

	if (sat->El >= 0) {
		geo_pos2_t* in_pts;
		egpws_terr_probe_t probe = {};

		vect2_t v = geo2fpp({ sat->llh.lat, sat->llh.lon }, &fpp);
		double dist = clamp(gc_distance(rec, { sat->llh.lat, sat->llh.lon }), MIN_DIST, MAX_DIST);
		v = vect2_scmul(vect2_unit(v, NULL), dist);

		probe.num_pts = clampi(dist / SPACING, 2, MAX_PTS);
		in_pts = (geo_pos2_t*)safe_calloc(probe.num_pts, sizeof(*probe.in_pts));
		probe.in_pts = in_pts;
		probe.out_elev = (double*)safe_calloc(probe.num_pts, sizeof(*probe.out_elev));
		probe.out_water = (double*)safe_calloc(probe.num_pts, sizeof(*probe.out_water));
		probe.filter_lin = B_TRUE;

		for (unsigned i = 0; i < probe.num_pts; i++) {
			vect2_t p = vect2_scmul(v, i / (double)(probe.num_pts - 1));
			in_pts[i] = fpp2geo(p, &fpp);
		}
		egpws_intf->terr_probe(&probe);

		double mask_angle = 0;
		for (unsigned i = 0; i < probe.num_pts; i++) {
			double elevation = probe.out_elev[i] - obs.elev;
			double distance = gc_distance(rec, probe.in_pts[i]);
			if (elevation > 0 && distance > 1) {
				double grade = atan(elevation / distance);

				if (grade >= mask_angle) {
					mask_angle = grade;
				}
			}
		}

		sat->terr_mask_angle = mask_angle;

		free(in_pts);
		free(probe.out_elev);
		free(probe.out_water);
	} else {
		sat->terr_mask_angle = 0;
	}
}

/* Private */
void SpaceLoop::DownloadAlmanac() {
	Log("Downloading GPS Almanac.");

	CURL* curl;
	//FILE* fp;
	std::stringstream out;
	CURLcode res;
	char* out_file = mkpathname(xp_dir, "Output", "caches", "COWS NXi", "current_almanac.al3", NULL);
	char* pem_file = mkpathname(plug_dir, "data", "ca-bundle.crt", NULL);

	curl = curl_easy_init();
	if (curl) {
		curl_easy_setopt(curl, CURLOPT_URL, ALMANAC_URL);
		curl_easy_setopt(curl, CURLOPT_CAINFO, pem_file);
		curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, WriteData);
		curl_easy_setopt(curl, CURLOPT_ACCEPT_ENCODING, "");
		curl_easy_setopt(curl, CURLOPT_FOLLOWLOCATION, 1L);
		curl_easy_setopt(curl, CURLOPT_WRITEDATA, &out);
		res = curl_easy_perform(curl);
		curl_easy_cleanup(curl);

		if (res != CURLE_OK) {
			std::stringstream ss;
			if (res == CURLE_HTTP_RETURNED_ERROR) {
				ss << "An Error has occured with cURL: HTTP Error >= 400";
			} else {
				ss << "An Error has occured with cURL: " << res << " " << curl_easy_strerror(res);
			}
			std::string str = ss.str();
			Log(str.c_str());
		} else {
			std::string str = out.str();

			if (!str.empty() && isdigit(str[0])) {
				std::ofstream almanac(out_file);
				almanac << str;
				almanac.close();
			} else {
				Log("The contents of the URL are not the GPS Almanac.");
			}
		}
	} else {
		Log("cURL failed to init.");
	}

	lacf_free(out_file);
	lacf_free(pem_file);
}

void SpaceLoop::FindAlmanac() {
	char* path;
	if (FileExists(mkpathname(xp_dir, "Output", "caches", "COWS NXi", "current_almanac.al3", NULL))) {
		path = mkpathname(xp_dir, "Output", "caches", "COWS NXi", "current_almanac.al3", NULL);
	} else {
		path = mkpathname(plug_dir, "data", "default_almanac.al3", NULL);
	}
	ReadAlmanac(path);
	lacf_free(path);
}

bool SpaceLoop::FileIsEmpty(std::ifstream& pFile) {
	return pFile.peek() == std::ifstream::traits_type::eof();
}
void SpaceLoop::ReadAlmanac(const char* file_path) {
	Log("Reading GPS Almanac.");

	std::ifstream cFile(file_path);

	if (FileIsEmpty(cFile)) {
		Log("Issue with almanac, reverting to default file.");
		cFile.close();
		file_path = mkpathname(plug_dir, "data", "default_almanac.al3", NULL);
		cFile.open(file_path);
	}

	if (cFile.is_open()) {
		std::string line;
		int line_number = 0;
		GNSS = {};
		Satellite temp_sat = {};
		std::vector<std::string>satellite_lines;

		while (!cFile.eof()) {
			getline(cFile, line);
			line_number++;
			if (line_number <= 2) {
				if (line_number == 1) {
					// num of sats
					std::vector<std::string> result = splitString(line, ' ');
					GNSS.sat_num = stoi(result[0]);
					GNSS.satellites.clear();
					GNSS.satellites.reserve(GNSS.sat_num);
				} else if (line_number == 2) {
					// week, seconds
					std::vector<std::string> result = splitString(line, ' ');
					GNSS.week_num = stoi(result[0]);
					GNSS.applic_secs = stoi(result[1]);
				}
			} else {
				std::string no_white = line;
				no_white.erase(remove_if(no_white.begin(), no_white.end(), isspace), no_white.end());
				if (!no_white.empty()) {
					satellite_lines.push_back(line);
				} else {
					if (satellite_lines.size() == 8) {
						for (size_t i = 0; i < satellite_lines.size(); i++) {
							std::vector<std::string> result = splitString(satellite_lines[i], ' ');
							if (i == 0) {
								temp_sat.prn = stoi(result[0]);
							} else if (i == 1) {
								temp_sat.svn = stoi(result[0]);
							} else if (i == 2) {
								temp_sat.ura = stoi(result[0]);
							} else if (i == 3) {
								std::stringstream ecc(result[0]);
								ecc >> temp_sat.ecc;
								std::stringstream incl(result[1]);
								incl >> temp_sat.incl;
								temp_sat.incl = temp_sat.incl * 3.1415926536;
								std::stringstream ascen(result[2]);
								ascen >> temp_sat.ascen;
								temp_sat.ascen = temp_sat.ascen * 3.1415926536;
							} else if (i == 4) {
								std::stringstream semi_maj(result[0]);
								semi_maj >> temp_sat.semi_maj;
								std::stringstream lng(result[1]);
								lng >> temp_sat.lng;
								temp_sat.lng = temp_sat.lng * 3.1415926536;
								std::stringstream perigee(result[2]);
								perigee >> temp_sat.perigee;
								temp_sat.perigee = temp_sat.perigee * 3.1415926536;
							} else if (i == 5) {
								std::stringstream anom(result[0]);
								anom >> temp_sat.anom;
								temp_sat.anom = temp_sat.anom * 3.1415926536;
								std::stringstream zero_clock(result[1]);
								zero_clock >> temp_sat.zero_clock;
								std::stringstream first_clock(result[2]);
								first_clock >> temp_sat.first_clock;
							} else if (i == 6) {
								temp_sat.health = stoi(result[0]);
							} else if (i == 7) {
								temp_sat.config = stoi(result[0]);
							}
						}
					}

					if (temp_sat.prn != 0) {
						GNSS.satellites.push_back(temp_sat);
						temp_sat = {};
						satellite_lines.clear();
					}
				}
			}
		}
	} else {
		Log("Couldn't read almanac.");
	}
	cFile.close();

	std::time_t t = std::time(nullptr);
	std::tm* const pTInfo = std::localtime(&t);
	year = 1900 + pTInfo->tm_year;

	// start of week in gps unix time
	GNSS.week_epoch = (GNSS.week_num * 7 * 86400) + 1554595200;
	GNSS.week_epoch += GNSS.applic_secs;
	GNSS.week_epoch = UnixToGPS(GNSS.week_epoch);

	std::sort(GNSS.satellites.begin(),
		GNSS.satellites.end(),
		[](const Satellite& lhs, const Satellite& rhs)
		{
			return lhs.prn < rhs.prn;
		});

	for (size_t i = 0; i < GNSS.satellites.size(); i++) {
		CalcSatPos(&GNSS.satellites[i], 0);
	}
}

inline time_t SpaceLoop::UnixToGPS(time_t time) {
	ASSERT3U(time, >= , GPS_EPOCH_UNIX);
	return (time - GPS_EPOCH_UNIX + 18);
}
double SpaceLoop::SolveE(float meanAnomaly, float eccentricity) {
	double E[4] = { 0.0 };
	E[0] = meanAnomaly;
	for (int i = 1; i <= 3; i++) {
		E[i] = E[i - 1] + ((meanAnomaly - E[i - 1] + (eccentricity * sin(E[i - 1]))) / (1 - (eccentricity * cos(E[i - 1]))));
	}
	return E[3];
}

std::vector<std::string> SpaceLoop::splitString(std::string str, char splitter) {
	std::vector<std::string> result;
	std::string current = "";
	for (size_t i = 0; i < str.size(); i++) {
		if (str[i] == splitter) {
			if (current != "") {
				result.push_back(current);
				current = "";
			}
			continue;
		}
		current += str[i];
	}
	if (current.size() != 0)
		result.push_back(current);
	return result;
}

void SpaceLoop::Log(const char* msg) {
	std::stringstream msg_stream;
	msg_stream << "[COWS SpaceLoop]: " << msg << "\n";
	std::string string = msg_stream.str();
	XPLMDebugString(string.c_str());
}

.Space.hpp

#ifndef SPACE_HPP
#define SPACE_HPP

#include <vector>
#include <XPLMDataAccess.h>
#include <acfutils/worker.h>
#include <acfutils/geom.h>

#define	ALMANAC_URL "https://navcen.uscg.gov/sites/default/files/gps/almanac/current_sem.al3"
#define GPS_EPOCH_UNIX (time_t)315964800

struct Satellite {
	int prn;           /* GPST number */
	int svn;           /* Space craft vehicle number*/
	int ura;           /* Average URA */

	float ecc;         /* Eccentricity                   - n/a */
	float incl;        /* Inclination offset             - radians */
	float ascen;       /* Rate of right ascension        - radians/sec */

	float semi_maj;    /* Square root of semi-major axis - meters */
	float lng;		   /* Longitude of orbital plane     - radians */
	float perigee;     /* Argument of perigee            - radians */

	float anom;        /* Mean anomaly                   - radians */
	float zero_clock;  /* Zero order clock correction    - seconds */
	float first_clock; /* First order clock correction   - seconds */

	int health;        /* Health of the satellite */
	int config;        /* Satellite configuration */

	/* Position in ECEF */
	vect3_t ecef;

	/* position in degrees decimal */
	geo_pos3_t llh;

	double Az;	            /* Azimuth degrees */
	double El;              /* Elevation degrees */
	double range;           /* Range to satellite m */

	double terr_mask_angle; /* Mask angle produced by Terrain */
};

struct Constellation {
	int sat_num = 0;                  /* Number of satellites in almanac */
	int week_num = 0;                 /* Week number */
	int applic_secs = 0;              /* GPS time of applicability */
	time_t week_epoch = 0;            /* GPS almanac epoch */

	std::vector<Satellite>satellites; /* contains satellites */
};

/* WAAS */
struct WaasStruct {
	Satellite satellites[3] = {
		{
			131, 131, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ -19080139, -37609008, 0},
			{ -116.9, 0, 35794000},
			0, 0, 0
		},
		{
			133, 133, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ -26596947, -32727535, 0},
			{ -129.1, 0, 35794000},
			0, 0, 0
		},
		{
			138, 138, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ -12538496, -40256477, 661693},
			{ -107.3, 0.9, 35791000},
			0, 0, 0
		}
	};
};

/* EGNOS */
struct EgnosStruct {
	Satellite satellites[2] = {
		{
			123, 123, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ 43144562, 26439021, 0},
			{ 31.5, 0, 44223000}
		},
		{
			136, 136, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ 48201017, 4217043, 0},
			{ 5, 0, 42007000}
		}
	};
};

/* MSAS */
struct MsasStruct {
	Satellite satellites[1] = {
		{
			137, 137, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0, 0,
			0, 0,
			{ -25382835, 33684159, 0},
			{ 127, 0, 35799000}
		}
	};
};

class SpaceLoop {
public:
	SpaceLoop(const char* xp, const char* plug);
	~SpaceLoop();

	int year;

	XPLMDataRef zhrs_dref = NULL;
	XPLMDataRef zmins_dref = NULL;
	XPLMDataRef zsecs_dref = NULL;
	XPLMDataRef days_dref = NULL;

	time_t gps_time;     /* GPS time in X-Plane */
	time_t d_time;       /* Delta time since epoch */
	vect3_t ecef_pos;    /* ecef position of reciever */
	mutex_t prop_mute;   /* Mutex used for propagating */

	Constellation GNSS;

	WaasStruct WAAS;
	EgnosStruct EGNOS;
	MsasStruct MSAS;

	inline time_t UnixToGPS(time_t time);

	void CalcSatPos(Satellite* sat, time_t delta_time);
	void SatLookAngle(Satellite* sat, vect3_t obs);
	void SatMaskAngle(Satellite* sat, geo_pos3_t obs);
private:
	const char* xp_dir;
	const char* plug_dir;

	worker_t worker;     /* Background worker thread */

	XPLMFlightLoopID flight_id;

	void DownloadAlmanac();

	void FindAlmanac();

	bool FileIsEmpty(std::ifstream& pFile);
	std::vector<std::string> splitString(std::string str, char splitter);
	void ReadAlmanac(const char* file_path);

	double SolveE(float meanAnomaly, float eccentricity);

	void Log(const char* msg);
};

#endif

README.md

Read me!

This code uses libacfutils, opengpws and relies being used in X-Plane with their SDK.

To use this class you declare it as a pointer to a SpaceLoop:

SpaceLoop* space = nullptr;

Then when you're ready to construct the class simply do:

space = new SpaceLoop(xp_dir, plug_dir);

When you're done with the class call:

delete space