diff --git a/math/conversion.h b/math/conversion.h new file mode 100644 index 0000000..fc992ab --- /dev/null +++ b/math/conversion.h @@ -0,0 +1,54 @@ +#include +#include + +std::tuple gpsToUTM(double lat, double lon) { + // Constants for UTM conversion + constexpr double a = 6378137.0; // WGS-84 major axis + constexpr double f = 1 / 298.257223563; // WGS-84 flattening + constexpr double k0 = 0.9996; // UTM scale factor + constexpr double e = sqrt(f * (2 - f)); // Eccentricity + constexpr double e2 = e * e; + + // UTM Zone + int zone = static_cast((lon + 180) / 6) + 1; + + // Convert latitude and longitude to radians + double lat_rad = lat * M_PI / 180.0; + double lon_rad = lon * M_PI / 180.0; + + // Central meridian of the UTM zone + double lon_origin = (zone - 1) * 6 - 180 + 3; + double lon_origin_rad = lon_origin * M_PI / 180.0; + + // Calculations for UTM coordinates + double N = a / sqrt(1 - e2 * sin(lat_rad) * sin(lat_rad)); + double T = tan(lat_rad) * tan(lat_rad); + double C = e2 / (1 - e2) * cos(lat_rad) * cos(lat_rad); + double A = cos(lat_rad) * (lon_rad - lon_origin_rad); + + double M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * lat_rad + - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * sin(2 * lat_rad) + + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * sin(4 * lat_rad) + - (35 * e2 * e2 * e2 / 3072) * sin(6 * lat_rad)); + + double easting = k0 * N * (A + (1 - T + C) * A * A * A / 6 + + (5 - 18 * T + T * T + 72 * C - 58 * e2) * A * A * A * A * A / 120) + 500000.0; + + double northing = k0 * (M + N * tan(lat_rad) * (A * A / 2 + + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + + (61 - 58 * T + T * T + 600 * C - 330 * e2) * A * A * A * A * A * A / 720)); + + // Correct for southern hemisphere + if (lat < 0) { + northing += 10000000.0; + } + + return std::make_tuple(easting, northing, zone); +} + +std::tuple gpsToGlobalCoord(double lat0, double lon0, double lat1, double lon1) { + // Convert latitude and longitude to utm coordinates + const auto& [e0, n0, zone0] = gpsToUTM(lat0, lon0); + const auto& [e1, n1, zone1] = gpsToUTM(lat1, lon1); + return {e1 - e0, n1 - n0}; +}