feedc0de/mp-units
1
0
Fork 0
forked from mpusz/mp-units
Code Pull Requests Activity

feat(example): spherical_distance() now uses quantity-based trigonometric functions

Browse Source
This commit is contained in:
Mateusz Pusz
2023-04-07 22:54:47 +02:00
parent 40f6ae2da7
commit 9666c1a065
1 changed files with 9 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
example/glide_computer/include/geographic.h
View File

@@ -24,6 +24,7 @@
#include "ranged_representation.h" #include "ranged_representation.h"
#include <mp_units/bits/fmt_hacks.h> #include <mp_units/bits/fmt_hacks.h>
#include <mp_units/math.h>
#include <mp_units/quantity.h> #include <mp_units/quantity.h>
#include <mp_units/systems/isq/space_and_time.h> #include <mp_units/systems/isq/space_and_time.h>
#include <mp_units/systems/si/units.h> #include <mp_units/systems/si/units.h>
@@ -139,30 +140,23 @@ template<typename T>
distance spherical_distance(position<T> from, position<T> to) distance spherical_distance(position<T> from, position<T> to)
{ {
using namespace mp_units; using namespace mp_units;
constexpr auto earth_radius = 6371 * isq::radius[si::kilo<si::metre>]; constexpr auto earth_radius = 6'371 * isq::radius[si::kilo<si::metre>];
constexpr auto p = std::numbers::pi_v<T> / 180; using isq::sin, isq::cos, isq::asin, isq::acos;
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;
// https://en.wikipedia.org/wiki/Great-circle_distance#Formulae // https://en.wikipedia.org/wiki/Great-circle_distance#Formulae
if constexpr (sizeof(T) >= 8) { if constexpr (sizeof(T) >= 8) {
// spherical law of cosines // spherical law of cosines
const auto central_angle = const auto central_angle = acos(sin(from.lat) * sin(to.lat) + cos(from.lat) * cos(to.lat) * cos(to.lon - from.lon));
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(to.lat - from.lat) / 2 + cos(from.lat) * cos(to.lat) * (1
// 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 - from.lon)) / 2));
// - cos(lon2_rad - lon1_rad)) / 2));
return quantity_cast<isq::distance>(earth_radius * central_angle); return quantity_cast<isq::distance>(earth_radius * central_angle);
} else { } else {
// the haversine formula // the haversine formula
const auto sin_lat = sin((lat2_rad - lat1_rad) / 2); const auto sin_lat = sin((to.lat - from.lat) / 2);
const auto sin_lon = sin((lon2_rad - lon1_rad) / 2); const auto sin_lon = sin((to.lon - from.lon) / 2);
const auto central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(lat1_rad) * cos(lat2_rad) * sin_lon * sin_lon)); const auto central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(from.lat) * cos(to.lat) * sin_lon * sin_lon));
return quantity_cast<isq::distance>(earth_radius * central_angle); return quantity_cast<isq::distance>(earth_radius * central_angle);
} }