diff --git a/example/glide_computer/include/geographic.h b/example/glide_computer/include/geographic.h index c31f2619..f1efa632 100644 --- a/example/glide_computer/include/geographic.h +++ b/example/glide_computer/include/geographic.h @@ -24,6 +24,7 @@ #include "ranged_representation.h" #include +#include #include #include #include @@ -139,30 +140,23 @@ template distance spherical_distance(position from, position to) { using namespace mp_units; - constexpr auto earth_radius = 6371 * isq::radius[si::kilo]; + constexpr auto earth_radius = 6'371 * isq::radius[si::kilo]; - constexpr auto p = std::numbers::pi_v / 180; - const auto lat1_rad = from.lat.number() * p; - const auto lon1_rad = from.lon.number() * p; - const auto lat2_rad = to.lat.number() * p; - const auto lon2_rad = to.lon.number() * p; - - using std::sin, std::cos, std::asin, std::acos, std::sqrt; + using isq::sin, isq::cos, isq::asin, isq::acos; // https://en.wikipedia.org/wiki/Great-circle_distance#Formulae if constexpr (sizeof(T) >= 8) { // spherical law of cosines - const auto central_angle = - acos(sin(lat1_rad) * sin(lat2_rad) + cos(lat1_rad) * cos(lat2_rad) * cos(lon2_rad - lon1_rad)); - // const auto central_angle = 2 * asin(sqrt(0.5 - cos(lat2_rad - lat1_rad) / 2 + cos(lat1_rad) * cos(lat2_rad) * (1 - // - cos(lon2_rad - lon1_rad)) / 2)); + const auto central_angle = acos(sin(from.lat) * sin(to.lat) + cos(from.lat) * cos(to.lat) * cos(to.lon - from.lon)); + // const auto central_angle = 2 * asin(sqrt(0.5 - cos(to.lat - from.lat) / 2 + cos(from.lat) * cos(to.lat) * (1 + // - cos(lon2_rad - from.lon)) / 2)); return quantity_cast(earth_radius * central_angle); } else { // the haversine formula - const auto sin_lat = sin((lat2_rad - lat1_rad) / 2); - const auto sin_lon = sin((lon2_rad - lon1_rad) / 2); - const auto central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(lat1_rad) * cos(lat2_rad) * sin_lon * sin_lon)); + const auto sin_lat = sin((to.lat - from.lat) / 2); + const auto sin_lon = sin((to.lon - from.lon) / 2); + const auto central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(from.lat) * cos(to.lat) * sin_lon * sin_lon)); return quantity_cast(earth_radius * central_angle); }