IPGP · alfoltran · Oct 22, 2018 · Oct 23, 2018
diff --git 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
@@ -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
@@ -0,0 +1,123 @@
+## Copyright (C) 2014 Alfredo Foltran <[email protected]>
+## 
+## 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
+## <http://www.gnu.org/licenses/>.
+
+## -*- 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 <[email protected]>
+
+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
@@ -0,0 +1,140 @@
+## Copyright (C) 2014 Alfredo Foltran <[email protected]>
+##
+## 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 <http://www.gnu.org/licenses/>.
+
+## -*- 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 <[email protected]>
+
+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