Merge tag 'upstream/1.4.2' into debian/latest

Upstream version 1.4.2

1 files changed, 51 insertions, 17 deletions

diff --git a/inst/meridianarc.m b/inst/meridianarc.m
index 85039e6..235e441 100644
--- a/inst/meridianarc.m
+++ b/inst/meridianarc.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2019-2020 Philip Nienhuis
+## Copyright (C) 2019-2022 Philip Nienhuis
 ##
 ## 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
@@ -15,20 +15,26 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {@var{s} =}  meridianarc (@var{phi}, @var{phi_2}, @var{spheroid}, @var{angleUnit}) 
-## Returns the meridian arch length given two latitudes @var{phi} and @var{phi_2}.
+## @deftypefn  {Function File} {@var{s} =}  meridianarc (@var{phi}, @var{phi_2})
+## @deftypefnx {Function File} {@var{s} =}  meridianarc (@var{phi}, @var{phi_2}, @var{spheroid})
+## @deftypefnx {Function File} {@var{s} =}  meridianarc (@dots{}, @var{angleUnit})
+## Returns the meridian arc length given two latitudes @var{phi} and @var{phi_2}.
 ##
-## If phi_2 is larger than phi the value will be negative.
+## @var{phi} and @var{phi_2} can be scalars, vectors or arrays of any desired
+## size and dimension; but if any is non-scalar, the other must be scalar or
+## have the exact same dimensions.  For any @var{phi_2} larger than @var{phi}
+## the output value will be negative.
 ##
 ## If no spheroid is given the default is wgs84.
-## The angleunit can be degrees or radians (the latter is default).
+##
+## @var{angleUnit} can be 'degrees' or 'radians' (the latter is default).
 ##
 ## Examples
 ## Full options:
 ## @example
 ## s = meridianarc (0, 56, "int24", "degrees")
 ## => s =
-## 6.2087e+06 
+## 6.2087e+06
 ## @end example
 ## Short version:
 ## @example
@@ -42,13 +48,13 @@
 ## => s =
 ## 6208.7
 ## @end example
-## @seealso{referenceEllipsoid} 
+## @seealso{referenceEllipsoid}
 ## @end deftypefn
 
 ## Function supplied by anonymous contributor, see:
 ## https://savannah.gnu.org/patch/index.php?9720
 
-function s = meridianarc (phi, phi_2, spheroid="wgs84", angleUnit="radians")
+function s = meridianarc (phi, phi_2, spheroid="", angleUnit="radians")
 
   persistent intv = "-pi/2, pi/2";
   persistent degintv = "-90, 90";
@@ -57,21 +63,43 @@ function s = meridianarc (phi, phi_2, spheroid="wgs84", angleUnit="radians")
     print_usage ();
   endif
 
-  if (strncmpi (lower (angleUnit), "d", 1) == 1)
+  if (strncmpi (lower (angleUnit), "degrees", numel (angleUnit)) == 1)
     phi = deg2rad (phi);
     phi_2 = deg2rad (phi_2);
     intv = degintv;
   endif
+  ## Check on identical input sizes or one scalar
+  if (! isscalar (phi))
+    if (isscalar (phi_2))
+      phi_2 = phi_2 * ones (size (phi));
+    elseif (! numel (size (phi) == numel (size (phi_2))) || ...
+        ! all (size(phi) == size (phi_2)))
+      error ("meridianarc: both inputs must have same size and dimensions");
+    endif
+  elseif (! isscalar (phi_2))
+    if (isscalar (phi))
+      phi = phi * ones (size (phi_2));
+    elseif (! numel (size (phi) == numel (size (phi_2))) || ...
+        ! all (size(phi) == size (phi_2)))
+      error ("meridianarc: both inputs must have same size and dimensions");
+    endif
+  endif
   if (abs (phi) > pi / 2 || abs (phi_2) > pi / 2)
     error ("meridianarc: latitudes must lie in interval [%s]", intv);
   endif
 
-  if isempty (spheroid)
-    spheroid = "wgs84";
-  end
+  if (isempty (spheroid))
+    E = referenceEllipsoid ("wgs84");
+  else
+    if (isnumeric (spheroid))
+      spheroid = num2str (spheroid);
+    endif
+    E = sph_chk (spheroid);
+  endif
+
 
   ## From: Algorithms for global positioning. Kai Borre and Gilbert Strang pg 373
-  ## Note: Using integral instead of Taylor Expansion 
+  ## Note: Using integral instead of Taylor Expansion
   if (isstruct (spheroid))
     E = spheroid;
   elseif (ischar (spheroid))
@@ -82,16 +110,22 @@ function s = meridianarc (phi, phi_2, spheroid="wgs84", angleUnit="radians")
 
   e_sq = E.Eccentricity ^ 2;
   F = @(x) ((1 - e_sq * sin(x) ^ 2) ^ (-3 / 2));
-  s = E.SemimajorAxis * (1 - e_sq) * quad ( F, phi, phi_2, 1.0e-12);
+  s = zeros (size (phi));
+  for ii=1:numel (phi)
+    s(ii) = E.SemimajorAxis * (1 - e_sq) * quad ( F, phi(ii), phi_2(ii), 1.0e-12);
+  end %for
 
 endfunction
 
 %!test
 %! s = meridianarc (0, 56, "int24", "degrees");
-%! assert (s, 6208700.08662672, 10e-6)
+%! assert (s, 6208700.08662672, 1e-6)
+
+%!test
+%! s = meridianarc ([-1/120; 45-1/120; 89+119/120], [1/120; 45+1/120; 90], ...
+%!                  "int24", "degrees");
+%! assert (s, [1842.92463205; 1852.25585828; 1861.66609497/2], 1e-6);
 
-%!error <spheroid> meridianarc ( 0, pi/4, 7)
-%!error <spheroid> meridianarc ( 0, pi/4, 7i)
 %!error <latitudes> meridianarc (-2, 2)
 %!error <latitudes> meridianarc (-91, 91, "", "d")