AI Assistance - QO-100
QO-100 Footprints
In my map software I've wanted to correctly display a set of footprints for the QO-100 / Es' Hail 2 satellite but could not find the algorithm anywhere, especially one which corrects for atmospheric refraction which is the reason why conact can be made when the satellite is just below the horizon. Microsoft's CoPilot took ~30 seconds to show me the algorithm and then deliver a nice C++ class!
Below is the code - very readable and fully commented. The only change I made is to play 1440 points (360 * 4).
#pragma once
#include <cmath>
#include <vector>
//
// Lat/Lon pair (in radians)
//
struct LatLon {
double lat; // latitude (radians)
double lon; // longitude (radians)
};
//
// Compute satellite footprint polygon with atmospheric refraction correction.
//
// This uses the standard radio‑propagation "k‑factor" model:
//
// R_eff = k * R_E
//
// where k ≈ 4/3 for normal atmospheric conditions.
//
// IMPORTANT:
// - R_eff is used ONLY inside the geometry to compute the central angle ψ.
// - The final footprint is still drawn on the REAL Earth (R_E).
//
class RefractedSatelliteFootprintPolygon
{
public:
// Physical Earth radius (mean WGS‑84)
static constexpr double EarthRadius = 6371008.8;
//
// Compute effective Earth radius for refraction.
// k_factor = 1.0 → no refraction
// k_factor = 4/3 → standard radio horizon model
//
static double EffectiveEarthRadius(double k_factor)
{
return k_factor * EarthRadius;
}
//
// Compute central angle ψ (radians) from satellite to footprint boundary.
//
// This is the key geometry:
//
// ψ = acos( (R_eff / (R_eff + h)) * cos(e) ) - e
//
// where:
// R_eff = effective Earth radius (refraction corrected)
// h = satellite altitude above surface (meters)
// e = minimum elevation angle (radians)
//
// ψ is the angular distance from the sub‑satellite point to the footprint edge.
//
static double CentralAngle(double altitude_m,
double elevation_deg,
double k_factor = 4.0 / 3.0)
{
// Convert elevation angle from degrees → radians
const double e = elevation_deg * M_PI / 180.0;
// Effective Earth radius (refraction corrected)
const double Re = EffectiveEarthRadius(k_factor);
// Satellite distance from Earth's center
const double r = Re + altitude_m;
// Compute geometric term inside acos()
const double term = (Re / r) * std::cos(e);
// Clamp to [-1, 1] to avoid numerical domain errors in acos()
const double x = std::max(-1.0, std::min(1.0, term));
// Full central angle formula
return std::acos(x) - e;
}
//
// Compute the full 360‑point footprint polygon.
//
// For each bearing θ = 0…359 degrees:
//
// - Move ψ radians away from the sub‑satellite point
// - Along bearing θ
// - Using great‑circle forward‑geodesic equations
//
// Output lat/lon are in radians.
//
static std::vector<LatLon> ComputePolygon(double satLat_deg,
double satLon_deg,
double altitude_m,
double elevation_deg,
double k_factor = 4.0 / 3.0)
{
// Convert satellite sub‑point to radians
const double lat0 = satLat_deg * M_PI / 180.0;
const double lon0 = satLon_deg * M_PI / 180.0;
// Compute refraction‑corrected central angle ψ
const double psi = CentralAngle(altitude_m, elevation_deg, k_factor);
// Precompute trig values for speed
const double sin_lat0 = std::sin(lat0);
const double cos_lat0 = std::cos(lat0);
const double sin_psi = std::sin(psi);
const double cos_psi = std::cos(psi);
std::vector<LatLon> poly;
poly.reserve(360);
//
// Loop over bearings 0° → 359°
//
for (int bearing_deg = 0; bearing_deg < 360; ++bearing_deg)
{
// Bearing θ in radians
const double theta = bearing_deg * M_PI / 180.0;
//
// Great‑circle forward‑geodesic:
//
// lat = asin( sin(lat0)*cos(ψ) + cos(lat0)*sin(ψ)*cos(θ) )
//
const double sin_lat =
sin_lat0 * cos_psi +
cos_lat0 * sin_psi * std::cos(theta);
const double lat = std::asin(sin_lat);
//
// lon = lon0 + atan2( sin(θ)*sin(ψ)*cos(lat0),
// cos(ψ) - sin(lat0)*sin(lat) )
//
const double y = std::sin(theta) * sin_psi * cos_lat0;
const double x = cos_psi - sin_lat0 * sin_lat;
const double lon = lon0 + std::atan2(y, x);
poly.push_back({ lat, lon });
}
return poly;
}
};
