diff --git a/README.md b/README.md index 3285d8a..00d4b35 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,9 @@ Latitude/longitude to and from UTM coordinates precise and vectorized conversion ### [GREATCIRCLE and LOXODROME](greatcircle) Shortest and rhumb line path, distance and bearing. +### [VINCENTY AND ELLIPSOID](vincenty) +Calculates the distance (in kilometers) using the formula devised by Thaddeus Vincenty. + ### [RADIOCOVER](radiocover) Radio link coverage map. diff --git a/vincenty/README.md b/vincenty/README.md new file mode 100644 index 0000000..22e56d3 --- /dev/null +++ b/vincenty/README.md @@ -0,0 +1,71 @@ +# Vincenty + +## vincenty.m +Calculates the distance (in kilometers) between **point 1** and **point 2** using the +formula devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the +earth. The default ellipsoidal model is **WGS-84**, which is the most globally accurate model. + +### Forms + +```matlab +dist = vincenty(pt1, pt2) +dist = vincenty(pt1, pt2, ellipsoid) +[dist, az] = vincenty(pt1, pt2) +[dist, az] = vincenty(pt1, pt2, ellipsoid) +``` + + `pt1` and `pt2` are two-column matrices of the form `[latitude longitude]`. + The units for the input coordinates angles must be *degrees*. + `ellipsoid` defines the reference ellipsoid to use. + +Sample values for `ellipsoid` are the following: + +* WGS\_84 (default) - referenceEllipsoid(7030) +* GRS\_80 - referenceEllipsoid(7019) +* Airy - referenceEllipsoid(7001) +* Intl - referenceEllipsoid(7022) +* Clarke - referenceEllipsoid(7012) +* GRS - referenceEllipsoid(7003) + +The sample values are the following: + +| Model | Major (km) | Minor (km) | 1 / f | +| ------------------ | ----------- | ------------- | ------------- | +| WGS 1984 | 6378.137 | 6356.7523142 | 298.257223563 | +| GRS 1980 | 6378.137 | 6356.7523141 | 298.257222101 | +| G.B. Airy 1830 | 6377.563396 | 6356.256909 | 299.3249646 | +| Internacional 1924 | 6378.388 | 6356.911946 | 297.0 | +| Clarke 1880 | 6378.249145 | 6356.51486955 | 293.465 | +| Australian Nat. | 6378.1600 | 6356.774719 | 298.25 | + +### Example + +```matlab +>> vincenty([37, -76], [37, -9]) +ans = 5830.081 +>> vincenty([37, -76], [67, -76], referenceEllipsoid(7019)) +ans = 3337.843 +``` + +## Reference Ellipsoid +This function returns a reference ellipsoid object corresponding to the +specified `code` (numerical EPSG). The values of the `SemimajorAxis` +and `SemiminorAxis` properties are in kilometers. The reference +ellipsoid has five properties: `Code`, `Name`, `SemimajorAxis`, `SemiminorAxis` +and `Flattening`. + +The form code can receive a valid EPSG code. 46 codes are currently +implemented between 7001 and 7053 (except for 7017, 7023, 7026 and +7037-7040). + +The valid values for name form are as follows: sphere, unitsphere, earth, +moon, mercury, venus, mars, jupiter, saturn, uranus, neptune and pluto. + +### Forms + +```matlab +ref = referenceEllipsoid(code) +ref = referenceEllipsoid(name) +``` + +Type `help vincenty` or `help referenceEllipsoid` for syntax, help and examples. diff --git a/vincenty/referenceEllipsoid.m b/vincenty/referenceEllipsoid.m new file mode 100644 index 0000000..0d1378d --- /dev/null +++ b/vincenty/referenceEllipsoid.m @@ -0,0 +1,123 @@ +## Copyright (C) 2014 Alfredo Foltran +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {Function File} {} @var{ref} = referenceEllipsoid(@var{code}) +## @deftypefnx {Function File} {} @var{ref} = referenceEllipsoid(@var{name}) +## +## This function returns a reference ellipsoid object corresponding to the +## specified @var{code} (numerical EPSG). The values of the SemimajorAxis +## and SemiminorAxis properties are in kilometers. The reference +## ellipsoid has five properties: Code, Name, SemimajorAxis, SemiminorAxis +## and Flattening. +## +## The form code can receive a valid EPSG code. 46 codes are currently +## implemented between 7001 and 7053 (except for 7017, 7023, 7026 and +## 7037-7040). +## +## The valid values for name form are as follows: sphere, unitsphere, earth, +## moon, mercury, venus, mars, jupiter, saturn, uranus, neptune and pluto. +## +## @end deftypefn + +## Author: Alfredo Foltran + +function ref = referenceEllipsoid(model) + ELLIPSOIDS.C7001 = {'Airy 1830' 6377563.396 299.3249646}; + ELLIPSOIDS.C7002 = {'Airy Modified 1849' 6377340.189 299.3249646}; + ELLIPSOIDS.C7003 = {'Australian National Spheroid' 6378160 298.25}; + ELLIPSOIDS.C7004 = {'Bessel 1841' 6377397.155 299.1528128}; + ELLIPSOIDS.C7005 = {'Bessel Modified' 6377492.018 299.1528128}; + ELLIPSOIDS.C7006 = {'Bessel Namibia' 6377483.865 299.1528128}; + ELLIPSOIDS.C7007 = {'Clarke 1858' 20926348 294.260676369}; + ELLIPSOIDS.C7008 = {'Clarke 1866' 6378206.4 294.978698213898}; + ELLIPSOIDS.C7009 = {'Clarke 1866 Michigan' 20926631.531 294.978697164674}; + ELLIPSOIDS.C7010 = {'Clarke 1880 (Benoit)' 6378300.789 293.466315538981}; + ELLIPSOIDS.C7011 = {'Clarke 1880 (IGN)' 6378249.2 293.466021293627}; + ELLIPSOIDS.C7012 = {'Clarke 1880 (RGS)' 6378249.145 293.465}; + ELLIPSOIDS.C7013 = {'Clarke 1880 (Arc)' 6378249.145 293.4663077}; + ELLIPSOIDS.C7014 = {'Clarke 1880 (SGA 1922)' 6378249.2 293.46598}; + ELLIPSOIDS.C7015 = {'Everest 1830 (1937 Adjustment)' 6377276.345 300.8017}; + ELLIPSOIDS.C7016 = {'Everest 1830 (1967 Definition)' 6377298.556 300.8017}; + ELLIPSOIDS.C7018 = {'Everest 1830 Modified' 6377304.063 300.8017}; + ELLIPSOIDS.C7019 = {'GRS 1980' 6378137 298.257222101}; + ELLIPSOIDS.C7020 = {'Helmert 1906' 6378200 298.3}; + ELLIPSOIDS.C7021 = {'Indonesian National Spheroid' 6378160 298.247}; + ELLIPSOIDS.C7022 = {'International 1924' 6378388 297}; + ELLIPSOIDS.C7024 = {'Krassowsky 1940' 6378245 298.3}; + ELLIPSOIDS.C7025 = {'NWL 9D' 6378145 298.25}; + ELLIPSOIDS.C7027 = {'Plessis 1817' 6376523 308.64}; + ELLIPSOIDS.C7028 = {'Struve 1860' 6378298.3 294.73}; + ELLIPSOIDS.C7029 = {'War Office' 6378300 296}; + ELLIPSOIDS.C7030 = {'WGS 84' 6378137 298.257223563}; + ELLIPSOIDS.C7031 = {'GEM 10C' 6378137 298.257223563}; + ELLIPSOIDS.C7032 = {'OSU86F' 6378136.2 298.257223563}; + ELLIPSOIDS.C7033 = {'OSU91A' 6378136.3 298.257223563}; + ELLIPSOIDS.C7034 = {'Clarke 1880' 20926202 293.465}; + ELLIPSOIDS.C7035 = {'Sphere' 6371000 Inf}; + ELLIPSOIDS.C7036 = {'GRS 1967' 6378160 298.247167427}; + ELLIPSOIDS.C7041 = {'Average Terrestrial System 1977' 6378135 298.257}; + ELLIPSOIDS.C7042 = {'Everest (1830 Definition)' 20922931.8 300.8017}; + ELLIPSOIDS.C7043 = {'WGS 72' 6378135 298.26}; + ELLIPSOIDS.C7044 = {'Everest 1830 (1962 Definition)' 6377301.243 300.8017255}; + ELLIPSOIDS.C7045 = {'Everest 1830 (1975 Definition)' 6377299.151 300.8017255}; + ELLIPSOIDS.C7046 = {'Bessel Namibia (GLM)' 6377397.155 299.1528128}; + ELLIPSOIDS.C7047 = {'GRS 1980 Authalic Sphere' 6370997 Inf}; + ELLIPSOIDS.C7048 = {'GRS 1980 Authalic Sphere' 6371007 Inf}; + ELLIPSOIDS.C7049 = {'Xian 1980' 6378140 298.257}; + ELLIPSOIDS.C7050 = {'GRS 1967 (SAD69)' 6378160 298.25}; + ELLIPSOIDS.C7051 = {'Danish 1876' 6377019.27 300}; + ELLIPSOIDS.C7052 = {'Clarke 1866 Authalic Sphere' 6370997 Inf}; + ELLIPSOIDS.C7053 = {'Hough 1960' 6378270 297}; + + ELLIPSOIDS.Nsphere = {'Sphere' 6371000 Inf 7035}; + ELLIPSOIDS.Nunitsphere = {'Unit Sphere' 1 Inf NaN}; + ELLIPSOIDS.Nearth = {'WGS 84' 6378137 298.257223563 7030}; + ELLIPSOIDS.Nmoon = {'Moon' 1738100 827.666666666702 NaN}; + ELLIPSOIDS.Nmercury = {'Mercury' 2439700 Inf NaN}; + ELLIPSOIDS.Nvenus = {'Venus' 6051800 Inf NaN}; + ELLIPSOIDS.Nmars = {'Mars' 3396200 169.81 NaN}; + ELLIPSOIDS.Njupiter = {'Jupiter' 71492000 15.4144027598103 NaN}; + ELLIPSOIDS.Nsaturn = {'Saturn' 60268000 10.2079945799458 NaN}; + ELLIPSOIDS.Nuranus = {'Uranus' 25559000 43.6160409556314 NaN}; + ELLIPSOIDS.Nneptune = {'Neptune' 24764000 58.5437352245863 NaN}; + ELLIPSOIDS.Npluto = {'Pluto' 1195 Inf NaN}; + + if nargin ~= 1 + error("Code or name must be specified!"); + endif + + if isnumeric(model) + try + model_values = eval(['ELLIPSOIDS.C' num2str(model)]); + ref.Code = model; + catch + error("Invalid EPSG code!"); + end_try_catch + else + try + model_values = eval(['ELLIPSOIDS.N' num2str(model)]); + ref.Code = model_values{4}; + catch + error("Invalid ellipsoid name!"); + end_try_catch + endif + + ref.Name = model_values{1}; + ref.SemimajorAxis = model_values{2} / 1000; + ref.SemiminorAxis = ref.SemimajorAxis - ref.SemimajorAxis / model_values{3}; + ref.Flattening = 1 / model_values{3}; +endfunction diff --git a/vincenty/vincenty.m b/vincenty/vincenty.m new file mode 100644 index 0000000..0fc9405 --- /dev/null +++ b/vincenty/vincenty.m @@ -0,0 +1,140 @@ +## Copyright (C) 2014 Alfredo Foltran +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see . + +## -*- texinfo -*- +## @deftypefn {Function File} {} @var{dist} = vincenty(@var{pt1}, @var{pt2}) +## @deftypefnx {Function File} {} @var{dist} = vincenty(@var{pt1}, @var{pt2}, @var{ellipsoid}) +## @deftypefnx {Function File} {[@var{dist}, @var{az}] = } {vincenty(@var{pt1}, @var{pt2})} +## @deftypefnx {Function File} {[@var{dist}, @var{az}] = } {vincenty(@var{pt1}, @var{pt2}, @var{ellipsoid})} +## +## Calculates the distance (in kilometers) between @var{pt1} and @var{pt2} using the +## formula devised by Thaddeus Vincenty, with an accurate ellipsoidal model of the +## earth (@var{ellipsoid}). The default ellipsoidal model is 'WGS-84', which is the most +## globally accurate model. +## +## @var{pt1} and @var{pt2} are two-column matrices of the form [latitude longitude]. +## The units for the input coordinates angles must be degrees. +## @var{ellipsoid} defines the reference ellipsoid to use. +## +## Sample values for @var{ellipsoid} are the following: +## +## @multitable @columnfractions .7 .3 +## @headitem Model @tab @var{ellipsoid} +## @item WGS 1984 (default) @tab referenceEllipsoid(7030) +## @item GRS 1980 @tab referenceEllipsoid(7019) +## @item G.B. Airy 1830 @tab referenceEllipsoid(7001) +## @item Internacional 1924 @tab referenceEllipsoid(7022) +## @item Clarke 1880 @tab referenceEllipsoid(7012) +## @item Australian Nat. @tab referenceEllipsoid(7003) +## @end multitable +## +## The sample model values are the following: +## +## @multitable @columnfractions .35 .20 .20 .25 +## @headitem Model @tab Major (km) @tab Minor (km) @tab 1 / f +## @item WGS 1984 @tab 6378.137 @tab 6356.7523142 @tab 298.257223563 +## @item GRS 1980 @tab 6378.137 @tab 6356.7523141 @tab 298.257222101 +## @item G.B. Airy 1830 @tab 6377.563396 @tab 6356.256909 @tab 299.3249646 +## @item Internacional 1924 @tab 6378.388 @tab 6356.911946 @tab 297.0 +## @item Clarke 1880 @tab 6378.249145 @tab 6356.51486955 @tab 293.465 +## @item Australian Nat. @tab 6378.1600 @tab 6356.774719 @tab 298.25 +## @end multitable +## +## Usage: +## @example +## >> vincenty([37, -76], [37, -9]) +## ans = 5830.081 +## >> vincenty([37, -76], [67, -76], referenceEllipsoid(7019)) +## ans = 3337.843 +## @end example +## +## @seealso{distance,referenceEllipsoid} +## @end deftypefn + +## Author: Alfredo Foltran + +function [dist, az] = vincenty(pt1, pt2, ellipsoid) + + if nargin < 3 + ellipsoid = referenceEllipsoid(7030); + endif + + major = ellipsoid.SemimajorAxis; + minor = ellipsoid.SemiminorAxis; + f = ellipsoid.Flattening; + + iter_limit = 20; + + pt1 = deg2rad(pt1); + pt2 = deg2rad(pt2); + + [lat1 lng1] = deal(pt1(1), pt1(2)); + [lat2 lng2] = deal(pt2(1), pt2(2)); + + delta_lng = lng2 - lng1; + + reduced_lat1 = atan((1 - f) * tan(lat1)); + reduced_lat2 = atan((1 - f) * tan(lat2)); + + [sin_reduced1 cos_reduced1] = deal(sin(reduced_lat1), cos(reduced_lat1)); + [sin_reduced2 cos_reduced2] = deal(sin(reduced_lat2), cos(reduced_lat2)); + + lambda_lng = delta_lng; + lambda_prime = 2 * pi; + + i = 0; + while abs(lambda_lng - lambda_prime) > 10e-12 && i <= iter_limit + i += 1; + [sin_lambda_lng cos_lambda_lng] = deal(sin(lambda_lng), cos(lambda_lng)); + sin_sigma = sqrt((cos_reduced2 * sin_lambda_lng) ** 2 + (cos_reduced1 * sin_reduced2 - sin_reduced1 * cos_reduced2 * cos_lambda_lng) ** 2); + + if sin_sigma == 0 + dist = 0; + return; + endif + + cos_sigma = (sin_reduced1 * sin_reduced2 + cos_reduced1 * cos_reduced2 * cos_lambda_lng); + sigma = atan2(sin_sigma, cos_sigma); + sin_alpha = (cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma); + cos_sq_alpha = 1 - sin_alpha ** 2; + + if cos_sq_alpha != 0 + cos2_sigma_m = cos_sigma - 2 * (sin_reduced1 * sin_reduced2 / cos_sq_alpha); + else + cos2_sigma_m = 0.0; # Equatorial line + endif + + C = f / 16.0 * cos_sq_alpha * (4 + f * (4 - 3 * cos_sq_alpha)); + + lambda_prime = lambda_lng; + lambda_lng = (delta_lng + (1 - C) * f * sin_alpha * (sigma + C * sin_sigma * (cos2_sigma_m + C * cos_sigma * (-1 + 2 * cos2_sigma_m ** 2)))); + endwhile + + if i > iter_limit + error("Inverse Vincenty's formulae failed to converge!"); + endif + + u_sq = cos_sq_alpha * (major ** 2 - minor ** 2) / minor ** 2; + A = 1 + u_sq / 16384.0 * (4096 + u_sq * (-768 + u_sq * (320 - 175 * u_sq))); + B = u_sq / 1024.0 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq))); + delta_sigma = (B * sin_sigma * (cos2_sigma_m + B / 4. * (cos_sigma * (-1 + 2 * cos2_sigma_m ** 2) - B / 6. * cos2_sigma_m * (-3 + 4 * sin_sigma ** 2) * (-3 + 4 * cos2_sigma_m ** 2)))); + dist = minor * A * (sigma - delta_sigma); + + if nargout() > 1 + alpha1 = atan2(cos_reduced2 * sin_lambda_lng, cos_reduced1 * sin_reduced2 - sin_reduced1 * cos_reduced2 * cos_lambda_lng); + alpha2 = atan2(cos_reduced1 * sin_lambda_lng, -sin_reduced1 * cos_reduced2 + cos_reduced1 * sin_reduced2 * cos_lambda_lng); + az = [rad2deg(alpha1) rad2deg(alpha2)]; + endif +endfunction