diff options
Diffstat (limited to 'inst')
132 files changed, 3035 insertions, 499 deletions
diff --git a/inst/aer2ecef.m b/inst/aer2ecef.m index 8f4deb6..f883374 100644 --- a/inst/aer2ecef.m +++ b/inst/aer2ecef.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020, Felipe Geremia Nievinski -## Copyright (C) 2019-2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2019-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -84,7 +84,7 @@ ## Note: aer2ecef is a mere wrapper for functions geodetic2ecef, aer2enu and ## enu2uvw. ## -## @seealso {ecef2aer, aer2enu, aer2geodetic, aer2ned, referenceEllipsoid} +## @seealso{ecef2aer, aer2enu, aer2geodetic, aer2ned, referenceEllipsoid} ## @end deftypefn ## Function adapted from patch by anonymous contributor, see: @@ -115,21 +115,19 @@ function [x,y,z] = aer2ecef (az, el, slantrange, lat0, lon0, alt0, ... if (! all (size (lat0) == size (az)) || ... ! all (size (lon0) == size (el)) || ... ! all (size (alt0) == size (slantrange))) - error ("aer2ecef: non-matching dimensions of observer points and \ -target points"); + error (["aer2ecef: non-matching dimensions of observer points and ", ... + "target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid) ||isnumeric (spheroid)) - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid) && isscalar (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); + %% Origin of the local system in geocentric coordinates. - [x0, y0, z0] = geodetic2ecef (spheroid, lat0, lon0, alt0, angleUnit); + [x0, y0, z0] = geodetic2ecef (E, lat0, lon0, alt0, angleUnit); %% Convert Local Spherical AER to ENU [e, n, u] = aer2enu (az, el, slantrange, angleUnit); %% Rotating ENU to ECEF diff --git a/inst/aer2enu.m b/inst/aer2enu.m index d347266..8c854b2 100644 --- a/inst/aer2enu.m +++ b/inst/aer2enu.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2019-2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2019-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without modification, are permitted ## provided that the following conditions are met: @@ -64,7 +64,7 @@ ## u = 866.03 ## @end example ## -## @seealso {enu2aer, aer2ecef, aer2geodetic, aer2ned} +## @seealso{enu2aer, aer2ecef, aer2geodetic, aer2ned} ## ## @end deftypefn diff --git a/inst/aer2geodetic.m b/inst/aer2geodetic.m index 36dba0c..fbdce05 100644 --- a/inst/aer2geodetic.m +++ b/inst/aer2geodetic.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020, Felipe Geremia Nievinski -## Copyright (C) 2019-2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2019-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -80,7 +80,7 @@ ## Note: aer2geodetic is a mere wrapper for functions aer2ecef followed by ## ecef2geodetic. ## -## @seealso {geodetic2aer, aer2ecef, aer2enu, aer2ned, referenceEllipsoid} +## @seealso{geodetic2aer, aer2ecef, aer2enu, aer2ned, referenceEllipsoid} ## @end deftypefn ## Function adapted by anonymous contributor, see: @@ -109,18 +109,14 @@ function [lat1, lon1, alt1] = aer2geodetic (az, el, slantrange, lat0, lon0, alt0 endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid) || isnumeric (spheroid)) - E = referenceEllipsoid (spheroid); - else - error ("aer2geodetic: illegal value for 'spheroid'."); + if (isnumeric (spheroid) && isscalar (spheroid)) + spheroid = num2str (spheroid); endif - [x, y, z] = aer2ecef (az, el, slantrange, lat0, lon0, alt0, spheroid, angleUnit); - [lat1, lon1, alt1] = ecef2geodetic (spheroid, x, y, z, angleUnit); + E = sph_chk (spheroid); + + [x, y, z] = aer2ecef (az, el, slantrange, lat0, lon0, alt0, E, angleUnit); + [lat1, lon1, alt1] = ecef2geodetic (E, x, y, z, angleUnit); endfunction diff --git a/inst/aer2ned.m b/inst/aer2ned.m index 5126968..dc0c01f 100644 --- a/inst/aer2ned.m +++ b/inst/aer2ned.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020, Felipe Geremia Nievinski -## Copyright (C) 2019-2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022, Felipe Geremia Nievinski +## Copyright (C) 2019-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -62,7 +62,7 @@ ## d = -866.03 ## @end example ## -## @seealso {ned2aer, aer2ecef, aer2enu, aer2geodetic} +## @seealso{ned2aer, aer2ecef, aer2enu, aer2geodetic} ## ## @end deftypefn diff --git a/inst/almanac.m b/inst/almanac.m new file mode 100644 index 0000000..61b03f5 --- /dev/null +++ b/inst/almanac.m @@ -0,0 +1,142 @@ +## Copyright (C) 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 +## 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 <https://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {@var{retval} =} almanac () +## @deftypefnx {} {@var{retval} =} almanac (@var{body}) +## @deftypefnx {} {@var{retval} =} almanac (@var{body}, @var{param}) +## @deftypefnx {} {@var{retval} =} almanac (@var{body}, @var{param}, @var{unit}) +## Return basic info about one parameter of a celestial body. +## +## The celestial bodies and parameters that almanac can return are shown when +## invoking almanac without any parameters. +## +## Input parameters: +## +## @itemize +## @item @var{body} - celestial body for which info is requested +## (case-insensitive). +## +## @item @var{param} (optional) - the selected parameter whose value is +## requested. If 'geoid' or 'ellipsoid' is specified a 1x2 vector comprising +## the semimajor axis and eccentricity of @var{body} is returned. If 'Earth' +## is selected, a specific geoid name can be specified and the semimajor +## axis and eccentricity of that geoid are returned. +## +## @item @var{unit} (optional) - the unit in which the requested parameter is +## expressed. Any unit recognized by function 'unitsratio' is accepted. +## @end itemize +## +## almamanc is merely a wrapper around, resp. based on, referenceEllipsoid. +## +## @seealso{referenceEllipsoid,referenceSphere,unitsratio} +## @end deftypefn + +function retval = almanac (varargin) + + persistent a geoids; + outerbodies = {"Sun", "Mercury", "Venus", "Moon", "Mars", "Jupiter", ... + "Saturn", "Neptune", "Uranus", "Pluto", "Unit Sphere"}; + unit = "meter"; + if (isempty (a) || isempty (geoids)) + a = referenceEllipsoid (0); + geoids = setdiff (a(:, 3), [outerbodies {"Earth"}])'; + endif + + if (nargin == 0) + ## Return info on celestial bodies + if (nargout == 0) + printf (["Implemented celestial bodies:\nSun\nMercury\nVenus\n", ... + "Earth\nMoon\nMars\n\Jupiter\nSaturn\nNeptune\nUranus\n", ... + "Pluto\nUnit Sphere\n"]); + else + retval = [outerbodies(1:3) {"Earth"} outerbodies(4:end)]; + endif + + elseif (nargin > 0 && ! iscellstr (varargin)) + error ("almanac: all input args should be text strings."); + + elseif (nargin == 1) + ## Return info on parameters of a specific celestial body + if (strcmpi (varargin{1}, 'earth')) + printf ("Parameters for Earth:\n ", varargin{1}); + printf (["radius, geoid, volume, surfacearea, or one of geoids ", ... + "listed below.\nUnits:\n degrees (deg), kilometers ", ... + "(km), nautical miles (nm), radians (rad) or statute ", ... + "miles (sm)\nReference bodies:\n sphere, geoid or actual\n"]); + idx = find (! ismember (a(:, 3), [outerbodies {"Earth"}])); + printf ("Available geoids:\n %s\n", ... + strjoin (cellfun (@(x) sprintf ("'%s'", x), a(idx, 3), "uni", 0), ", ")); + elseif (ismember (varargin{1}, lower (outerbodies))) + printf ("Parameters for %s:\n ", varargin{1}); + printf (["radius, geoid, volume or surfacearea.\nUnits:\n degrees ", ... + "(deg), kilometers (km), nauticalmiles (nm), radians (rad) ", ... + "or statutemiles (sm)\nReference bodies:\n sphere, geoid ", ... + "or actual\n"]); + else + error ("almanac: unknown celestial body."); + end + if (nargout > 0) + retval = []; + endif + + elseif (nargin >= 2) + ## Return info on a specific celestial body with units + param = varargin{2}; + if (nargin > 2) + unit = varargin{3}; + endif + idx = find (strncmpi (param, geoids, numel (param))); + if (! isempty (idx)) + b = referenceEllipsoid (geoids{idx(1)}, unit); + retval = [b.SemimajorAxis, b.Eccentricity]; + else + b = referenceEllipsoid (varargin{1}, unit); + switch lower (param) + case "radius" + retval = b.MeanRadius; + case "volume" + retval = b.Volume; + case {"surfarea", "surfacearea"} + retval = b.SurfaceArea; + case {"geoid", "ellipsoid"} + retval = [b.SemimajorAxis, b.Eccentricity]; + otherwise + error ("almanac: unknown parameter requested - %s", param); + endswitch + endif + + endif + +endfunction + + +%!test +%! assert (strcmpi (almanac (){1}, "Sun"), true); +%! assert (strcmpi (almanac (){12}, "Unit Sphere"), true); + +%!test +%! assert (almanac ("sun", "radius", "sm"), 432285.77700111, 1e-6); + +%!test +%! assert (almanac ("earth", "everest", "nm"), [3443.456768421318, 0.0814729809], 1e-9); + +%!test +%! assert (almanac ("jupiter", "ellipsoid", "km"), [71492.0 0.3543164], 1e-7); + +%!error <unknown> almanac ("UFO") +%!error <all input args> almanac ("Moon", 12) +%!error <unknown parameter> almanac ("Mars", "flattening") + diff --git a/inst/angl2str.m b/inst/angl2str.m index 9a2e8aa..eaaa344 100644 --- a/inst/angl2str.m +++ b/inst/angl2str.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Ricardo Fantin da Costa +## Copyright (C) 2018-2022 Ricardo Fantin da Costa ## ## 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 diff --git a/inst/angltostr.m b/inst/angltostr.m index e3874e2..4da3dd8 100644 --- a/inst/angltostr.m +++ b/inst/angltostr.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 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 @@ -113,7 +113,7 @@ function str = angltostr (ang, hemcode="none", unit="degrees", acc=-2) ## Split up in degrees, minutes and seconds degs = sign (ang) .* floor (abs (ang)); - mins = abs (ang .- degs) * 60; + mins = abs (ang - degs) * 60; secs = (mins - floor (mins)) * 60; is = sign (ang); ## Set zero values to positive diff --git a/inst/antipode.m b/inst/antipode.m index 1476693..2959076 100644 --- a/inst/antipode.m +++ b/inst/antipode.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 diff --git a/inst/areaquad.m b/inst/areaquad.m new file mode 100644 index 0000000..e90faae --- /dev/null +++ b/inst/areaquad.m @@ -0,0 +1,262 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{aq} =} areaquad (@var{lat1}, @var{lon1}, @var{lat2}, @var{lon2}) +## @deftypefnx {} {@var{aq} =} areaquad (@var{lat1}, @var{lon1}, @var{lat2}, @var{lon2}, @var{spheroid}) +## @deftypefnx {} {@var{aq} =} areaquad (@var{lat1}, @var{lon1}, @var{lat2}, @var{lon2}, @var{angleUnit}) +## @deftypefnx {} {@var{aq} =} areaquad (@var{lat1}, @var{lon1}, @var{lat2}, @var{lon2}, @var{spheroid}, @var{angleUnit}) +## Returns the area of a quadrilateral given two points. +## +## If no ellipsoid is given the result will be a fraction of a unit sphere +## with a radius of one meter. If an ellipsoid struct is supplied the result's +## unit will be the squared unit of that ellipsoid; otherwise, if just an +## ellipsoid name is given result will be in standard units squared (e.g., m^2). +## +## Input +## @itemize +## @item +## @var{lat1}, @var{lon1}: the first point of the quadrilateral. +## +## @item +## @var{lat2}, @var{lon2}: the second point of the quadrilateral. +## @end itemize +## +## @indentedblock +## These coordinate inputs be scalars, vectors or 2D arrays. If any of these +## is non-scalar, all other LAT/LON inputs must either have the same size or +## be scalars. +## @end indentedblock +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid parameter struct; default +## is 'Unit Sphere'. A string value or EPSG code describing the +## spheroid is also accepted. Alternatively, a numerical vector +## @code{[semimajoraxis eccentricity]} can be supplied. If specified, +## @var{spheroid} must be the 5th input argument. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or +## 'radians', case-insensitive, just the first charachter will do). Default +## is 'degrees'. If specified, @var{angleUnit} must be the last input argument. +## @end itemize +## +## Output: +## @itemize +## @item +## @var{aq}: the area of the quadrilateral. If one or more of the inputs were +## vectors or arrays, areaquad's output will have the same size. +## @end itemize +## +## Example +## No ellipsoid (fraction of a sphere) +## @example +## aq = areaquad (0, 0, 90, 360) +## aq = 0.50000 +## @end example +## +## With radians +## @example +## aq = areaquad (0, 0, pi / 2 , 2 * pi, "radians") +## aq = 0.50000 +## @end example +## +## With ellipsoid in m^2 +## @example +## aq = areaquad (0, 0, 90, 360, "wgs84") +## aq = 2.5503e+14 +## @end example +## +## With vector +## @example +## aq = areaquad(-90,0,90,360,[6000 0]); +## aq = 4.5239e+08 +## @end example +## +## @seealso{referenceEllipsoid,referenceSphere} +## @end deftypefn + +function [aq] = areaquad (varargin) + + spheroid = ""; + angleUnit = "degrees"; + insq = 0; # When a spheriod is given the area is in units squared + + if (nargin < 4 || nargin > 6) + print_usage(); + elseif (nargin == 5) + if (isnumeric (varargin{5})) + ## EPSG spheroid code + spheroid = varargin{5}; + elseif (ischar (varargin{5})) + if (! isempty (varargin{5}) && ismember (lower (varargin{5}(1)), {"r", "d"})) + angleUnit = varargin{5}; + else + spheroid = varargin{5}; + endif + elseif (isstruct (varargin{5})) + spheroid = varargin{5}; + else + error ("areaquad: spheroid or angleUnit expected for arg. #5"); + endif + elseif (nargin == 6) + spheroid = varargin{5}; + angleUnit = varargin{6}; + endif + if (isnumeric (spheroid) && isscalar (spheroid)) + spheroid = num2str (spheroid); + endif + + if (! (all (cellfun ("isnumeric", varargin(1:4))) && ... + all (cellfun ("isreal", varargin(1:4))))) + error ("areaquad: numeric values expected for first four inputs"); + endif + + ism = ! cellfun ("isscalar", varargin(1:4)); + if (! isempty (find (ism))) + ## Some of the location inputs are vector or matrix. Check those sizes + sz = reshape (cell2mat (cellfun (@(x) size (x), varargin(ism), "uni", 0)), 2, [])'; + ## Check if matrix sizes are all identical + if (numel (find (ism)) > 1 && any (diff (sz))) + error ("areaquad: sizes of input location vectors/matrices must match."); + endif + ## Make sure all scalar inputs become matrices of same size + varargin(! ism) = cellfun (@(x) repmat (x, sz(ism(1), 1), sz(ism(1), 2)), ... + varargin(! ism), "uni", 0); + else + nv = 1; + endif + lat1 = varargin{1}; + lon1 = varargin{2}; + lat2 = varargin{3}; + lon2 = varargin{4}; + + if (! ischar (angleUnit)) + error ("areaquad: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs ([lat1 lat2]) > 90)) + error ("areaquad: azimuth value(s) out of acceptable range [-90, 90]") + endif + lat1 = deg2rad (lat1); + lon1 = deg2rad (lon1); + lat2 = deg2rad (lat2); + lon2 = deg2rad (lon2); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] as azimuth isn't defined outside + if (any (abs ([lat1 lat2]) > pi / 2)) + error("areaquad: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("areaquad: illegal input for 'angleUnit'"); + endif + + if (isempty (spheroid)) + E = referenceSphere (); + insq = 1; # Returns the fraction of the surface area + else + E = sph_chk (spheroid); + endif + + a = E.SemimajorAxis; + e = E.Eccentricity; + s1 = sin (lat1); + s2 = sin (lat2); + del = lon1 - lon2; + + if (e < eps) + aq = abs ((del .* a .^ 2) .* (s2 - s1)); + else + ## From Earl Burkholder: 3d Global Spatial Data Model 2nd Ed. pg 168. + ## Make the equation more readable + e2 = e .^ 2; + f = 1 ./ (2 .* e); + e2m1 = 1 - e2; + + s21 = s1 .^ 2; + + s22 = s2 .^ 2; + se1 = 1 - e2 .* s21; + se2 = 1 - e2 .* s22; + + c = (del .* a .^2 .* e2m1) ./ 2; + t1 = 1 + e .* s1; + t2 = 1 + e .* s2; + + b1 = 1 - e .* s1; + b2 = 1 - e .* s2; + + g = f .* (log ( t2 ./ b2) - (log (t1 ./ b1))); + + aq = abs (c .* ((s2 ./ (se2)) - (s1 ./ (se1)) + g)); + endif + + if (insq == 1) + aq = aq ./ E.SurfaceArea; + endif + +endfunction + + +%!test +%! aq = areaquad (-90, 0, 90, 360); +%! assert (aq, 1, 1e-8); + +%!test +%! aq = areaquad (-pi / 2, 0, pi / 2, 2 * pi,"r"); +%! assert (aq, 1, 1e-8); + +%!test +%! aq = areaquad (0, 0, 90, 360, "wgs84"); +%! assert (aq, 2.5503281086204421875e+14, 1e-8); + +%!test +%! aq = areaquad(-90,0,90,360,[0 6000]); +%! assert (aq, 4 * pi * 6000 ^ 2, 1e-8) + +%!test +%! aq = areaquad(-90,0,90,360,[6000 0]); +%! assert (aq, 4 * pi * 6000 ^ 2, 1e-8) + +%!test +%! aq = areaquad(-90,0,90,360,[1 0]); +%! assert (aq, 4 * pi, 1e-8) + +%!test +%! aq = areaquad(-90,0,90,360,[0 1]); +%! assert (aq, 4 * pi, 1e-8) + +%!test +%! assert (areaquad ([-60 -45; -30 0], 0, 45, 180), ... +%! [0.393283046242, 0.353553390593; ... +%! 0.301776695296, 0.176776695296], 1e-11); + +%!error <numeric> areaquad ("s", 0, 90, 360) +%!error <numeric> areaquad (5i, 0, 90, 360) +%!error <numeric> areaquad (0, "s", 90, 360) +%!error <numeric> areaquad (0, 5i, 90, 360) +%!error <numeric> areaquad (0, 0, "s", 360) +%!error <numeric> areaquad (0, 0, 5i, 360) +%!error <numeric> areaquad (0, 0, 90, "s") +%!error <numeric> areaquad (0, 0, 90, 5i) + +%!error <ecc> areaquad (-90, 0, 90, 360, [2.9 6000]) +%!error <ellipsoid> areaquad (-90, 0, 90, 360, [2.9]) +%!error <element> areaquad (-90, 0, 90, 360, [2.9 6000 70]) +%!error <azimuth value> areaquad (-91, 0, 90, 360); +%!error <azimuth value> areaquad (-90, 0, 91, 360); +%!error <sizes of> areaquad ([1 2], 3, [4; 5], 6); + diff --git a/inst/auth2geo.m b/inst/auth2geo.m new file mode 100644 index 0000000..b2d5a33 --- /dev/null +++ b/inst/auth2geo.m @@ -0,0 +1,116 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{phi} =} auth2geo (@var{xi}) +## @deftypefnx {} {@var{phi} =} auth2geo (@var{xi}, @var{spheroid}) +## @deftypefnx {} {@var{phi} =} auth2geo (@var{xi}, @var{spheroid}, @var{angleUnit}) +## Returns the geodetic latitude given authalic latitude @var{xi}. +## +## Input +## @itemize +## @item +## @var{xi}: the authalic latitude(s); scalar, vector or nD array. +## +## @item +## (optional) @var{spheroid}: referenceEllipsoid name, EPSG code, or parameter +## struct: the default is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first charachter will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @var{phi}: the geodetic latitude(s), shape matching that of @var{xi}. +## @item +## @end itemize +## +## Example +## @example +## auth2geo(44.872) +## ans = 45.000 +## @end example +## @seealso{geo2auth, geo2con, geo2iso, geo2rect} +## @end deftypefn + + +function [phi] = auth2geo (xi, spheroid = "", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (xi) || ! isreal (xi)) + error ("auth2geo: numeric input expected"); + endif + + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif + E = sph_chk (spheroid); + + if (! ischar (angleUnit)) + error ("auth2geo: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (xi) > 90)) + error ("auth2geo: azimuth value(s) out of acceptable range [-90, 90]") + endif + xi = deg2rad (xi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (phi) > pi / 2)) + error("auth2geo: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("auth2geo: illegal input for 'angleUnit'"); + endif + + + ecc = E.Eccentricity; + ## From Snyder's "Map Projections-A Working Manual" [pg 16]. + e2 = ecc ^ 2; + e4 = ecc ^ 4; + e6 = ecc ^ 6; + phi = xi + ... + (((e2 / 3) + (31 / 180) * e4 + (517 /5040) * e6) * sin(2 * xi)) + ... + (((23 / 360) * e4 + (251 / 3780) * e6) * sin(4 * xi)) + ... + (((761 / 45360) * e6) * sin(6 * xi)); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + phi = rad2deg (phi); + endif + +endfunction + + +%!test +%! xi = [0:15:90]; +%! phi = auth2geo (xi); +%! Z = degrees2dm(phi - xi); +%! check = [0 +%! 3.8575173 +%! 6.6751024 +%! 7.6978157 +%! 6.6579355 +%! 3.8403504 +%! 0]; +%! assert (Z(:,2), check, 1e-6) + +%!error <numeric> auth2geo ("s") +%!error <numeric> auth2geo (5i) + diff --git a/inst/axes2ecc.m b/inst/axes2ecc.m index aacda22..85e2fe5 100644 --- a/inst/axes2ecc.m +++ b/inst/axes2ecc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 @@ -34,13 +34,13 @@ ## axes2ecc (6378137, 6356752.314245) ## => 0.0818191908429654 ## @end example -## +## ## Row vector (semimajor, semiminor): ## @example ## axes2ecc ([6378137, 6356752.314245]) ## => 0.0818191908429654 ## @end example -## +## ## Multivectors: ## @example ## axes2ecc ([ 71492, 66854; ... @@ -49,8 +49,8 @@ ## 0.3543163789650412 ## 0.0818191908429654 ## @end example -## -## @seealso(ecc2flat,flat2ecc) +## +## @seealso{ecc2flat,flat2ecc} ## @end deftypefn ## Function supplied by anonymous contributor, see: @@ -69,10 +69,10 @@ function ecc = axes2ecc (semimajor, semiminor=[]) if (s(2) != 2) error ("axes2ecc: Nx2 matrix expected for arg. #1"); endif - ecc = sqrt ((semimajor(:, 1) .^ 2 .- semimajor(:, 2) .^ 2) ./ ... + ecc = sqrt ((semimajor(:, 1) .^ 2 - semimajor(:, 2) .^ 2) ./ ... (semimajor(:, 1) .^ 2)); else - ecc = sqrt ((semimajor .^ 2 .- semiminor .^ 2) ./ (semimajor .^ 2)); + ecc = sqrt ((semimajor .^ 2 - semiminor .^ 2) ./ (semimajor .^ 2)); endif endfunction diff --git a/inst/azimuth.m b/inst/azimuth.m index 698e0fd..a898738 100644 --- a/inst/azimuth.m +++ b/inst/azimuth.m @@ -1,5 +1,5 @@ -## Copyright (C) 2004-2020 Andrew Collier <abcollier@users.sourceforge.net> -## Copyright (C) 2006-2020 Alexander Barth <abarth93@users.sourceforge.net> +## Copyright (C) 2004-2022 Andrew Collier <abcollier@users.sourceforge.net> +## Copyright (C) 2006-2022 Alexander Barth <abarth93@users.sourceforge.net> ## ## 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 diff --git a/inst/closePolygonParts.m b/inst/closePolygonParts.m index 9a60863..fe8f789 100644 --- a/inst/closePolygonParts.m +++ b/inst/closePolygonParts.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 diff --git a/inst/con2geo.m b/inst/con2geo.m new file mode 100644 index 0000000..36aea90 --- /dev/null +++ b/inst/con2geo.m @@ -0,0 +1,114 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{phi} =} con2geo (@var{chi}) +## @deftypefnx {} {@var{phi} =} con2geo (@var{chi}, @var{spheroid}) +## @deftypefnx {} {@var{phi} =} con2geo (@var{chi}, @var{spheroid}, @var{angleUnit}) +## Returns the geodetic latitude given latitude conformal @var{chi} +## +## Input +## @itemize +## @item +## @var{chi}: the conformal latitude(s). Scalar, vector or nD arrays are +## accepted. +## +## @item +## (optional) @var{spheroid}: referenceEllipsoid parameter struct: the default +## is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first charachter will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{phi}: the geodetic latitude(s), shape similar to @var{chi}. +## @end itemize +## +## Example +## @example +## con2geo(44.8077) +## ans = 45.000 +## @end example +## @seealso{geo2auth, geo2con, geo2iso, geo2rect, rect2geo} +## @end deftypefn + +function [phi] = con2geo (chi, spheroid ="", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (chi) || ! isreal (chi)) + error ("con2geo: numeric input expected"); + endif + + E = sph_chk (spheroid); + + if (! ischar (angleUnit)) + error ("con2geo: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (chi) > 90)) + error ("con2geo: azimuth value(s) out of acceptable range [-90, 90]") + endif + chi = deg2rad (chi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (chi) > pi / 2)) + error("con2geo: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("con2geo: illegal input for 'angleUnit'"); + endif + + ecc = E.Eccentricity; + ## From Snyder's "Map Projections-A Working Manual" [pg 15]. + e2 = ecc ^ 2; + e4 = ecc ^ 4; + e6 = ecc ^ 6; + e8 = ecc ^ 8; + phi = chi + ... + ((e2 / 2 + (5 /24) * e4 + e6 / 12 + (13 / 360) * e8) * sin (2 * chi)) + ... + ((7/48 * e4 + (29 / 240) * e6 + (811 / 11520) * e8) * sin (4 * chi)) + ... + (((7 / 120) * e6 + (81 / 1120) * e8) * sin (6 * chi)) + ... + (((4279 / 161280) * e8) * sin (8 * chi)); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + phi = rad2deg (phi); + endif + +endfunction + + +%!test +%! chi = [0:15:90]; +%! phi = con2geo (chi); +%! Z = degrees2dm (phi - chi); +%! check = [0 +%! 5.7891134 +%! 10.01261 +%! 11.538913 +%! 9.9734791 +%! 5.7499818 +%! 0]; +%! assert (Z(:,2), check, 1e-6) + +%!error <numeric> con2geo ("s") +%!error <numeric> con2geo (5i) + diff --git a/inst/data/egm96_plot.m b/inst/data/egm96_plot.m index 2d92d23..b40f04e 100644 --- a/inst/data/egm96_plot.m +++ b/inst/data/egm96_plot.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2020-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 diff --git a/inst/deg2km.m b/inst/deg2km.m index d50d6f5..bc71a98 100644 --- a/inst/deg2km.m +++ b/inst/deg2km.m @@ -1,4 +1,4 @@ -## Copyright (C) 2013-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2013-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/deg2nm.m b/inst/deg2nm.m index 728d266..56fb94f 100644 --- a/inst/deg2nm.m +++ b/inst/deg2nm.m @@ -1,5 +1,5 @@ -## Copyright (C) 2013-2020 Alexander Barth -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2013-2022 Alexander Barth +## Copyright (C) 2018-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 the Free Software diff --git a/inst/deg2sm.m b/inst/deg2sm.m index ea60b2a..a95565e 100644 --- a/inst/deg2sm.m +++ b/inst/deg2sm.m @@ -1,5 +1,5 @@ -## Copyright (C) 2013-2020 Alexander Barth -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2013-2022 Alexander Barth +## Copyright (C) 2018-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 the Free Software diff --git a/inst/degrees2dm.m b/inst/degrees2dm.m index 1bfb713..a990c42 100644 --- a/inst/degrees2dm.m +++ b/inst/degrees2dm.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/degrees2dms.m b/inst/degrees2dms.m index 3e9d9f3..31f4fed 100644 --- a/inst/degrees2dms.m +++ b/inst/degrees2dms.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/degtorad.m b/inst/degtorad.m index 7ac81f8..614e8d9 100644 --- a/inst/degtorad.m +++ b/inst/degtorad.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/departure.m b/inst/departure.m index 4da98d5..aa630c7 100644 --- a/inst/departure.m +++ b/inst/departure.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 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 @@ -107,7 +107,9 @@ function [dist] = departure (varargin) spheroid = varargin{4}; angleUnit = varargin{5}; endif - + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif long1 = varargin{1}; long2 = varargin{2}; @@ -133,7 +135,7 @@ function [dist] = departure (varargin) lat = rad2deg (lat); endif - delta = abs (long2 .- long1); + delta = abs (long2 - long1); if (nargin == 3 || angck) dist = delta .* cosd (lat); @@ -141,13 +143,7 @@ function [dist] = departure (varargin) dist = deg2rad (dist); endif else - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid)) - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); ## To include spheroid use ## https://en.wikipedia.org/wiki/Longitude#Length_of_a_degree_of_longitude diff --git a/inst/distance.m b/inst/distance.m index 674c088..41523e6 100644 --- a/inst/distance.m +++ b/inst/distance.m @@ -1,6 +1,6 @@ -## Copyright (C) 2004-2020 Andrew Collier <abcollier@users.sourceforge.net> -## Copyright (C) 2011-2020 Alexander Barth <abarth93@users.sourceforge.net> -## Copyright (C) 2019-2020 Philip Nienhuis <prnienhuis@users.sf.net> +## Copyright (C) 2004-2022 Andrew Collier <abcollier@users.sourceforge.net> +## Copyright (C) 2011-2022 Alexander Barth <abarth93@users.sourceforge.net> +## Copyright (C) 2019-2022 Philip Nienhuis <prnienhuis@users.sf.net> ## ## 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 @@ -40,7 +40,7 @@ ## ans = 3.1416 ## @end example ## -## @seealso{azimuth,elevation} +## @seealso{azimuth, elevation, geodeticarc} ## @end deftypefn ## Author: Andrew Collier <abcollier@users.sourceforge.net> @@ -89,12 +89,12 @@ function [dist, az] = distance (varargin) C = deg2rad (C); endif - id = abs (a .- b) < 0.4 & abs (C) < 0.4; + id = abs (a - b) < 0.4 & abs (C) < 0.4; ## Plain spherical cosine formula - dist(! id) = acos (sin (b(! id)) .* sin (a(! id)) .+ ... + dist(! id) = acos (sin (b(! id)) .* sin (a(! id)) + ... cos (b(! id)) .* cos (a(! id)) .* cos (C(! id))); ## Haversine formula for small distances - dist(id) = 2 * asin (sqrt ((sin (b(id) .- a(id)) / 2) .^ 2 .+ ... + dist(id) = 2 * asin (sqrt ((sin (b(id) - a(id)) / 2) .^ 2 + ... cos (b(id)) .* cos (a(id)) .* sin (C(id) / 2) .^ 2)); if (strncmpi (units, "degrees", luni)) diff --git a/inst/dm2degrees.m b/inst/dm2degrees.m index 7d560d6..58f3a18 100644 --- a/inst/dm2degrees.m +++ b/inst/dm2degrees.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/dms2degrees.m b/inst/dms2degrees.m index 69bc5cd..594b2d7 100644 --- a/inst/dms2degrees.m +++ b/inst/dms2degrees.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/dxfdraw.m b/inst/dxfdraw.m index d70da84..898157c 100644 --- a/inst/dxfdraw.m +++ b/inst/dxfdraw.m @@ -1,4 +1,4 @@ -## Copyright (C) 2016-2020 Philip Nienhuis +## Copyright (C) 2016-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 @@ -83,8 +83,8 @@ function [h] = dxfdraw (dxf, clr=[0.7 0.7 0.7; 0.8 0.9 0.99], varargin) dxf = dxfread (dxf); else - error ("dxfdraw: DXF cell array, DXF file name, or DXF drawing struct \ -expected for arg #1 "); + error (["dxfdraw: DXF cell array, DXF file name, or DXF drawing ", ... + "struct expected for arg #1 "]); endif ## parse DXF cell array into a DXF drawing struct diff --git a/inst/dxfparse.m b/inst/dxfparse.m index 15e3f39..cd45bce 100644 --- a/inst/dxfparse.m +++ b/inst/dxfparse.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 @@ -196,7 +196,7 @@ function [dxfo] = dxfparse (dxfi, outstyle=0) y1 = dxfn(find (dxfn(idx(ii):idx(ii+1), 1) == 20) + idx(ii)-1, 2); z1 = dxfn(find (dxfn(idx(ii):idx(ii+1), 1) == 30) + idx(ii)-1, 2); lb = dxfi(find (dxfn(idx(ii):idx(ii+1), 1) == 8) + idx(ii)-1, 2){1}; - if (outsyle == 0) + if (outstyle == 0) ## X, Y [,Z] coordinates XYZp(++is, 1) = x1; XYZ(is, 2) = y1; diff --git a/inst/dxfread.m b/inst/dxfread.m index 6dd0eb3..0238342 100644 --- a/inst/dxfread.m +++ b/inst/dxfread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2016-2020 Philip Nienhuis +## Copyright (C) 2016-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 diff --git a/inst/earthRadius.m b/inst/earthRadius.m index 37f1f11..926b6df 100644 --- a/inst/earthRadius.m +++ b/inst/earthRadius.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 @@ -26,7 +26,7 @@ ## => ans = 6371 ## @end example ## -## @seealso {validateLengthUnit,unitsratio} +## @seealso{validateLengthUnit,unitsratio} ## @end deftypefn function R = earthRadius (unit) diff --git a/inst/ecc2flat.m b/inst/ecc2flat.m index 820b99b..b986014 100644 --- a/inst/ecc2flat.m +++ b/inst/ecc2flat.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/ecc2n.m b/inst/ecc2n.m index 522acfe..b420731 100644 --- a/inst/ecc2n.m +++ b/inst/ecc2n.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 diff --git a/inst/ecef2aer.m b/inst/ecef2aer.m index 1d68623..f694ec8 100644 --- a/inst/ecef2aer.m +++ b/inst/ecef2aer.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -114,6 +114,9 @@ function [az, el, slantRange] = ecef2aer (varargin) spheroid = varargin{7}; angleUnit = varargin{8}; endif + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif x = varargin{1}; y = varargin{2}; @@ -137,8 +140,8 @@ function [az, el, slantRange] = ecef2aer (varargin) if (! all (size (lat0) == size (x)) || ... ! all (size (lon0) == size (y)) || ... ! all (size (alt0) == size (z))) - error ("ecef2aer: non-matching dimensions of observer points and \ -target points"); + error (["ecef2aer: non-matching dimensions of observer points and ", ... + "target points"]); endif endif @@ -146,15 +149,9 @@ target points"); error ("ecef2aer: angleUnit should be one of 'degrees' or 'radians'") endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); - [e, n, u] = ecef2enu (x, y, z, lat0, lon0, alt0, spheroid, angleUnit); + [e, n, u] = ecef2enu (x, y, z, lat0, lon0, alt0, E, angleUnit); [az, el, slantRange] = enu2aer (e, n, u, angleUnit); endfunction diff --git a/inst/ecef2enu.m b/inst/ecef2enu.m index 4d49f19..d67b385 100644 --- a/inst/ecef2enu.m +++ b/inst/ecef2enu.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -69,7 +69,7 @@ ## u = 939.10 ## @end example ## -## @seealso {enu2ecef, ecef2enuv, ecef2geodetic, ecef2ned, ecef2enu, +## @seealso{enu2ecef, ecef2enuv, ecef2geodetic, ecef2ned, ecef2enu, ## referenceEllipsoid} ## @end deftypefn @@ -100,6 +100,9 @@ function [e,n,u] = ecef2enu (varargin) spheroid = varargin{7}; angleUnit = varargin{8}; endif + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif x = varargin{1}; y = varargin{2}; @@ -124,24 +127,18 @@ function [e,n,u] = ecef2enu (varargin) if (! all (size (lat0) == size (x)) || ... ! all (size (lon0) == size (y)) || ... ! all (size (alt0) == size (z))) - error ("ecef2enu: non-matching dimensions of observer points and \ -target points"); + error (["ecef2enu: non-matching dimensions of observer points and ", ... + "target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("ecef2enu: angleUnit should be one of 'degrees' or 'radians'") endif - [x0, y0, z0] = geodetic2ecef (spheroid, lat0, lon0, alt0, angleUnit); + [x0, y0, z0] = geodetic2ecef (E, lat0, lon0, alt0, angleUnit); [e, n, u] = ecef2enuv (x - x0, y - y0, z - z0, lat0, lon0, angleUnit); endfunction diff --git a/inst/ecef2enuv.m b/inst/ecef2enuv.m index 2151304..77b688a 100644 --- a/inst/ecef2enuv.m +++ b/inst/ecef2enuv.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -67,7 +67,7 @@ ## u = 657.11 ## @end example ## -## @seealso {enu2ecefv, ecef2enu, ecef2geodetic, ecef2ned, ecef2enu, +## @seealso{enu2ecefv, ecef2enu, ecef2geodetic, ecef2ned, ecef2enu, ## referenceEllipsoid} ## @end deftypefn diff --git a/inst/ecef2geodetic.m b/inst/ecef2geodetic.m index 4868b38..a9bb337 100644 --- a/inst/ecef2geodetic.m +++ b/inst/ecef2geodetic.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software @@ -25,9 +25,11 @@ ## @itemize ## @item ## @var{spheroid} can be a referenceEllipsoid name or struct (see help for -#3 referenceEllipsoid.m). Unfortunately an EPSG number cannot be accepted. If -## omitted or if an empty string or empty array ('[]') is supplied the WGS84 -## ellipsoid (EPSG 7030) will be selected. +## referenceEllipsoid.m). Unfortunately an EPSG number as first argument +## input as a numeric value cannot be accepted UNLESS ecef2geodetic is called +## with five (5) input arguments. Rather input the number as a character string +## (between quotes). If omitted or if an empty string or empty array ('[]') +## is supplied the WGS84 ellipsoid (EPSG 7030) will be selected. ## ## Inputting @var{spheroid} as 4th argument is accepted but not recommended; ## in that case the @var{lat} and @var{lon} outputs are returned in radians. @@ -109,6 +111,9 @@ function [lat, lon, alt] = ecef2geodetic (varargin) elseif (nargin == 5) ip = 1; spheroid = varargin{1}; + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif angleUnit = varargin{5}; endif X = varargin{ip + 1}; @@ -128,13 +133,7 @@ function [lat, lon, alt] = ecef2geodetic (varargin) error ("ecef2geodetic: angleUnit should be one of 'degrees' or 'radians'") endif - if isempty(spheroid) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); ## Insight from: http://wiki.gis.com/wiki/index.php/Geodetic_system lon = atan2 (Y, X); diff --git a/inst/ecef2ned.m b/inst/ecef2ned.m index ad0198a..23bdfc6 100644 --- a/inst/ecef2ned.m +++ b/inst/ecef2ned.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -69,7 +69,7 @@ ## d = -939.10 ## @end example ## -## @seealso {ned2ecef, ecef2aer, ecef2enu, ecef2enuv, ecef2geodetic, ecef2nedv, +## @seealso{ned2ecef, ecef2aer, ecef2enu, ecef2enuv, ecef2geodetic, ecef2nedv, ## referenceEllipsoid} ## @end deftypefn @@ -87,7 +87,7 @@ function [n, e, d] = ecef2ned (varargin) elseif (nargin == 7) if (isnumeric (varargin{7})) ## EPSG spheroid code - spheroid = varargin{7}; + spheroid = num2str (varargin{7}); elseif (ischar (varargin{7})) if (! isempty (varargin{7}) && ismember (varargin{7}(1), {"r", "d"})) angleUnit = varargin{7}; @@ -126,24 +126,18 @@ function [n, e, d] = ecef2ned (varargin) if (! all (size (lat0) == size (x)) || ... ! all (size (lon0) == size (y)) || ... ! all (size (alt0) == size (z))) - error ("ecef2ned: non-matching dimensions of observer points and \ -target points"); + error (["ecef2ned: non-matching dimensions of observer points and ", ... + "target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("ecef2ned: angleUnit should be one of 'degrees' or 'radians'") endif - [x0, y0, z0] = geodetic2ecef (spheroid, lat0, lon0, alt0, angleUnit); + [x0, y0, z0] = geodetic2ecef (E, lat0, lon0, alt0, angleUnit); [n, e, d] = ecef2nedv (x - x0, y - y0, z - z0, lat0, lon0, angleUnit); endfunction diff --git a/inst/ecef2nedv.m b/inst/ecef2nedv.m index 8caf1cb..826fb99 100644 --- a/inst/ecef2nedv.m +++ b/inst/ecef2nedv.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -67,7 +67,7 @@ ## d = -657.11 ## @end example ## -## @seealso {ned2ecefv, ecef2aer, ecef2enu, ecef2enuv, ecef2ned, ecef2geodetic} +## @seealso{ned2ecefv, ecef2aer, ecef2enu, ecef2enuv, ecef2ned, ecef2geodetic} ## @end deftypefn ## Function adapted by anonymous contributor, see: diff --git a/inst/egm96geoid.m b/inst/egm96geoid.m index 1eb57a0..f15eafe 100644 --- a/inst/egm96geoid.m +++ b/inst/egm96geoid.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 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 @@ -22,8 +22,8 @@ ## ## The input values @var{lat}, @var{lon} are in (decimal) degrees. Scalar ## values as well as vectors/2D arrays are accepted are accepted, in the -## latter case the dimensions of @var{lon}, @var{lat} should match. The -## are wrapped in the latitude = (-90:90) and longitude (-80:180) intervals. +## latter case the dimensions of @var{lon} and @var{lat} should match. They +## are wrapped in the latitude = [-90:90] and longitude [-180:180] intervals. ## ## The optional third input @var{method} defines the interpolation method ## and can be one of the following: "nearest" (default), "linear", @@ -40,7 +40,7 @@ ## (721, 721) of the base grid correspond to (Lon, Lat) = (0, 90), ## (180, 90), (360, 90) = (0, 90), and (180, -90), respectively. ## -## @var{hgt}, the egm96geoid output value(s) are in meters relative to +## @var{hgt}, the egm96geoid output value(s), are in meters relative to ## the WGS84 standard ellipsoid. ## ## @example @@ -57,7 +57,7 @@ ## while near the poles the grid cells are effectively more needle-like ## triangles with a base of merely 121 m and a height of 27.76 km. ## -## @seealso{interp2,referenceEllipsoid} +## @seealso{interp2, referenceEllipsoid} ## @end deftypefn ## Author: Philip Nienhuis <prnienhuis at users.sf.net> diff --git a/inst/enu2aer.m b/inst/enu2aer.m index 27d3a8f..5dcb3c5 100644 --- a/inst/enu2aer.m +++ b/inst/enu2aer.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without modification, are permitted ## provided that the following conditions are met: diff --git a/inst/enu2ecef.m b/inst/enu2ecef.m index 73b0744..39407d3 100644 --- a/inst/enu2ecef.m +++ b/inst/enu2ecef.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -77,7 +77,7 @@ ## z = 4.4884e+06 ## @end example ## -## @seealso {ecef2enu, enu2aer, enu2ecefv, enu2geodetic, enu2uvw} +## @seealso{ecef2enu, enu2aer, enu2ecefv, enu2geodetic, enu2uvw} ## @end deftypefn ## Function adapted by anonymous contributor, see: @@ -94,7 +94,7 @@ function [x, y, z] = enu2ecef (varargin) elseif (nargin == 7) if (isnumeric (varargin{7})) ## EPSG spheroid code - spheroid = varargin{7}; + spheroid = num2str (varargin{7}); elseif (ischar (varargin{7})) if (! isempty (varargin{7}) && ismember (varargin{7}(1), {"r", "d"})) angleUnit = varargin{7}; @@ -135,21 +135,15 @@ function [x, y, z] = enu2ecef (varargin) if (! all (size (lat) == size (e)) || ... ! all (size (lon) == size (n)) || ... ! all (size (alt) == size (u))) - error ("enu2ecef: non-matching dimensions of observer points and \ -target points"); + error (["enu2ecef: non-matching dimensions of observer points and ", ... + "target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid)) - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); ## Origin of the local system in geocentric coordinates. - [x0, y0, z0] = geodetic2ecef (spheroid, lat, lon, alt, angleUnit); + [x0, y0, z0] = geodetic2ecef (E, lat, lon, alt, angleUnit); ## Rotating ENU to ECEF [dx, dy, dz] = enu2uvw (e, n, u, lat, lon, angleUnit); ## Origin + offset from origin equals position in ECEF diff --git a/inst/enu2ecefv.m b/inst/enu2ecefv.m index 5665606..0b0fc6a 100644 --- a/inst/enu2ecefv.m +++ b/inst/enu2ecefv.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -69,7 +69,7 @@ ## w = 1000.0 ## @end example ## -## @seealso {ecef2enuv, enu2aer, enu2ecef, enu2geodetic, enu2uvw} +## @seealso{ecef2enuv, enu2aer, enu2ecef, enu2geodetic, enu2uvw} ## @end deftypefn ## Function adapted by anonymous contributor, see: diff --git a/inst/enu2geodetic.m b/inst/enu2geodetic.m index 060b276..11d87f5 100644 --- a/inst/enu2geodetic.m +++ b/inst/enu2geodetic.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -83,7 +83,7 @@ ## alt = 1139.7 ## @end example ## -## @seealso {geodetic2enu, enu2aer, enu2ecef, enu2ecefv, enu2uvw, egm96geoid, +## @seealso{geodetic2enu, enu2aer, enu2ecef, enu2ecefv, enu2uvw, egm96geoid, ## referenceEllipsoid} ## @end deftypefn @@ -99,7 +99,7 @@ function [lat, lon, alt] = enu2geodetic (varargin) elseif (nargin == 7) if (isnumeric (varargin{7})) ## EPSG spheroid code - spheroid = varargin{7}; + spheroid = num2str (varargin{7}); elseif (ischar (varargin{7})) if (! isempty (varargin{7}) && ismember (varargin{7}(1), {"r", "d"})) angleUnit = varargin{7}; @@ -140,21 +140,15 @@ function [lat, lon, alt] = enu2geodetic (varargin) if (! all (size (lat0) == size (e)) || ... ! all (size (lon0) == size (n)) || ... ! all (size (alt0) == size (u))) - error ("enu2geodetic: non-matching dimensions of non-scalar observer \ -points and target points"); + error (["enu2geodetic: non-matching dimensions of non-scalar ", ... + "observer points and target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid)) - E = referenceEllipsoid (spheroid); - endif + E = sph_chk (spheroid); - [x, y, z] = enu2ecef(e, n, u, lat0, lon0, alt0, spheroid, angleUnit); - [lat, lon, alt] = ecef2geodetic(spheroid, x, y, z, angleUnit); + [x, y, z] = enu2ecef (e, n, u, lat0, lon0, alt0, E, angleUnit); + [lat, lon, alt] = ecef2geodetic (E, x, y, z, angleUnit); endfunction diff --git a/inst/enu2uvw.m b/inst/enu2uvw.m index 9cc286e..9dbf801 100644 --- a/inst/enu2uvw.m +++ b/inst/enu2uvw.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (C) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2019-2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (C) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2019-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -41,7 +41,7 @@ ## ## @var{u}, @var{v}, @var{w}: coordinates of point(s) (meters). ## -## @seealso {enu2aer, enu2ecef, enu2ecefv, enu2geodetic} +## @seealso{enu2aer, enu2ecef, enu2ecefv, enu2geodetic} ## @end deftypefn ## Function adapted by anonymous contributor, see: @@ -67,8 +67,8 @@ function [u, v, w] = enu2uvw (east, n, up, lat0, lon0, angleUnit = "degrees") ## Check if for each target point a matching obsrver point is given if (! all (size (lat0) == size (east)) || ... ! all (size (lon0) == size (n))) - error ("enu2uvw: non-matching dimensions of observer points and \ -target points"); + error (["enu2uvw: non-matching dimensions of observer points and ", ... + "target points"]); endif endif diff --git a/inst/extractfield.m b/inst/extractfield.m index b18aa29..48954e6 100644 --- a/inst/extractfield.m +++ b/inst/extractfield.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/flat2ecc.m b/inst/flat2ecc.m index 6e3137d..cda57da 100644 --- a/inst/flat2ecc.m +++ b/inst/flat2ecc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/fromDegrees.m b/inst/fromDegrees.m index 01575b6..127850b 100644 --- a/inst/fromDegrees.m +++ b/inst/fromDegrees.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/fromRadians.m b/inst/fromRadians.m index 7725068..af32ff2 100644 --- a/inst/fromRadians.m +++ b/inst/fromRadians.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/gc2sc.m b/inst/gc2sc.m index c9e037e..ba54c4f 100644 --- a/inst/gc2sc.m +++ b/inst/gc2sc.m @@ -1,5 +1,4 @@ -## Copyright (C) 2020 The Octave Project Developers -## Copyright (C) 2014 Alfredo Foltran <alfoltran@gmail.com> +## Copyright (C) 2020-2022 The Octave Project Developers ## ## 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 @@ -119,7 +118,7 @@ function [lat, lon, radius] = gc2sc (varargin); range = pi / 2; radius = repmat (range, size (vect(:, 1))); - [lat, lon] = reckon (vect(:, 1), vect(:, 2), range, vect(:, 3) .+ range, "radians"); + [lat, lon] = reckon (vect(:, 1), vect(:, 2), range, vect(:, 3) + range, "radians"); if (abs (abs (lat) - pi / 2) < (4 * eps)) ## NOTE: You are at the pole many longitudes so NaN should be used diff --git a/inst/gcxgc.m b/inst/gcxgc.m index db43df6..5205697 100644 --- a/inst/gcxgc.m +++ b/inst/gcxgc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 The Octave Project Developers +## Copyright (C) 2022 The Octave Project Developers ## ## 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 @@ -142,7 +142,7 @@ function [lat, lon, idl] = gcxgc (varargin) ## 1. Find circles with azimuth = lat == 0 (as those ARE on equator) iaz0 = double (abs (rem (vect(:, [3 6]) + pi/2, pi)) < 2 * eps & ... abs (vect(:, [1 4])) < 2 * eps); - iaz0 = (iaz0(1:nv) .+ iaz0(nv+1:end))'; # Sum iaz0 by rows + iaz0 = (iaz0(1:nv) + iaz0(nv+1:end))'; # Sum iaz0 by rows ## vect(iaz0==2, :) => two great circles equaling equators => set to NaN & skip lat(iaz0 > 1.5, :) = NaN; lon(iaz0 > 1.5, :) = NaN; @@ -162,7 +162,7 @@ function [lat, lon, idl] = gcxgc (varargin) id2 = loni2(:, 1) <= 0; loni2(id2, 1) = loni2(id2, 2); ## 4. Find out which loni's coincide - idl = (abs (loni1(:, 1) .- loni2(:, 1)) < 2 * eps)(:, 1) & abs (id1 .- id2); + idl = (abs (loni1(:, 1) - loni2(:, 1)) < 2 * eps)(:, 1) & abs (id1 - id2); ## 5. Set output relating to coinciding great circles to NaN, NaN lat(iazx(idl), :) = NaN; lon(iazx(idl), :) = NaN; @@ -189,18 +189,18 @@ endfunction function [lat, lon] = get_intscs (vect) ## Algorithm from https://www.movable-type.co.uk/scripts/latlong-vectors.html#intersection - c1(:, 1) = sin (vect(:, 2)) .* cos (vect(:, 3)) .- sin (vect(:, 1)) .* ... + c1(:, 1) = sin (vect(:, 2)) .* cos (vect(:, 3)) - sin (vect(:, 1)) .* ... cos (vect(:, 2)) .* sin (vect(:, 3)); - c1(:, 2) = -cos (vect(:, 2)) .* cos (vect(:, 3)) .- sin (vect(:, 1)) .* ... + c1(:, 2) = -cos (vect(:, 2)) .* cos (vect(:, 3)) - sin (vect(:, 1)) .* ... sin (vect(:, 2)) .* sin (vect(:, 3)); c1(:, 3) = cos (vect(:, 1)) .* sin (vect(:, 3)); - c2(:, 1) = sin (vect(:, 5)) .* cos (vect(:, 6)) .- sin (vect(:, 4)) .* ... + c2(:, 1) = sin (vect(:, 5)) .* cos (vect(:, 6)) - sin (vect(:, 4)) .* ... cos (vect(:, 5)) .* sin (vect(:, 6)); - c2(:, 2) = -cos (vect(:, 5)) .* cos (vect(:, 6)) .- sin (vect(:, 4)) .* ... + c2(:, 2) = -cos (vect(:, 5)) .* cos (vect(:, 6)) - sin (vect(:, 4)) .* ... sin (vect(:, 5)) .* sin (vect(:, 6)); c2(:, 3) = cos (vect(:, 4)) .* sin (vect(:, 6)); diff --git a/inst/gcxsc.m b/inst/gcxsc.m index 556762e..0c56da8 100644 --- a/inst/gcxsc.m +++ b/inst/gcxsc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 The Octave Project Developers +## Copyright (C) 2022 The Octave Project Developers ## ## 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 diff --git a/inst/geo2auth.m b/inst/geo2auth.m new file mode 100644 index 0000000..f47720f --- /dev/null +++ b/inst/geo2auth.m @@ -0,0 +1,123 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{xi} =} geo2auth (@var{phi}) +## @deftypefnx {} {@var{xi} =} geo2auth (@var{phi}, @var{spheroid}) +## @deftypefnx {} {@var{xi} =} geo2auth (@var{phi}, @var{spheroid}, @var{angleUnit}) +## Returns the authalic latitude given geodetic latitude @var{phi} +## +## Input +## @itemize +## @item +## @var{phi}: the geodetic latitude +## @end itemize +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid name, EPSG number, or +## parameter struct: the default is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first charachter will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{xi}: the authalic latitude +## @end itemize +## +## Example +## @example +## geo2auth(45) +## ans = 44.872 +## @end example +## @seealso{geo2con, geo2iso, geo2rect} +## @end deftypefn + + +function [xi] = geo2auth (phi, spheroid ="", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (phi) || ! isreal (phi)) + error ("geo2auth: numeric input expected"); + endif + + if (isempty (spheroid)) + E = referenceEllipsoid ("wgs84"); + else + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif + E = sph_chk (spheroid); + endif + + if (! ischar (angleUnit)) + error ("geo2auth: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (phi) > 90)) + error ("geo2auth: azimuth value(s) out of acceptable range [-90, 90]") + endif + phi = deg2rad (phi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (phi) > pi / 2)) + error ("geo2auth: azimuth value(s) out of acceptable range (-pi/2, pi/2)"); + endif + else + error ("geo2auth: illegal input for 'angleUnit'"); + endif + + + ecc = E.Eccentricity; + + em = 1 - ecc .^ 2; + sp = sin (phi); + + q_p = 1 + ((em ./ ecc) .* atanh(ecc)); + + q = (em .* sp) ./ (1 - ecc .^ 2 .* sp .^2) + ... + ((em ./ ecc) .* atanh(ecc .* sp)); + + xi = asin (q ./ q_p); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + xi = rad2deg (xi); + endif + +endfunction + +%!test +%! phi = [0:15:90]; +%! xi = geo2auth (phi); +%! Z = degrees2dm(xi - phi); +%! check = [0; ... +%! -3.84258303680237; ... +%! -6.66017793242659; ... +%! -7.69782759396364; ... +%! -6.67286580689009; ... +%! -3.85527095139878; ... +%! 0]; +%! assert (Z(:,2), check, 1e-6) + +%!error <numeric> geo2auth ("s") +%!error <numeric> geo2auth (5i) + diff --git a/inst/geo2con.m b/inst/geo2con.m new file mode 100644 index 0000000..9603daa --- /dev/null +++ b/inst/geo2con.m @@ -0,0 +1,114 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{chi} =} geo2con (@var{phi}) +## @deftypefnx {} {@var{chi} =} geo2con (@var{phi}, @var{spheroid}) +## @deftypefnx {} {@var{chi} =} geo2con (@var{phi}, @var{spheroid}, @var{angleUnit}) +## Returns the conformal latitude given geodetic latitude @var{phi} +## +## Input +## @itemize +## @item +## @var{phi}: the geodetic latitude +## @end itemize +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid name, EPSG number, or +## parameter struct: the deualt is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first charachter will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{chi}: the conformal latitude +## @end itemize +## +## Example +## @example +## geo2con(45) +## ans = 44.808 +## @end example +## @seealso{geo2auth, geo2iso, geo2rect} +## @end deftypefn + +function [chi] = geo2con (phi, spheroid ="", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (phi) || ! isreal (phi)) + error ("geo2con: numeric input expected"); + endif + + if (isempty (spheroid)) + E = referenceEllipsoid ("wgs84"); + else + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif + E = sph_chk (spheroid); + endif + + if (! ischar (angleUnit)) + error ("geo2con: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (phi) > 90)) + error ("geo2con: azimuth value(s) out of acceptable range [-90, 90]") + endif + phi = deg2rad (phi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (phi) > pi / 2)) + error("geo2con: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("geo2con: illegal input for 'angleUnit'"); + endif + + ecc = E.Eccentricity; + + chi = atan (sinh (asinh (tan (phi)) - ecc .* atanh (ecc .* sin (phi)))); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + chi = rad2deg (chi); + endif + +endfunction + + +%!test +%! phi = [0:15:90]; +%! chi = geo2con (phi); +%! Z = degrees2dm(chi - phi); +%! check = [0; ... +%! -5.755543956550260; ... +%! -9.979077451140980; ... +%! -11.53895663467140; ... +%! -10.00703049899710; ... +%! -5.783497189944170; ... +%! 0]; +%! assert (Z(:,2), check, 1e-6) + +%!error <numeric> geo2con ("s") +%!error <numeric> geo2con (5i) + diff --git a/inst/geo2iso.m b/inst/geo2iso.m new file mode 100644 index 0000000..9777889 --- /dev/null +++ b/inst/geo2iso.m @@ -0,0 +1,120 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{psi} =} geo2iso (@var{phi}) +## @deftypefnx {} {@var{psi} =} geo2iso (@var{phi}, @var{spheroid}) +## @deftypefnx {} {@var{psi} =} geo2iso (@var{phi}, @var{spheroid}, @var{angleUnit}) +## Return the isometric latitude given geodetic latitude @var{phi}. +## +## Input +## @itemize +## @item +## @var{phi}: the geodetic latitude. Can be a scalar value, a vector or an +## ND-array. +## @end itemize +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid. For admissible values see +## `referenceEllipsoid.m`. The default ellipsoid is WGS84. +## is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first character will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{psi}: the isometric latitude(s), same shape as @var{phi}. +## @end itemize +## +## Example +## @example +## geo2iso (45) +## ans = 50.227 +## @end example +## @seealso{geo2auth, geo2con, geo2rect, iso2geo} +## @end deftypefn + +function [psi] = geo2iso (phi, spheroid = "", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (phi) || ! isreal (phi)) + error ("geo2iso: numeric input expected"); + endif + + E = sph_chk (spheroid); + + if (! ischar (angleUnit)) + error ("geo2iso: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (phi) > 90)) + error ("geo2iso: azimuth value(s) out of acceptable range [-90, 90]") + endif + phi = deg2rad (phi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (phi) > pi / 2)) + error("geo2iso: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("geo2iso: illegal input for 'angleUnit'"); + endif + + ## Some previously tried solutions: + + ## https://en.wikipedia.org/wiki/Latitude#Isometric_latitude + ## ecc = E.Eccentricity; + ## psi = asinh (tan (phi)) - ecc .* atanh (ecc .* sin (phi)); + + ## From Snyder's "Map Projections - A Working Manual" [pg 15]. Eq (3-7) + ## psi = log(tan(pi/4+phi/2)*((1-ecc.*sin(phi))/(1+ecc.*sin(phi)))^(ecc/2)); + + ## From Snyder's "Map Projections - A Working Manual" [pg 15]. Eq (3-8) + chi = geo2con (phi, E, "radians"); + psi = log (tan (pi / 4 + chi / 2)); + + psi(psi < -5.24) = -Inf; + psi(psi > 5.24) = Inf; + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + psi = rad2deg (psi); + endif + +endfunction + + +%!test +%! psi = geo2iso (45); +%! assert (psi, 50.227465817, 1e-6) + +%!test +%! psi = geo2iso ([-90 90]); +%! assert (psi, [-Inf Inf]) + +%!test +%! chi = geo2iso (iso2geo ([-90 : 10: 0; 0 : 10 : 90]), "wgs84"); +%! assert (chi, [-90 : 10: 0; 0 : 10 : 90], 1e-6); + +%!error <numeric> geo2iso ("s") +%!error <numeric> geo2iso (5i) + diff --git a/inst/geo2rect.m b/inst/geo2rect.m new file mode 100644 index 0000000..de4019c --- /dev/null +++ b/inst/geo2rect.m @@ -0,0 +1,107 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{mu} =} geo2rect (@var{phi}) +## @deftypefnx {} {@var{mu} =} geo2rect (@var{phi}, @var{spheroid}) +## @deftypefnx {} {@var{mu} =} geo2rect (@var{phi}, @var{spheroid}, @var{angleUnit}) +## Returns the rectifying latitude given geodetic latitude @var{phi} +## +## Input +## @itemize +## @item +## @var{phi}: the geodetic latitude +## @end itemize +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid name, EPSG code, or parameter +## struct: default is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first character will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{mu}: the rectifying latitude +## @end itemize +## +## Example +## @example +## geo2rect(45) +## ans = 44.856 +## @end example +## @seealso{geo2auth, geo2con, geo2iso} +## @end deftypefn + +function [mu] = geo2rect (phi, spheroid ="", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (phi) || ! isreal (phi)) + error ("geo2rect: numeric input expected"); + endif + + if (isempty (spheroid)) + E = referenceEllipsoid ("wgs84"); + else + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif + E = sph_chk (spheroid); + endif + + if (! ischar (angleUnit)) + error ("geo2rect: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (phi) > 90)) + error ("geo2rect: azimuth value(s) out of acceptable range [-90, 90]") + endif + phi = deg2rad (phi); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (phi) > pi / 2)) + error("geo2rect: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("geo2rect: illegal input for 'angleUnit'"); + endif + + m = meridianarc (0, phi, E); + m_p = meridianarc (0, pi / 2, E); + + mu = pi / 2 .* m ./ m_p; + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + mu = rad2deg (mu); + endif + +endfunction + + +%!test +%! mu = geo2rect (45); +%! Z = degrees2dm(mu - 45); +%! assert (Z(:,2), -8.65908066558504, 1e-6); + +%!error <numeric> geo2rect ("s") +%!error <numeric> geo2rect (5i) + diff --git a/inst/geocentricLatitude.m b/inst/geocentricLatitude.m index 6f79b14..2d90d0e 100644 --- a/inst/geocentricLatitude.m +++ b/inst/geocentricLatitude.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/geodetic2aer.m b/inst/geodetic2aer.m index 6d8df9a..f1772c8 100644 --- a/inst/geodetic2aer.m +++ b/inst/geodetic2aer.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2018 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -73,7 +73,7 @@ ## slantRange = 846.65 ## @end example ## -## @seealso {aer2geodetic, geodetic2ecef, geodetic2enu, geodetic2ned, +## @seealso{aer2geodetic, geodetic2ecef, geodetic2enu, geodetic2ned, ## referenceEllipsoid} ## @end deftypefn @@ -129,24 +129,21 @@ function [az, el, slantRange] = geodetic2aer (varargin) if (! all (size (lat0) == size (lat)) || ... ! all (size (lon0) == size (lon)) || ... ! all (size (alt0) == size (alt))) - error ("geodetic2aer: non-matching dimensions of geodetic locations \ -and AER target points"); + error (["geodetic2aer: non-matching dimensions of geodetic locations ", ... + "and AER target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("geodetic2aer: angleUnit should be one of 'degrees' or 'radians'") endif - [e, n, u] = geodetic2enu (lat, lon, alt, lat0, lon0, alt0, spheroid, angleUnit); + [e, n, u] = geodetic2enu (lat, lon, alt, lat0, lon0, alt0, E, angleUnit); [az, el, slantRange] = enu2aer (e, n, u, angleUnit); endfunction diff --git a/inst/geodetic2ecef.m b/inst/geodetic2ecef.m index af45692..132402b 100644 --- a/inst/geodetic2ecef.m +++ b/inst/geodetic2ecef.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software @@ -128,13 +128,10 @@ function [X, Y, Z] = geodetic2ecef (varargin) error ("geodetic2ecef: angleUnit should be one of 'degrees' or 'radians'") endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); if (strncmpi (lower (angleUnit), "r", 1) == 1) c_p = cos (lat); diff --git a/inst/geodetic2enu.m b/inst/geodetic2enu.m index e591019..3d82082 100644 --- a/inst/geodetic2enu.m +++ b/inst/geodetic2enu.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -73,7 +73,7 @@ ## u = 799.99 ## @end example ## -## @seealso {enu2geodetic, geodetic2aer, geodetic2ecef, geodetic2ned, +## @seealso{enu2geodetic, geodetic2aer, geodetic2ecef, geodetic2ned, ## egm96geoid, referenceEllipsoid} ## @end deftypefn @@ -124,13 +124,10 @@ function [e, n, u] = geodetic2enu (varargin) error ("geodetic2enu: non-matching dimensions of ECEF inputs."); endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("geodetic2enu: angleUnit should be one of 'degrees' or 'radians'") @@ -139,9 +136,9 @@ function [e, n, u] = geodetic2enu (varargin) [x1, y1, z1] = geodetic2ecef (E, lat, lon, alt, angleUnit); [x2, y2, z2] = geodetic2ecef (E, lat0, lon0, alt0, angleUnit); - dx = x1 .- x2; - dy = y1 .- y2; - dz = z1 .- z2; + dx = x1 - x2; + dy = y1 - y2; + dz = z1 - z2; [e, n, u] = ecef2enuv (dx, dy, dz, lat0, lon0, angleUnit); diff --git a/inst/geodetic2ned.m b/inst/geodetic2ned.m index 214fef7..8c60699 100644 --- a/inst/geodetic2ned.m +++ b/inst/geodetic2ned.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -46,7 +46,7 @@ ## geoid and the WGS84 reference ellipsoid. ## ## @item -## @var{spheroid} (optional): a user-specified sheroid (see referenceEllipsoid); +## @var{spheroid} (optional): a user-specified spheroid (see referenceEllipsoid); ## it can be omitted or given as an empty string or empty numeric array('[]'), ## in which cases WGS84 will be selected as default spheroid. ## @@ -74,7 +74,7 @@ ## u = -799.99 ## @end example ## -## @seealso {ned2geodetic, geodetic2aer, geodetic2ecef, geodetic2enu, +## @seealso{ned2geodetic, geodetic2aer, geodetic2ecef, geodetic2enu, ## egm96geoid, referenceEllipsoid} ## @end deftypefn @@ -125,13 +125,10 @@ function [n, e, d] = geodetic2ned (varargin) error ("geodetic2ned: non-matching dimensions of ECEF inputs."); endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("geodetic2ned: angleUnit should be one of 'degrees' or 'radians'") @@ -140,9 +137,9 @@ function [n, e, d] = geodetic2ned (varargin) [x1, y1, z1] = geodetic2ecef (E, lat, lon, alt, angleUnit); [x2, y2, z2] = geodetic2ecef (E, lat0, lon0, alt0, angleUnit); - dx = x1 .- x2; - dy = y1 .- y2; - dz = z1 .- z2; + dx = x1 - x2; + dy = y1 - y2; + dz = z1 - z2; [n, e, d] = ecef2nedv (dx, dy, dz, lat0, lon0, angleUnit); diff --git a/inst/geodeticLatitudeFromGeocentric.m b/inst/geodeticLatitudeFromGeocentric.m index 04aad96..dad7114 100644 --- a/inst/geodeticLatitudeFromGeocentric.m +++ b/inst/geodeticLatitudeFromGeocentric.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/geodeticLatitudeFromParametric.m b/inst/geodeticLatitudeFromParametric.m index ef7203f..7c91927 100644 --- a/inst/geodeticLatitudeFromParametric.m +++ b/inst/geodeticLatitudeFromParametric.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/geodeticarc.m b/inst/geodeticarc.m new file mode 100644 index 0000000..ea9e119 --- /dev/null +++ b/inst/geodeticarc.m @@ -0,0 +1,244 @@ +## Copyright (C) 2014-2022 Alfredo Foltran <alfoltran@gmail.com> +## Copyright (C) 2021-2022 Philip Nienhuis <prnienhuis@users.sf.net> +## +## 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} = geodeticarc(@var{pt1}, @var{pt2}) +## @deftypefnx {Function File} {} @var{dist} = geodeticarc(@var{pt1}, @var{pt2}, @var{ellipsoid}) +## @deftypefnx {Function File} {[@var{dist}, @var{az}] = } {geodeticarc(@var{pt1}, @var{pt2})} +## @deftypefnx {Function File} {[@var{dist}, @var{az}] = } {geodeticarc(@var{pt1}, @var{pt2}, @var{ellipsoid})} +## Calculates the distance (in meters) between two (sets of) locations on +## an ellipsoid. +## +## The formula devised by Thaddeus Vincenty (1975) is used with an accurate +## ellipsoidal model of the earth (@var{ellipsoid}). @* +## Note: for antipodal points (within 0.5 degree) Vincenty's formulae are +## known to be inaccurate and may even break down. +## +## Inputs: +## +## @itemize +## @item +## @var{pt1} and @var{pt2} are two-column matrices of the form @* +## [latitude longitude]. +## The units for the input coordinates angles must be degrees. +## +## @item +## Optional argument @var{ellipsoid} defines the reference ellipsoid to use. +## The default ellipsoidal model is 'WGS84', which is the globally most +## accurate model. +## @end itemize +## +## Outputs: +## +## @itemize +## @item +## @var{dist} is the computed distance between @var{pt1} and @var{pt2} in +## meters, computed along the shortest geodesic. +## +## @item +## @var{az} is a 2-column array of starting and ending azimuths of the +## geodesics in the direction from @var{pt1} to @var{pt2}, in degrees relative +## to the North (clockwise). +## @end itemize +## +## Examples: +## @example +## >> geodeticarc ([37, -76], [37, -9]) +## ans = 5830081.06 +## >> geodeticarc ([37, -76], [67, -76], referenceEllipsoid (7019)) +## ans = 3337842.87 +## @end example +## +## @seealso{distance, geodeticfwd, referenceEllipsoid, vincenty} +## @end deftypefn + +function [dist, az] = geodeticarc (varargin) + + if (nargin < 2) + error ("geodeticarc: too few arguments"); + elseif (! (all (cellfun ("isnumeric", varargin(1:2))) && ... + cellfun ("isreal", varargin(1:2)))) + error ("geodeticarc: first 2 arguments should be real numeric nx2 arrays"); + elseif (nargin == 2) + ellipsoid = referenceEllipsoid ("wgs84"); + elseif (isstruct (varargin{3})) + ellipsoid = varargin{3}; + elseif (isnumeric (varargin{3}) && numel (varargin{3}) == 2) + ## Might be SemimajorAxis and Flattening + ellipsoid.SemimajorAxis = varargin{3}(1); + ellipsoid.Flattening = varargin{3}(2); + ecc = flat2ecc (varargin{3}(2)); + ellipsoid.SemiminorAxis = minaxis (varargin{3}(1), ecc); + else + ellipsoid = referenceEllipsoid (varargin{3}); + endif + + ## Check & process numeric input data sizes + isnv = ! cellfun ("isvector", varargin(1:2)); + if (any (isnv)) + ## Check dimensions + cdim = cellfun (@(x) numel (size (x)), varargin(isnv)); + if (any (cdim != 2)) + error ("geodeticarc: first 2 arguments should be real numeric nx2 arrays"); + endif + ## Check for proper orientation and expand any scalar args, if required + for ii=1:2 + vsz = size (varargin{ii}); + if (vsz(2) != 2) ## orientation + error ("geodeticarc: first 2 arguments should be real numeric nx2 arrays"); + else ## nr.of points + npt(ii) = vsz(1); + endif + endfor + if (any (npt > 1)) + if (any (npt == 1)) + varargin{npt == 1} = repmat (varargin{npt == 1}, npt(npt != 1), 1); + elseif (! (npt(1) - npt(2) == 0)) + error (["geodeticarc: at least one of first 2 inputs must be a 1x2 vector\n" ... + " or both must have the same dimensions\n"]); + endif + endif + elseif (! all (cellfun ("isrow", varargin(1:2)))) + error ("geodeticarc: first 2 arguments should be real numeric nx2 arrays"); + endif + + ## Do the actual work. Based on Alfredo Foltran's vincenty.m + ## https://github.com/alfoltranteam/octave-map/blob/master/src/vincenty.m + ## (In turn based on Vicenty.T (1975) Direct and inverse solutions of + ## geodesics on the ellipsoid with application of nested equations) + + major = ellipsoid.SemimajorAxis; + minor = ellipsoid.SemiminorAxis; + f = ellipsoid.Flattening; + ## Avoid confusion of length units imposed by ellipsoid, standardize on meters + if (isfield (ellipsoid, "LengthUnit") && ! isempty (ellipsoid.LengthUnit)) + major *= unitsratio ("meters", ellipsoid.LengthUnit); + minor *= unitsratio ("meters", ellipsoid.LengthUnit); + endif + + iter_limit = 20; + + pt1 = deg2rad (varargin{1}); + pt2 = deg2rad (varargin{2}); + + [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; + ## Keep track of which sigmas have converged + lit = true (size (lambda_lng)); + dist = NaN (size (lambda_lng)); + ## Preallocate to avoid dimension changes avoid shadowing function "sigma" + sin_sigma = cos_sigma = sin_alpha = cos_sq_alpha = cos2_sigma_m = C = ... + gsigma = zeros (size (lambda_lng)); + do + i++; + [sin_lambda_lng cos_lambda_lng] = deal (sin (lambda_lng), cos (lambda_lng)); + sin_sigma(lit) = sqrt ((cos_reduced2(lit) .* sin_lambda_lng(lit)) .^ 2 + ... + (cos_reduced1(lit) .* sin_reduced2(lit) - ... + sin_reduced1(lit) .* cos_reduced2(lit) .* ... + cos_lambda_lng(lit)) .^ 2); + + dist(find (abs (sin_sigma < eps))) = 0; + lit(find (abs (sin_sigma < eps))) = false; + + cos_sigma(lit) = (sin_reduced1(lit) .* sin_reduced2(lit) + ... + cos_reduced1(lit) .* cos_reduced2(lit) .* cos_lambda_lng(lit)); + gsigma(lit) = atan2 (sin_sigma(lit), cos_sigma(lit)); + sin_alpha(lit) = (cos_reduced1(lit) .* cos_reduced2(lit) .* ... + sin_lambda_lng(lit) ./ sin_sigma(lit)); + cos_sq_alpha(lit) = 1 - sin_alpha(lit) .^ 2; + + if (abs (cos_sq_alpha(lit) > eps)) + cos2_sigma_m(lit) = cos_sigma(lit) - 2 .* (sin_reduced1(lit) ... + .* sin_reduced2(lit) ./ cos_sq_alpha(lit)); + else + cos2_sigma_m(lit) = 0.0; ## Equatorial line + endif + + C(lit) = f / 16.0 * cos_sq_alpha(lit) .* (4 + f * (4 - 3 * cos_sq_alpha(lit))); + + lambda_prime = lambda_lng; + lambda_lng(lit) = (delta_lng(lit) + (1 - C(lit)) * f .* sin_alpha(lit) .* (gsigma(lit) + ... + C(lit) .* sin_sigma(lit) .* (cos2_sigma_m(lit) + C(lit) .* cos_sigma(lit) .* ... + (-1 + 2 * cos2_sigma_m(lit) .^ 2)))); + ## Which inputs haven't converged yet + lit(lit) = abs (lambda_lng(lit) - lambda_prime(lit)) > 10e-12 | abs (sin_sigma(lit) < eps); + ## printf ("%2d ", i, lit); printf ("\n"); ## Debug, remove when tested + until (all (! lit) || i > iter_limit); + + if (i > iter_limit) + warning ("Inverse vincenty's formulae failed to converge for some inputs!"); + sin_sigma(lit) = NaN; + 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(isnan (dist)) = (minor * A .* (gsigma - delta_sigma))(isnan (dist)); + + 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 alpha2]); + endif + +endfunction + + +%!demo +%! lgts = geodeticarc ([[0:1:89]', zeros(90, 1)], [[1:1:90]', zeros(90, 1)]) / 60; +%! plot (0:89, lgts); +%! axis tight; +%! grid on; +%! hold on; +%! plot ([0 89], [1852 1852], "k", "linestyle", "-."); +%! plot ([45 45], [min(lgts) max(lgts)], "m", "linestyle", "-."); +%! title ("Arcminute length vs. Latitude"); +%! xlabel ("Latitude (degrees)", "FontWeight", "bold"); +%! ylabel ("Arcminute length (m)", "FontWeight", "bold"); +%! text (5, 1852.5, "Nautical mile (1852 m)"); + + +%!error <should be real numeric nx2> geodeticarc ("a", [0 0]); +%!error <should be real numeric nx2> geodeticarc ([0 0], [1; 1]); +%!error <should be real numeric nx2> geodeticarc ([0 0; 0 0; 0 0], [1; 1]); +%!error <should be real numeric nx2> geodeticarc ([1 2 3; 2 2 2], [0 0]); +%!error <at least one> geodeticarc ([1 3 2; 2 2 2]', [0 0; 5 6]) + +%!test +%! [dst, az] = geodeticarc ([37, -76], [37, -9; 67 -76; 37, -76]); +%! assert (dst, [5830081.06; 3337842.871; 0], 1e-2); +%! assert (mean (az'), [90, 0, 0,], 1e-7); + diff --git a/inst/geodeticfwd.m b/inst/geodeticfwd.m new file mode 100644 index 0000000..a2e98cb --- /dev/null +++ b/inst/geodeticfwd.m @@ -0,0 +1,344 @@ +## Copyright (C) 2021-2022 Philip Nienhuis <prnienhuis at users.sf.net> +## Copyright (C) 2014-2022 Alfredo Foltran <alfoltran at gmail.com> +## +## 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{lato}, @var{lono}, @var{azo}] = } {geodeticfwd(@var{lat}, @var{lon}, @var{range}, @var{azi})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {geodeticfwd(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {geodeticfwd(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{angleUnit})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {geodeticfwd(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{ellipsoid})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {geodeticfwd(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim}, @var{angleUnit}, @var{ellipsoid})} +## Compute the coordinates of the end-point of a displacement along a geodesic. +## +## Inputs: +## +## @itemize +## @item +## @var{lat}, @var{lon}: coordinates (latitude, longitude) of the starting point. +## +## @item +## @var{range}: displacement along a specified geodesic (angle or length). +## +## @item +## @var{azi}: direction of the displacement relative (clockwise) to the North. +## +## @indent +## All these inputs can be scalars, vectors or 2D/ND arrays. If any input is +## a vector or 2D or ND array, all other inputs MUST be either scalars, OR +## vectors or arrays of the exact same size, OR scalars. Scalars will be +## automatically expanded to the size of the input vectors/arrays. +## +## The following optional arguments can be specified in any desired order: +## @end indent +## +## @item +## @var{dim}: (char, case-insensitive, can be shortened to just the first +## letter) unit of @var{range}: either "length" or "angle" (default). The +## "length" unit is supposed to be "meters" but in case of a length input, +## one can also specify any length unit accepted by validateLengthUnit. +## +## @item +## @var{angleUnit}: angle unit for all input angles: "degrees" (default) or +## "radians" (case-insensitive, can be shortened to just the first letter). +## +## @item +## @var{ellipsoid}: reference ellipsoid. Can be either an ellipsoid name, +## an ellipsoid code (entered as numerical or character string), of a vector +## of SemimajorAxis and Flattening. +## @end itemize +## +## Output arguments: +## +## @itemize +## @item +## @var{lato}, @var{lono}: computed latitude and longitude after displacement. +## For displacement inputs speciied as lengths, geodeticfwd needs to iterate +## to get a satisfactory solution. If the maximum number number of iterations +## is exceeded a suitable warning is emitted and the related output(s) is/are +## set to NaN. +## +## @item +## @var{azo}: computed azimuth at (@var{lato}, @var{lono}). +## @end itemize +## +## geodeticfwd is based on vincentyDirect.m by Alfredo Foltran, in turn based +## on Vicenty.T (1975) "Direct and inverse solutions of geodesics on the +## ellipsoid with application of nested equations". +## +## @seealso{geodeticarc, meridianfwd, reckon, referenceEllipsoid, +## validateLengthUnit, vincentyDirect} +## @end deftypefn + +function [lato, lono, azo] = geodeticfwd (varargin) + + ## Basic input checks + if (nargin < 4) + error ("geodeticfwd: too few arguments"); + elseif (nargin > 7) + error ("geodeticfwd: too many arguments"); + elseif (! (all (cellfun ("isnumeric", varargin(1:4))) && ... + cellfun ("isreal", varargin(1:4)))) + error ("geodeticfwd: all first 4 arguments should be real numeric"); + endif + + ## Check & process numeric input data sizes + isv = ! cellfun ("isscalar", varargin(1:4)); + if (any (isv)) + ## Check dimensions + cdim = cellfun (@(x) numel (size (x)), varargin(isv)); + if (any (diff (cdim))) + error ("geodeticfwd: input arrays of different dimensions"); + endif + for ii=1:cdim(1) + if (any (diff (cellfun (@(x) size (x, ii), varargin(isv))))) + error ("geodeticfwd: input arrays of different dimensions"); + endif + endfor + ## Expand any scalar args, if required + jsv = find (! isv); + for ii=1:numel (jsv) + varargin{jsv(ii)} = varargin{jsv(ii)} * ones (size (varargin{find (isv)(1)})); + endfor + endif + + lat = varargin{1}; + lon = varargin{2}; + rng = varargin{3}; + azi = varargin{4}; + varargin(1:4) = []; + + ## Check & process optional inputs. + ## First set provisional defaults + dim = "angle"; + angleUnit = "degrees"; + ## Default ellipsoid = WGS84 + ellipsoid = referenceEllipsoid (7030); + + ## Keep track of presence of extra inputs + i_length = i_angle = i_ell = false; + for ii=1:numel (varargin) + ## Only accept first mention of one extra input + if (iscellstr (varargin(ii))) + + ## Check for angle units + if (! i_angle) + if (strncmpi (varargin{ii}, "degrees", numel (varargin{ii}))) + angleUnit = "degrees"; + i_angle = true; + continue; + elseif (strncmpi (varargin{ii}, "radians", numel (varargin{ii}))) + angleUnit = "radians"; + i_angle = true; + continue; + endif + endif + + ## Check for length units (or angle) + if (! i_length) + if (strncmpi (varargin{ii}, "length", numel (varargin{ii}))) + dim = lengthUnit = "length"; + i_length = true; + continue; + elseif (strncmpi (varargin{ii}, "angle", numel (varargin{ii}))) + dim = "angle"; + i_length = true; + continue; + else + try + ## Try to convert rng to meters + rng = rng * unitsratio ("m", varargin{ii}); + ## If control gets here, the length argument was a valid length unit + dim = lengthUnit = "length"; + i_length = true; + continue; + catch + end_try_catch + endif + endif + + ## Check reference ellipsoid + if (! i_ell) + try + ellipsoid = referenceEllipsoid (varargin{ii}); + i_ell = true; + continue; + catch + end_try_catch + endif + + error ("geodeticfwd: input argument #%d not recognized", ii+4); + + elseif (cellfun ("isnumeric", varargin(ii)) && ! i_ell) + if (isscalar (varargin{ii})) + ## Assume it's a reference ellipsoid code + ellipsoid = referenceEllipsoid (varargin{ii}); + elseif (isvector (varargin{ii}) && numel (varargin{ii}) == 2) + ## Semimajoraxis and Flattening; compute Semiminoraxis + ellipsoid.SemimajorAxis = a = varargin{ii}(1); + ellipsoid.Flattening = b = varargin{ii}(2); + ecc = flat2ecc (b); + ellipsoid.SemiminorAxis = minaxis (a, ecc); + else + error ("geodeticfwd: invalid ellipsoid input (arg. #%d)", ii+4); + endif + i_ell = true; + continue; + + elseif (isstruct (varargin{ii}) && ! i_ell) + ## FIXME - check on proper ellipsoid struct + ellipsoid = varargin{ii}; + if (! all (isfield (ellipsoid, {"SemimajorAxis", "SemiminorAxis", "Flattening"}))) + error ("geodeticfwd: invalid ellipsoid struct"); + endif + continue; + + else + ## Couldn't be recognized + error ("geodeticfwd: input argument #%d not recognized", ii+4); + endif + endfor + + if (strcmpi (angleUnit, "degrees")) + ## Rest of function based on radians + lat = deg2rad (lat); + lon = deg2rad (lon); + azi = deg2rad (azi); + ## If range was given as angle, convert it as well + if (strcmpi (dim, "angle")) + rng = deg2rad (rng); + endif + endif + + ## Do the actual work. Based on Alfredo Foltran's vincentyDirect.m + ## https://github.com/alfoltranteam/octave-map/blob/master/src/vincentyDirect.m + ## (In turn based on Vicenty.T (1975) Direct and inverse solutions of + ## geodesics on the ellipsoid with application of nested equations) + + major = ellipsoid.SemimajorAxis; + minor = ellipsoid.SemiminorAxis; + f = ellipsoid.Flattening; + + iter_limit = 20; + + tanU1 = (1 - f) * tan (lat); + U1 = atan (tanU1); + sigma1 = atan2 (tanU1, cos (azi)); + cosU1 = cos (U1); + sinAlpha = cosU1 .* sin (azi); + cos2Alpha = (1 - sinAlpha) .* (1 + sinAlpha); + u2 = cos2Alpha * (major ^ 2 - minor ^ 2) / minor ^ 2; + + A = 1 + u2 / 16384 .* (4096 + u2 .* (-768 + u2 .* (320 - 175 * u2))); + B = u2 / 1024 .* (256 + u2 .* (-128 + u2 .* (74 - 47 * u2))); + + if (strcmpi (dim, "length")) + ## Be sure all lengths are in meters + if (isfield (ellipsoid, "LengthUnit") && ! isempty (ellipsoid.LengthUnit)) + rng *= unitsratio ("meter", ellipsoid.LengthUnit); + major *= unitsratio ("meter", ellipsoid.LengthUnit); + minor *= unitsratio ("meter", ellipsoid.LengthUnit); + else + ## Maybe dimensionless as in Unit Sphere, or not given. Assume "meters" + endif + sigma = rng ./ (minor * A); + lastSigma = sigma + 1; + i = 0; + ## Keep track of which sigmas have converged + sgmit = true (size (sigma)); + ## Preallocate to avoid dimension changes + lastSigma = doubleSigmaM = deltaSigma = zeros (size (sigma)); + do + i++; + lastSigma(sgmit) = sigma(sgmit); + doubleSigmaM(sgmit) = 2 * sigma1(sgmit) + sigma(sgmit); + deltaSigma(sgmit) = ... + B(sgmit) .* sin (sigma(sgmit)) .* (cos (doubleSigmaM(sgmit)) + ... + 0.25 * B(sgmit) .* (cos (sigma(sgmit)) .* (-1 + 2 * cos (doubleSigmaM(sgmit)) .^ 2) ... + - 1/6 .* B(sgmit) .* cos (doubleSigmaM(sgmit)) .* (-3 + 4 * sin (sigma(sgmit)) .^ 2) ... + .* (-3 * 4 * cos (doubleSigmaM(sgmit)) .^ 2))); + sigma(sgmit) = rng(sgmit) ./ (minor * A(sgmit)) + deltaSigma(sgmit); + ## Which inputs haven't converged yet + sgmit = abs (lastSigma - sigma) > 10e-12; + ## printf ("%2d ", i, sgmit); printf ("\n"); ## Debug, remove when tested + until (all (! sgmit) || i > iter_limit) + if (i > iter_limit) + warning ("Direct Vincenty's formulae failed to converge for some inputs!"); + sigma(sgmit) = NaN; + endif + elseif (strcmpi (dim, "angle")) + sigma = rng; + else + error ("Parameter \"dim\" must be \"angle\", \"length\" or a valid length unit!"); + endif + + doubleSigmaM = 2 * sigma1 + sigma; + sinU1 = sin (U1); + lato = atan2 (sinU1 .* cos (sigma) + cosU1 .* sin (sigma) .* cos (azi), ... + (1 - f) * sqrt (sinAlpha .^ 2 + (sinU1 .* sin (sigma) - ... + cosU1 .* cos (sigma) .* cos (azi)) .^ 2)); + lambda = atan2 (sin (sigma) .* sin (azi), ... + cosU1 .* cos (sigma) - sinU1 .* sin (sigma) .* cos(azi)); + C = f/16 * cos2Alpha .* (4 + f * (4 - 3 * cos2Alpha)); + L = lambda - (1 - C) * f .* sinAlpha .* (sigma + C .* sin (sigma) .* ... + (cos (doubleSigmaM) + C .* cos (sigma) .* (-1 + 2 * cos (doubleSigmaM) .^ 2))); + lono = L + lon; + lono = wrapToPi (lono); + + if (nargout () > 2) + azo = atan2 (sinAlpha, -sinU1 .* sin (sigma) + cosU1 .* cos (sigma) .* cos (azi)); + else + azo = []; + endif + + if (strcmpi (angleUnit, "degrees")) + lato = rad2deg (lato); + lono = rad2deg (lono); + azo = rad2deg (azo); + endif + +endfunction + + +%!error <too few arguments> geodeticfwd (1, 1, 1); +%!error <should be real numeric> geodeticfwd (1, 2, 'a', 4); +%!error <should be real numeric> geodeticfwd ({2}, 2, 'a', 4); +%!error <should be real numeric> geodeticfwd (1+2i, 2, 3, 4); +%!error <arrays of different dimensions> geodeticfwd ([1 2], [2; 1], 0, 0); +%!error <not recognized> geodeticfwd (1, 2, 3, 4, 'b'); +%!error <not recognized> geodeticfwd (1, 2, 3, 4, 'm', 'b'); +%!error <not recognized> geodeticfwd (1, 2, 3, 4, 'd', 'wgs84', 'f'); +%!error <not recognized> geodeticfwd (1, 2, 3, 4, 'd', 'wgs84', {5}); +%!error <invalid ellipsoid input> geodeticfwd (1, 2, 3, 4, [1, 2, 3]); +%!error <invalid ellipsoid struct> geodeticfwd (1, 2, 3, 4, struct ("a", "b")); +%!error <invalid ellipsoid> geodeticfwd (pi/4, -pi/4, pi/2, 10, 'nm', 'r', struct ("a", [])); + +%!test +%! [lato, lono, azo] = geodeticfwd (0, 0, pi, pi/4, 'a', 'Unit Sphere', 'r'); +%! assert (rad2deg ([lato lono azo]), [0, 180, 135], 1e-10); + +%!test +%! [lato, lono, azo] = geodeticfwd (0, 0, pi/2, pi/4, 'angle', 'Unit Sphere', 'radians'); +%! assert (rad2deg ([lato lono azo]), [45, 90, 90], 1e-10); + +%!test +%! [lato, lono, azo] = geodeticfwd (0, pi/4, pi/2, 3*pi/4, 'a', 'Unit Sphere', 'r'); +%! assert (rad2deg ([lato lono azo]), [ -45, 135, 90], 1e-10); + +%!test +%! [lato, lono, azo] = geodeticfwd (45, -45, 180, 90, 'Unit Sphere', 'd'); +%! assert ([lato lono azo], [-45, 135, 90], 1e-10); + + diff --git a/inst/geoshow.m b/inst/geoshow.m index b851a99..40249dd 100644 --- a/inst/geoshow.m +++ b/inst/geoshow.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 diff --git a/inst/gmlread.m b/inst/gmlread.m index 6f3308c..ade3b2a 100644 --- a/inst/gmlread.m +++ b/inst/gmlread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 diff --git a/inst/gpxread.m b/inst/gpxread.m index 4fbc91a..16f915e 100644 --- a/inst/gpxread.m +++ b/inst/gpxread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 @@ -142,8 +142,8 @@ function [outp] = gpxread (fname, varargin) warning ("gpxread: unknown option '%s' - ignored\n", varargin{ii}); endswitch else - error ("gpxread: wrong argument type for arg. %d - 'FeatureType' or \ -'Index' expected", ii+1); + error (["gpxread: wrong argument type for arg. %d - 'FeatureType' or ", ... + "'Index' expected"], ii+1); endif endfor diff --git a/inst/isShapeMultipart.m b/inst/isShapeMultipart.m index 3db52d3..20a3e46 100644 --- a/inst/isShapeMultipart.m +++ b/inst/isShapeMultipart.m @@ -1,4 +1,4 @@ -## Copyright (C) 2016-2020 Philip Nienhuis +## Copyright (C) 2016-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 diff --git a/inst/iso2geo.m b/inst/iso2geo.m new file mode 100644 index 0000000..572ca07 --- /dev/null +++ b/inst/iso2geo.m @@ -0,0 +1,92 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{psi} =} iso2geo (@var{phi}) +## @deftypefnx {} {@var{psi} =} iso2geo (@var{phi}, @var{spheroid}) +## @deftypefnx {} {@var{psi} =} iso2geo (@var{phi}, @var{spheroid}, @var{angleUnit}) +## Returns the isometric latitude given geodetic latitude @var{phi} +## +## Input +## @itemize +## @item +## @var{phi}: the geodetic latitude(s). Can be a scalar, vector or ND-array. +## @end itemize +## +## @itemize +## @item +## (optional) @var{spheroid}: referenceEllipsoid. For admissible values see +## `referenceEllipsoid`. The default value is 'WGS84'. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first charachter will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{psi}: the isometric latitude(s), same shape as @var{phi}. +## @end itemize +## +## Example +## @example +## iso2geo (45) +## ans = 41.170 +## @end example +## @seealso{geo2auth, geo2con, geo2iso, geo2rect} +## @end deftypefn + +function [phi] = iso2geo (psi, spheroid = "", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (psi) || ! isreal (psi)) + error ("iso2geo: numeric input expected"); + endif + + E = sph_chk(spheroid); + + if (! ischar (angleUnit)) + error ("iso2geo: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + psi = deg2rad (psi); + elseif (! strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + error ("iso2geo: illegal input for 'angleUnit'"); + endif + + ## From Snyder's "Map Projections-A Working Manual" [pg 15]. + chi = 2 * atan (exp (psi)) - pi / 2; + phi = con2geo (chi, E, "radians"); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + phi = rad2deg (phi); + endif + +endfunction + +%!test +%! psi = iso2geo (50.227465817); +%! assert (psi, 45, 1e-6) + +%!test +%! assert(iso2geo (geo2iso ([-90 : 10: 0; 0 : 10 : 90]), "wgs84"), ... +%! [-90 : 10: 0; 0 : 10 : 90], 1e-8); + +%!error <numeric> iso2geo ("s") +%!error <numeric> iso2geo (5i) + diff --git a/inst/km2deg.m b/inst/km2deg.m index 80bca2b..494078e 100644 --- a/inst/km2deg.m +++ b/inst/km2deg.m @@ -1,5 +1,5 @@ -## Copyright (C) 2008-2020 Alexander Barth <abarth93@users.sourceforge.net> -## Copyright (C) 2013-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2008-2022 Alexander Barth <abarth93@users.sourceforge.net> +## Copyright (C) 2013-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/km2nm.m b/inst/km2nm.m index 82e66b9..3e73eeb 100644 --- a/inst/km2nm.m +++ b/inst/km2nm.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/km2rad.m b/inst/km2rad.m index 3f81b35..8fbc43d 100644 --- a/inst/km2rad.m +++ b/inst/km2rad.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/km2sm.m b/inst/km2sm.m index 48bde15..d04260d 100644 --- a/inst/km2sm.m +++ b/inst/km2sm.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/kmlread.m b/inst/kmlread.m index c96f97b..c1ccd02 100644 --- a/inst/kmlread.m +++ b/inst/kmlread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 @@ -61,7 +61,7 @@ function outp = kmlread (fname) '<Placemark.*?>.*?<name>(.+?)</name>.*?<Point', "tokens"); if (! isempty (pnam)) pname = cell2mat (pnam){end}; - outp(end).Name = pnam; + outp(end).Name = pname; else outp(end).Name = ["Point #" num2str(ipt)]; endif diff --git a/inst/kmzread.m b/inst/kmzread.m index 7dcc4ad..8543f49 100644 --- a/inst/kmzread.m +++ b/inst/kmzread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 diff --git a/inst/majaxis.m b/inst/majaxis.m index a4a4c9a..b3e02f9 100644 --- a/inst/majaxis.m +++ b/inst/majaxis.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/makesymbolspec.m b/inst/makesymbolspec.m index 7e535cb..15e0cb2 100644 --- a/inst/makesymbolspec.m +++ b/inst/makesymbolspec.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Philip Nienhuis +## Copyright (C) 2015-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 diff --git a/inst/mapshow.m b/inst/mapshow.m index bf001e2..f336de2 100644 --- a/inst/mapshow.m +++ b/inst/mapshow.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 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") diff --git a/inst/meridianfwd.m b/inst/meridianfwd.m new file mode 100644 index 0000000..0d0a431 --- /dev/null +++ b/inst/meridianfwd.m @@ -0,0 +1,69 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{lat2} =} meridianfwd (@var{lat}, @var{s}) +## @deftypefnx {} {@var{lat2} =} meridianfwd (@var{lat}, @var{s}, @var{spheroid}) +## @deftypefnx {} {@var{lat2} =} meridianfwd (@var{lat}, @var{s}, @var{spheroid}, @var{angleUnit}) +## Retuns the new latitude given a starting latitude and distance travelled +## along a meridian. +## +## Inputs +## +## @itemize +## @item +## @var{lat1}: the starting latitude. +## +## @item +## @var{s} the distance travelled. The units are based on the ellipsoid. +## The default is in meters but should match that of the ellipsoid (if any). +## +## @item +## (optional) @var{spheroid}: referenceEllipsoid (parameter struct, name or +## code): the default is 'wgs84'. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first character will do). Default is 'radians'. +## +## Output +## @item +## @var{lat2}: the final latitude after travelling a distance of @var{s} +## @end itemize +## +## Example +## @example +## lat = meridianfwd (40, 1e6) +## lat = 48.983 +## @end example +## +## @seealso{meridianarc, geodeticfwd} +## @end deftypefn + +function [lat2] = meridianfwd (lat, s, spheroid = "wgs84", angleUnit = "radians") + + if (nargin < 2 || nargin > 4) + print_usage(); + endif + + lat2 = geodeticfwd (lat, 0, s, 0, spheroid, angleUnit, "length"); + +endfunction + + +%!test +%! a = rad2deg (meridianfwd (-pi/4, [60 120] * 45 * 1846.2756967249574, 'wgs84')); +%! assert (a, [0, 45], 5e-7); + diff --git a/inst/minaxis.m b/inst/minaxis.m index c2650cf..902d2e6 100644 --- a/inst/minaxis.m +++ b/inst/minaxis.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/n2ecc.m b/inst/n2ecc.m index 2bd37c8..edf06e7 100644 --- a/inst/n2ecc.m +++ b/inst/n2ecc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 diff --git a/inst/ned2aer.m b/inst/ned2aer.m index 307b572..d05a23b 100644 --- a/inst/ned2aer.m +++ b/inst/ned2aer.m @@ -1,6 +1,6 @@ -## Copyright (C) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (C) 2013-2020, Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (C) 2013-2022, Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -67,7 +67,7 @@ ## slantrange = 1000.0 ## @end example ## -## @seealso {aer2ned, ned2ecef, ned2ecefv, ned2geodetic} +## @seealso{aer2ned, ned2ecef, ned2ecefv, ned2geodetic} ## @end deftypefn ## Function adapted by anonymous contributor, see: diff --git a/inst/ned2ecef.m b/inst/ned2ecef.m index ac63daa..7a073bb 100644 --- a/inst/ned2ecef.m +++ b/inst/ned2ecef.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 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 @@ -127,18 +127,15 @@ function [x, y, z] = ned2ecef (varargin) if (! all (size (lat) == size (e)) || ... ! all (size (lon) == size (n)) || ... ! all (size (alt) == size (u))) - error ("ned2ecef: non-matching dimensions of observer points and \ -target points"); + error (["ned2ecef: non-matching dimensions of observer points and ", ... + "target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid)) - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); [x, y, z] = enu2ecef (e, n, -u, lat, lon, alt, E, angleUnit); diff --git a/inst/ned2ecefv.m b/inst/ned2ecefv.m index 7bdc137..f94e4f1 100644 --- a/inst/ned2ecefv.m +++ b/inst/ned2ecefv.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -69,7 +69,7 @@ ## w = 1000.0 ## @end example ## -## @seealso {ecef2nedv, ned2aer, ned2ecef, ned2geodetic} +## @seealso{ecef2nedv, ned2aer, ned2ecef, ned2geodetic} ## @end deftypefn ## Function adapted by anonymous contributor, see: @@ -106,11 +106,11 @@ function [u, v, w] = ned2ecefv (n, e, d, lat, lon, angleUnit = "degrees") error ("ned2ecefv: illegal input for 'angleUnit'"); endif - t = cos(lat) .* -d .- sin(lat) .* n; - w = sin(lat) .* -d .+ cos(lat) .* n; + t = cos(lat) .* -d - sin(lat) .* n; + w = sin(lat) .* -d + cos(lat) .* n; - u = cos(lon) .* t .- sin(lon) .* e; - v = sin(lon) .* t .+ cos(lon) .* e; + u = cos(lon) .* t - sin(lon) .* e; + v = sin(lon) .* t + cos(lon) .* e; endfunction diff --git a/inst/ned2geodetic.m b/inst/ned2geodetic.m index 72bb60c..7582a94 100644 --- a/inst/ned2geodetic.m +++ b/inst/ned2geodetic.m @@ -1,6 +1,6 @@ -## Copyright (c) 2014-2020 Michael Hirsch, Ph.D. -## Copyright (c) 2013-2020 Felipe Geremia Nievinski -## Copyright (C) 2020 Philip Nienhuis +## Copyright (c) 2014-2022 Michael Hirsch, Ph.D. +## Copyright (c) 2013-2022 Felipe Geremia Nievinski +## Copyright (C) 2020-2022 Philip Nienhuis ## ## Redistribution and use in source and binary forms, with or without ## modification, are permitted provided that the following conditions are met: @@ -84,7 +84,7 @@ ## alt = 1139.7 ## @end example ## -## @seealso {geodetic2ned, ned2aer, ned2ecef, ned2ecefv, referenceEllipsoid} +## @seealso{geodetic2ned, ned2aer, ned2ecef, ned2ecefv, referenceEllipsoid} ## @end deftypefn ## Function adapted by anonymous contributor, see: @@ -139,18 +139,15 @@ function [lat, lon, alt] = ned2geodetic (varargin) if (! all (size (lat0) == size (e)) || ... ! all (size (lon0) == size (n)) || ... ! all (size (alt0) == size (d))) - error ("ned2geodetic: non-matching dimensions of non-scalar observer \ -points and target points"); + error (["ned2geodetic: non-matching dimensions of non-scalar observer ", ... + "points and target points"]); endif endif - if (isempty (spheroid)) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - elseif (ischar (spheroid) || isnumeric (spheroid)) - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); [x, y, z] = ned2ecef (n, e, d, lat0, lon0, alt0, E, angleUnit); [lat, lon, alt] = ecef2geodetic (E, x, y, z, angleUnit); diff --git a/inst/nm2deg.m b/inst/nm2deg.m index 332d912..7f52f0c 100644 --- a/inst/nm2deg.m +++ b/inst/nm2deg.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/nm2km.m b/inst/nm2km.m index 596ccb2..a651885 100644 --- a/inst/nm2km.m +++ b/inst/nm2km.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/nm2rad.m b/inst/nm2rad.m index f043706..261a1c1 100644 --- a/inst/nm2rad.m +++ b/inst/nm2rad.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/nm2sm.m b/inst/nm2sm.m index 91ac8de..ed64c21 100644 --- a/inst/nm2sm.m +++ b/inst/nm2sm.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/parametricLatitude.m b/inst/parametricLatitude.m index 0ac0c2e..52697f8 100644 --- a/inst/parametricLatitude.m +++ b/inst/parametricLatitude.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 the Free Software diff --git a/inst/polycut.m b/inst/polycut.m index 4423a93..ddb8b14 100644 --- a/inst/polycut.m +++ b/inst/polycut.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 @@ -46,8 +46,8 @@ function [X, Y, Z] = polycut (varargin) if (isempty (which ("polygon2patch"))) - error ("function polygon2patch not found. OF geometry package \ -installed and loaded?"); + error (["function polygon2patch not found. OF geometry package ", ... + "installed and loaded?"]); endif [X, Y, Z] = polygon2patch (varargin{:}); diff --git a/inst/private/__dbl2int64__.m b/inst/private/__dbl2int64__.m index 0fa8478..43fd1d0 100644 --- a/inst/private/__dbl2int64__.m +++ b/inst/private/__dbl2int64__.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 diff --git a/inst/private/clippln.m b/inst/private/clippln.m index cfcf32c..c8ad73d 100644 --- a/inst/private/clippln.m +++ b/inst/private/clippln.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 diff --git a/inst/private/spheres_radius.m b/inst/private/spheres_radius.m index 10f1778..6c5ab37 100644 --- a/inst/private/spheres_radius.m +++ b/inst/private/spheres_radius.m @@ -1,4 +1,4 @@ -## Copyright (C) 2013-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2013-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/rad2km.m b/inst/rad2km.m index 3225f3d..9cb6628 100644 --- a/inst/rad2km.m +++ b/inst/rad2km.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/rad2nm.m b/inst/rad2nm.m index 06e4a88..4b91aae 100644 --- a/inst/rad2nm.m +++ b/inst/rad2nm.m @@ -1,5 +1,5 @@ -## Copyright (C) 2014-2020 Pooja Rao -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2014-2022 Pooja Rao +## Copyright (C) 2018-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 the Free Software diff --git a/inst/rad2sm.m b/inst/rad2sm.m index 917e758..ccfa050 100644 --- a/inst/rad2sm.m +++ b/inst/rad2sm.m @@ -1,5 +1,5 @@ -## Copyright (C) 2014-2020 Pooja Rao -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2014-2022 Pooja Rao +## Copyright (C) 2018-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 the Free Software diff --git a/inst/radtodeg.m b/inst/radtodeg.m index ad9b7e0..14fe6ed 100644 --- a/inst/radtodeg.m +++ b/inst/radtodeg.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/rasterclip.m b/inst/rasterclip.m index d2ea08f..fdb3f88 100644 --- a/inst/rasterclip.m +++ b/inst/rasterclip.m @@ -1,4 +1,4 @@ -## Copyright (C) 2017-2020 Philip Nienhuis +## Copyright (C) 2017-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 diff --git a/inst/rasterdraw.m b/inst/rasterdraw.m index 2b93f30..074ad8d 100644 --- a/inst/rasterdraw.m +++ b/inst/rasterdraw.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Philip Nienhuis +## Copyright (C) 2015-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 @@ -117,16 +117,17 @@ function [hr] = rasterdraw (varargin) if (! isnumeric (ibands)) error ("rasterdraw: expecting numeric value(s) for 'bands' value"); elseif (any (ibands > numel (bands))) - warning ("rasterdraw: data contains %d band(s), requested band(s) [%s\ - ] ignored.\n", numel (bands), sprintf (" %d", ibands(find (ibands > numel (bands))))); + warning (["rasterdraw: data contains %d band(s), requested ", ... + "band(s) [%s ] ignored.\n"], numel (bands), ... + sprintf (" %d", ibands(find (ibands > numel (bands))))); ibands (find (ibands > numel (bands))) = []; endif if (! (numel (ibands) == 1 || numel (ibands) == 3)) error ("rasterdraw: either one or three bands must be selected.\n"); elseif (rinfo.BitDepth >= 32) ## Assume floating point data - warning ("rasterdraw: floating point raster data. only band %d\ -will be drawn\n", ibands(1)); + warning (["rasterdraw: floating point raster data. only band ", ... + "%d will be drawn\n"], ibands(1)); endif case "colo" clrmap = varargin{ii+1}; diff --git a/inst/rasterinfo.m b/inst/rasterinfo.m index 050e48f..4f5719d 100644 --- a/inst/rasterinfo.m +++ b/inst/rasterinfo.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Philip Nienhuis +## Copyright (C) 2015-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 diff --git a/inst/rasterread.m b/inst/rasterread.m index 18b06a6..460ee90 100644 --- a/inst/rasterread.m +++ b/inst/rasterread.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Philip Nienhuis +## Copyright (C) 2015-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 diff --git a/inst/rcurve.m b/inst/rcurve.m index 31ca788..21b2e8b 100644 --- a/inst/rcurve.m +++ b/inst/rcurve.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 @@ -121,13 +121,10 @@ function r = rcurve (varargin) type = "parallel"; endif - if isempty (spheroid) - E = wgs84Ellipsoid; - elseif (isstruct (spheroid)) - E = spheroid; - else - E = referenceEllipsoid (spheroid); + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); endif + E = sph_chk (spheroid); if (! ischar (angleUnit) || ! ismember (lower (angleUnit(1)), {"d", "r"})) error ("rcurve: angleUnit should be one of 'degrees' or 'radians'") diff --git a/inst/reckon.m b/inst/reckon.m index 090ec57..c8f7ed7 100644 --- a/inst/reckon.m +++ b/inst/reckon.m @@ -1,4 +1,4 @@ -## Copyright (C) 2008-2020 Alexander Barth <abarth93@users.sourceforge.net> +## Copyright (C) 2008-2022 Alexander Barth <abarth93@users.sourceforge.net> ## ## 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 @@ -25,6 +25,8 @@ ## ## This function can also be used to define a spherical coordinate system ## with rotated poles. +## +## @seealso{geodeticfwd, vincentyDirect} ## @end deftypefn ## Author: Alexander Barth <barth.alexander@gmail.com> diff --git a/inst/rect2geo.m b/inst/rect2geo.m new file mode 100644 index 0000000..85d2cda --- /dev/null +++ b/inst/rect2geo.m @@ -0,0 +1,124 @@ +## Copyright (C) 2022 The Octave Project Developers +## +## 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 {} {@var{phi} =} rect2geo (@var{mu}) +## @deftypefnx {} {@var{phi} =} rect2geo (@var{mu}, @var{spheroid}) +## @deftypefnx {} {@var{phi} =} rect2geo (@var{mu}, @var{spheroid}, @var{angleUnit}) +## Returns the geodetic latitude given rectifying latitude @var{mu} +## +## Input +## @itemize +## @item +## @var{mu}: the rectifying latitude(s). Can be a scalar value, a vector or +## an nD array. +## +## @item +## (optional) @var{spheroid}: referenceEllipsoid parameter struct: the +## default is wgs84. +## +## @item +## (optional) @var{angleUnit}: string for angular units ('degrees' or 'radians', +## case-insensitive, just the first character will do). Default is 'degrees'. +## @end itemize +## +## Output +## @itemize +## @item +## @var{phi}: the geodetic latitude(s), same shape as @var{mu}. +## @end itemize +## +## Example +## @example +## rect2geo(44.856) +## ans = 45.000 +## @end example +## @seealso{geo2auth, geo2con, geo2iso, geo2rect} +## @end deftypefn + +function [phi] = rect2geo (mu, spheroid ="", angleUnit = "degrees") + + if (nargin < 1 || nargin > 3) + print_usage(); + endif + + if (! isnumeric (mu) || ! isreal (mu)) + error ("rect2geo: numeric input expected"); + endif + + if (isnumeric (spheroid)) + spheroid = num2str (spheroid); + endif + E = sph_chk (spheroid); + + if (! ischar (angleUnit)) + error ("rect2geo: character value expected for 'angleUnit'"); + elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) + ## Latitude must be within [-90 ... 90] + if (any (abs (mu) > 90)) + error ("rect2geo: azimuth value(s) out of acceptable range [-90, 90]") + endif + mu = deg2rad (mu); + elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) + ## Latitude must be within [-pi/2 ... pi/2] + if (any (abs (mu) > pi / 2)) + error("rect2geo: azimuth value(s) out of acceptable range (-pi/2, pi/2)") + endif + else + error ("rect2geo: illegal input for 'angleUnit'"); + endif + + if (isfield (E, "ThirdFlattening") == 1) + n = E.ThirdFlattening; + else + ecc = E.Eccentricity; + ## From Snyder's "Map Projections-A Working Manual" [pg 17]. + ecc1 = (1 - ecc ^ 2) ^ ( 1 / 2); + n = (1 - ecc1) / (1 + ecc1); + endif + + ## From R.E. Deakin "A FRESH LOOK AT THE UTM PROJECTION" [pg 5] + n2 = n ^ 2; + n3 = n ^ 3; + n4 = n ^ 4; + phi = mu + ... + (((3 / 2) * n - (27 / 32) * n3) * sin (2 * mu)) + ... + (((21 / 16) * n2 - (55 / 32) * n4) * sin (4 * mu)) + ... + (((151 / 48) * n3) * sin (6 * mu)) + ... + (((1097 / 512) * n4) * sin (8 * mu)); + + if (strncmpi (angleUnit, "degrees", length (angleUnit))) + phi = rad2deg (phi); + endif + +endfunction + + +%!test +%! mu = [0:15:90]; +%! phi = rect2geo (mu); +%! Z = degrees2dm (phi - mu); +%! check = [0, +%! 4.3406136, +%! 7.5100085, +%! 8.6590367, +%! 7.4879718, +%! 4.3185768, +%! 0]; +%! assert (Z(:,2), check, 1e-6) + +%!error <numeric> rect2geo ("s") +%!error <numeric> rect2geo (5i) + diff --git a/inst/referenceEllipsoid.m b/inst/referenceEllipsoid.m index 3920a60..229271b 100644 --- a/inst/referenceEllipsoid.m +++ b/inst/referenceEllipsoid.m @@ -1,5 +1,5 @@ -## Copyright (C) 2018-2020 Philip Nienhuis -## Copyright (C) 2014 Alfredo Foltran <alfoltran@gmail.com> +## Copyright (C) 2018-2022 Philip Nienhuis +## Copyright (C) 2014-2022 Alfredo Foltran <alfoltran@gmail.com> ## ## 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 @@ -105,7 +105,7 @@ ## Function supplied by anonymous contributor, see: ## https://savannah.gnu.org/patch/index.php?9634 -function Ell = referenceEllipsoid (name="wgs84", unit="meter") +function Ell = referenceEllipsoid (name="unit sphere", unit="meter") ## List of implemented codes. To be updated if codes are added. persistent codes; @@ -113,14 +113,14 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") codes = {"7001", "Airy 1830", "Airy30"; ... "7002", "Airy Modified 1849", "Airy49"; ... "7003", "Australian National Spheroid", " "; ... - "7004", "Bessel 1841", "Bessel "; ... + "7004", "Bessel 1841", "Bessel"; ... "7005", "Bessel Modified", " "; ... "7006", "Bessel Namibia", " "; ... "7007", "Clarke 1858", " "; ... - "7008", "Clarke 1866", " "; ... + "7008", "Clarke 1866", "Clarke66"; ... "7009", "Clarke 1866 Michigan", " "; ... " ", "Clarke 1878" " "; ... - "7010", "Clarke 1880 (Benoit)", " "; ... + "7010", "Clarke 1880 (Benoit)", "Clarke80"; ... "7011", "Clarke 1880 (IGN)", " "; ... "7012", "Clarke 1880 (RGS)", " "; ... "7013", "Clarke 1880 (Arc)", " "; ... @@ -173,13 +173,14 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") " ", "Saturn", " "; ... " ", "Uranus", " "; ... " ", "Neptune", " "; ... - " ", "Pluto", " "}; + " ", "Pluto", " "; ... + " ", "Unit Sphere", " "}; endif if (isnumeric (name) && isreal (name)) name = num2str (fix (name)); elseif (! ischar (name)) - error ("referenceEllipsoid: value must be a string or real number"); + error ("referenceEllipsoid: value must be a string or integer number"); elseif (strcmp (name, " ")) error ("referenceEllipsoid: name required"); endif @@ -252,7 +253,7 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") InverseFlattening = 294.978697164674; case lower (codes(10, :)) - Code = " "; + Code = []; Name = "Clarke 1878"; SemimajorAxis = 6378190; InverseFlattening = 293.4659980; @@ -312,7 +313,7 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") InverseFlattening = 298.257222101; case lower (codes(20, :)) - Code = " "; + Code = []; Name = "Hayford"; SemimajorAxis = 6378388; InverseFlattening = 297; @@ -336,7 +337,7 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") InverseFlattening = 297; case lower (codes(24, :)) - Code = " "; + Code = []; Name = "New International 1967"; SemimajorAxis = 6378157.5; InverseFlattening = 298.24961539; @@ -348,7 +349,7 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") InverseFlattening = 298.3; case lower (codes(26, :)) - Code = " "; + Code = []; Name = "Maupertius"; SemimajorAxis = 6397300; InverseFlattening = 191; @@ -428,11 +429,12 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") case lower (codes(39, :)) Code = 7042; Name = "Everest (1830 Definition)"; - SemimajorAxis = 20922931.8; + ## SemimajorAxis = 20922931.8; # Indian feet (= 0.99999566 British foot) + SemimajorAxis = 6377281.935116282; InverseFlattening = 300.8017; case lower (codes(40, :)) - Code = " "; + Code = []; Name = "World Geodetic System 1966"; SemimajorAxis = 6378145; InverseFlattening = 298.25; @@ -504,93 +506,100 @@ function Ell = referenceEllipsoid (name="wgs84", unit="meter") InverseFlattening = 297; case lower (codes(52, :)) - Code = " "; + Code = []; Name = "IERS 1989"; SemimajorAxis = 6378136; InverseFlattening = 298.257; case lower (codes(53, :)) - Code = " "; + Code = []; Name = "IERS 2003"; SemimajorAxis = 6378136.6; InverseFlattening = 298.25642; case lower (codes(54, :)) - Code = " "; + Code = []; Name = "Sun"; SemimajorAxis = 695700000; InverseFlattening = 111111; case lower (codes(55, :)) - Code = " "; + Code = []; Name = "Mercury"; SemimajorAxis = 2440530; InverseFlattening = 1075; case lower (codes(56, :)) - Code = " "; + Code = []; Name = "Venus"; SemimajorAxis = 6051800; InverseFlattening = Inf; case lower (codes(57, :)) - Code = " "; + Code = []; Name = "Earth"; SemimajorAxis = 6378137; InverseFlattening = 298.2572235630; case lower (codes(58, :)) - Code = " "; + Code = []; Name = "Moon"; SemimajorAxis = 1738100; InverseFlattening = 833.33; case lower (codes(59, :)) - Code = " "; + Code = []; Name = "Mars"; SemimajorAxis = 3396190; InverseFlattening = 169.894; case lower (codes(60, :)) - Code = " "; + Code = []; Name = "Jupiter"; SemimajorAxis = 71492000; InverseFlattening = 15.4144; case lower (codes(61, :)) - Code = " "; + Code = []; Name = "Saturn"; SemimajorAxis = 60268000; InverseFlattening = 10.208; case lower (codes(62, :)) - Code = " "; + Code = []; Name = "Uranus"; SemimajorAxis = 25559000; InverseFlattening = 43.616; case lower (codes(63, :)) - Code = " "; + Code = []; Name = "Neptune"; SemimajorAxis = 24764000; InverseFlattening = 58.5437; case lower (codes(64, :)) - Code = " "; + Code = []; Name = "Pluto"; SemimajorAxis = 1188300; InverseFlattening = Inf; + case lower (codes(65, :)) + Code = []; + Name = "Unit Sphere"; + SemimajorAxis = 1; + InverseFlattening = Inf; + case "0" - ## Show list of codes - printf ("\n referenceEllipsoid.m:\n list of implemented ellipsoids:\n"); - printf (" Code Alias 1 Alias 2\n"); - printf (" ==== ======= =======\n"); - for ii=1:size (codes, 1) - printf ("%5s %20s %15s\n", codes (ii, :){:}); - endfor if (nargout > 0) Ell = [(num2cell (1:size (codes, 1))') codes] ; + else + ## Show list of codes + printf ("\n referenceEllipsoid.m:\n list of implemented ellipsoids:\n"); + printf (" Code Alias 1 Alias 2\n"); + printf (" ==== ======= =======\n"); + for ii=1:size (codes, 1) + printf ("%5s %20s %15s\n", codes (ii, :){:}); + endfor endif return diff --git a/inst/referenceSphere.m b/inst/referenceSphere.m new file mode 100644 index 0000000..702c7af --- /dev/null +++ b/inst/referenceSphere.m @@ -0,0 +1,196 @@ +## Copyright (C) 2022 John W. Eaton +## +## 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 PURPOSEll. 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; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} referenceSphere (@var{name}, @var{unit}) +## Return the parameters of the named spherical object. +## +## Valid names are "Unit Sphere", "Sun", "Sun", "Mercury", "Venus", +## "Earth", "Moon", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", +## and "Pluto". Case is not important. +## +## @var{unit} can be the name of any unit accepted by function +## validateLengthUnit.m. Also here case is not important. +## +## The output consists of a scalar struct with fields "Name", "LengthUnit", +## "Radius", "SemimajorAxis", "SemiminorAxis", "InverseFlattening", +## "Eccentricity", "Flattening", "ThirdFlattening", "MeanRadius", +## "SurfaceArea", and "Volume". +## +## Examples: +## +## @example +## >> E = referenceSphere () +## E = +## +## scalar structure containing the fields: +## +## Name = Unit Sphere +## LengthUnit = +## Radius = 1 +## SemimajorAxis = 1 +## SemiminorAxis = 1 +## InverseFlattening = Inf +## Eccentricity = 0 +## Flattening = 0 +## ThirdFlattening = 0 +## MeanRadius = 1 +## SurfaceArea = 12.566 +## Volume = 4.1888 +## @end example +## +## @example +## >> E = referenceSphere ("Earth", "km") +## E = +## +## scalar structure containing the fields: +## +## Name = Earth +## LengthUnit = km +## Radius = 6371 +## SemimajorAxis = 6371 +## SemiminorAxis = 6371 +## InverseFlattening = Inf +## Eccentricity = 0 +## Flattening = 0 +## ThirdFlattening = 0 +## MeanRadius = 6371.0 +## SurfaceArea = 5.1006e+08 +## Volume = 1.0832e+12 +## @end example +## +## @seealso{validateLengthUnit, referenceEllipsiod} +## @end deftypefn + +## Parts of this function were adapted from referenceEllipsoid. + +function Sph = referenceSphere (name = "unit", unit = "") + + if (! ischar (name)) + error ("referenceSphere: NAME must be a string"); + endif + + if (! ischar (unit)) + error ("referenceSphere: units name expected for input arg. #2"); + endif + + switch (lower (name)) + case "sun" + Name = "Sun"; + Radius = 694460000; + + case "mercury" + Name = "Mercury"; + Radius = 2439000; + + case "venus" + Name = "Venus"; + Radius = 6051000; + + case "earth" + Name = "Earth"; + Radius = 6371000; + + case "moon" + Name = "Moon"; + Radius = 1738000; + + case "mars" + Name = "Mars"; + Radius = 3390000; + + case "jupiter" + Name = "Jupiter"; + Radius = 69882000; + + case "saturn" + Name = "Saturn"; + Radius = 58235000; + + case "uranus" + Name = "Uranus"; + Radius = 25362000; + + case "neptune" + Name = "Neptune"; + Radius = 24622000; + + case "pluto" + Name = "Pluto"; + Radius = 1151000; + + case "unit" + Name = "Unit Sphere"; + Radius = 1.0; + + otherwise + error ("referenceSphere: unrecognized sphere: %s", name); + + endswitch + + if (! strcmpi (name, "unit") && isempty (unit)) + if (nargin == 2) + error ("referenceSphere: UNIT may not be empty for %s", Name); + else + unit = "Meters"; + endif + endif + + Sph = param_calc (Name, Radius, unit); + +endfunction + + +function sph = param_calc (Name, Radius, unit) + + if (! isempty (unit)) + ratio = unitsratio (unit, "Meters"); + Radius = Radius * ratio; + endif + + sph.Name = Name; + sph.LengthUnit = unit; + sph.Radius = Radius; + sph.SemimajorAxis = Radius; + sph.SemiminorAxis = Radius; + sph.InverseFlattening = Inf; + sph.Eccentricity = 0; + sph.Flattening = 0; + sph.ThirdFlattening = 0; + sph.MeanRadius = Radius; + sph.SurfaceArea = 4 * pi * Radius^2; + sph.Volume = (4 * pi) / 3 * Radius^3; + +endfunction + + +%!test +%! U = referenceSphere ("Unit"); +%! assert (U.LengthUnit, ""); +%! assert (U.Radius, 1); +%! assert (U.SemimajorAxis, 1); +%! assert (U.SemiminorAxis, 1); +%! assert (U.InverseFlattening, Inf); +%! assert (U.Eccentricity, 0); +%! assert (U.Flattening, 0); +%! assert (U.ThirdFlattening, 0); +%! assert (U.MeanRadius, 1); +%! assert (U.SurfaceArea, 4*pi); +%! assert (U.Volume, 4*pi/3); + +%!error <NAME must be a string> referenceSphere (7i) +%!error <unrecognized sphere> referenceSphere ("yy") +%!error <unknown unit> referenceSphere ("Unit", "##not a unit@@") diff --git a/inst/removeExtraNanSeparators.m b/inst/removeExtraNanSeparators.m index 618af11..02108a9 100644 --- a/inst/removeExtraNanSeparators.m +++ b/inst/removeExtraNanSeparators.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/roundn.m b/inst/roundn.m index a8a7bec..923e54f 100644 --- a/inst/roundn.m +++ b/inst/roundn.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Markus Bergholz <markuman@gmail.com> +## Copyright (C) 2015-2022 Markus Bergholz <markuman@gmail.com> ## ## 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 diff --git a/inst/scxsc.m b/inst/scxsc.m index cb29f9c..ad713be 100644 --- a/inst/scxsc.m +++ b/inst/scxsc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 The Octave Project Developers +## Copyright (C) 2022 The Octave Project Developers ## ## 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 @@ -179,40 +179,40 @@ function [lat, lon, idl] = scxsc (varargin) ## Explore spherical distances between circle centers. ## Use haversine formula for approx.distances < 25..50 degrees - id = abs ((vect(:, 1) + 360) .- (vect(:, 4) + 360)) > 25; + id = abs ((vect(:, 1) + 360) - (vect(:, 4) + 360)) > 25; id = id | abs ((vect(:, 2) + 180) - (vect(:, 5) + 180)) > 25; ## Plain spherical cosine formula sphd = NaN (nv, 1); - sphd(id) = acosd (sind (vect(id, 1)) .* sind (vect(id, 4)) .+ ... + sphd(id) = acosd (sind (vect(id, 1)) .* sind (vect(id, 4)) + ... cosd (vect(id, 1)) .* cosd (vect(id, 4)) .* ... - cosd (vect(id, 2) .- vect(id, 5))); + cosd (vect(id, 2) - vect(id, 5))); ## Haversine formula sphd(! id) = 2 * asind (sqrt ( ... - (sind (abs (vect(! id, 1) .- vect(! id, 4)) / 2)) .^ 2 .+ ... + (sind (abs (vect(! id, 1) - vect(! id, 4)) / 2)) .^ 2 + ... cosd (vect(! id, 1)) .* cosd (vect(! id, 4)) .* ... - sind (abs (vect(! id, 2) .- vect(! id, 5)) / 2) .^2)); + sind (abs (vect(! id, 2) - vect(! id, 5)) / 2) .^2)); ## Init tracker for various reasons for issues with intersections idl = zeros (nv, 1); ## Find circle pairs that have no intersections. Case 1: sphd > sum of radii - idl(sphd .- (vect(:, 3) .+ vect(:, 6)) > 0) = 1; + idl(sphd - (vect(:, 3) + vect(:, 6)) > 0) = 1; ## Find circle pairs that have no intersections. Case 2: one circle ## completely lying in the other w/o tangent points - idl(vect(:, 3) .- sphd > vect(:, 6)) = 1; - idl(vect(:, 6) .- sphd > vect(:, 3)) = 1; + idl(vect(:, 3) - sphd > vect(:, 6)) = 1; + idl(vect(:, 6) - sphd > vect(:, 3)) = 1; ## Find circle pairs with coinciding (or antipodal) centers ==> ## Some may have no intersections (one enclosed in the other) [alat, alon] = antipode (vect(:, 4), vect(:, 5)); - idlc = find ((abs (vect(:, 1) .- vect (:, 4)) < 1e-13 & ... - abs (vect(:, 2) .- vect (:, 5)) < 1e-13) | ... - (abs (vect(:, 1) .- alat) < 1e-13 & ... - abs (vect(:, 2) .- alon) < 1e-13)); + idlc = find ((abs (vect(:, 1) - vect (:, 4)) < 1e-13 & ... + abs (vect(:, 2) - vect (:, 5)) < 1e-13) | ... + (abs (vect(:, 1) - alat) < 1e-13 & ... + abs (vect(:, 2) - alon) < 1e-13)); idl(idlc) = 1; ## Some of these may have coinciding circles - idlca = abs (vect(idlc, 3) .- vect(idlc, 6)) < 2 * eps | ... - abs (vect(idlc, 3) - 180 .+ vect(idlc, 6)) < 2 * eps; + idlca = abs (vect(idlc, 3) - vect(idlc, 6)) < 2 * eps | ... + abs (vect(idlc, 3) - 180 + vect(idlc, 6)) < 2 * eps; idl(idlc(idlca)) = 3; ## Set up pointer to valid pairs idv = idl < 1; @@ -255,17 +255,17 @@ function [lat, lon, idl] = scxsc (varargin) r1 = vect(idv, 3); r2 = vect(idv, 6); - a = (cosd (r1) .- q .* cosd (r2)) ./ (1 .- q2); - b = (cosd (r2) .- q .* cosd (r1)) ./ (1 .- q2); + a = (cosd (r1) - q .* cosd (r2)) ./ (1 - q2); + b = (cosd (r2) - q .* cosd (r1)) ./ (1 - q2); - x0 = a .* c1 .+ b .* c2; + x0 = a .* c1 + b .* c2; x02 = dot (x0, x0, 2); ## Dot product cannot be larger than 1 here x02 = sign (x02) .* min (1, abs (x02)); - t = sqrt ((1 .- x02) ./ n2); # <=== No chance n2 can be zero? (div by zero !!!) - p1 = x0 .+ (t .* n); - p2 = x0 .- (t .* n); + t = sqrt ((1 - x02) ./ n2); # <=== No chance n2 can be zero? (div by zero !!!) + p1 = x0 + (t .* n); + p2 = x0 - (t .* n); lat(idv, 1) = atan2d (p1(:, 3), hypot (p1(:, 1), p1(:, 2))); lon(idv, 1) = atan2d (p1(:, 2), p1(:, 1)); diff --git a/inst/shapedraw.m b/inst/shapedraw.m index e97a080..b7b3830 100644 --- a/inst/shapedraw.m +++ b/inst/shapedraw.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 @@ -201,6 +201,8 @@ function [h] = shapedraw (shpf, varargin) [tmp(2:2:2*numel(shpf)-1)] = NaN; Y = [tmp{:}]; endif + ## Infer max. nr. of vertices for patches + mxverts = max (cellfun ("numel", tmp)) + 1; ## For "extended" (not strictly ML-compatible) mapstructs: if (isfield (shpf, "Z")) ## Give it a try. Could be one Z per shape, rather than one per vertex @@ -454,7 +456,7 @@ function [h] = shapedraw (shpf, varargin) if (isempty (h_idx)) mp(ii) = []; else - hdx(mp) = h_idx; + hdx(mp(ii)) = h_idx; endif endfor endif @@ -465,7 +467,7 @@ function [h] = shapedraw (shpf, varargin) idx = [ jdx (numel(X)+2) ]; endif ## Prepare 'faces' argument for patch - faces = NaN (numel(idx)-1, max (diff (find (isnan ([NaN; X; NaN]))))); + faces = NaN (numel(idx)-1, mxverts); for ii=1:numel (idx) - 1 XT = X(idx(ii):idx(ii+1)-2); YT = Y(idx(ii):idx(ii+1)-2); diff --git a/inst/shapeinfo.m b/inst/shapeinfo.m index e7416e6..34cac1d 100644 --- a/inst/shapeinfo.m +++ b/inst/shapeinfo.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 diff --git a/inst/shaperead.m b/inst/shaperead.m index efce946..6fcd358 100644 --- a/inst/shaperead.m +++ b/inst/shaperead.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 @@ -102,11 +102,16 @@ ## and returned in the output struct(s). To limit this to just some ## attributes, enter a value consisting of a cell array of attribute names to ## be read. To have no attributes read at all, specify @{@}, an empty cell -## array. +## array. @* +## Attributes "Geometry", "BoundingBox", "X", "Y", "Lat", and "Lon" are +## contained within the .shp file and will be returned in any case. If +## specified these values will be ignored. In the unlikely case that the +## associated .dbf file also contains these attributes (truncated to 10 +## characters) they will be prepended with an underscore (@code {_}). ## ## @item BoundingBox ## Select only those shape items (features) whose bounding box lies within, or -## intersets in at least one point with the limits of the BoundingBox value (a +## intersects in at least one point with the limits of the BoundingBox value (a ## 2 X 2 double array [Minx, MinY; MaxX, MaxY]). ## No intersection or clipping with the BoundingBox value will be done by ## default! @@ -221,8 +226,8 @@ function [ outs, oatt ] = shaperead (fname, varargin); elseif (ischar (outopts)) outopts = lower (outopts); if (! any (strncmp (outopts, {"ml", "ext", "oct", "dat"}, 1))) - error ("shaperead: arg. #2 char value should be one of 'ml', 'ext', \ -'oct' or 'dat'\n"); + error (["shaperead: arg. #2 char value should be one of 'ml', 'ext', ", ... + "'oct' or 'dat'\n"]); endif switch outopts case {"ml", "m"} @@ -234,9 +239,8 @@ function [ outs, oatt ] = shaperead (fname, varargin); case {"dat", "d"} outopts = 3; otherwise - error ("shaperead: illegal value for arg. #2: '%s' - expected 'ml', \ -'ext', 'oct' of 'dat'", ... - outopts); + error (["shaperead: illegal value for arg. #2: '%s' - expected ", ... + "'ml', 'ext', 'oct' of 'dat'"], outopts); endswitch else error ("shaperead: numeric or character type expected for arg. #2\n"); @@ -290,6 +294,9 @@ function [ outs, oatt ] = shaperead (fname, varargin); case "att" ## Select records based on attribute values s_atts = varargin{ii+1}; + ## Weed out .shp file attributes + s_atts = setdiff (s_atts, {"Geometry", "BoundingBox", ... + "X", "Y", "Lat", "Lon"}); case "bou" ## Select whether record/shape features partly lie inside or outside limits try @@ -322,8 +329,8 @@ function [ outs, oatt ] = shaperead (fname, varargin); try s_recs = sort (unique (double (varargin{ii+1}(:)))); if (have_shx && any (s_recs > nrec)) - printf ("shaperead: requ. record nos. > nr. of records (%d) ignored\n", ... - nrec); + printf (["shaperead: requested record nos. > nr. of records ", ... + "(%d) ignored\n"], nrec); s_recs (find (srecs > nrec)) = []; endif catch @@ -613,8 +620,9 @@ function [ outs, oatt ] = shaperead (fname, varargin); nullsh = [ nullsh ir ]; rincl = 0; elseif (abs(val(1, 6)) > 31) - error ("shaperead: unexpected shapetype value %f for feature # %d\n\ - Looks like a faulty shape file.", val(1, 6), ir); + error (["shaperead: unexpected shapetype value %f for feature ", ... + "# %d\n Looks like a faulty shape file."], ... + val(1, 6), ir); endif endswitch @@ -736,8 +744,8 @@ function [ outs, oatt ] = shaperead (fname, varargin); endif else if (! unsupp) - warning ("shaperead: shapefile contains unsupported shape \ -types\n"); + warning (["shaperead: shapefile contains unsupported ", ... + "shape types\n"]); outs = oatt = []; return endif @@ -908,44 +916,52 @@ types\n"); ## Try to read the .dbf file try atts = dbfread ([ bname ".dbf" ], s_recs, s_atts); - if (outopts < 2) - ## First check if any attributes match fieldnames; if so, pre-/append "_" - tags = {"Geometry", "BoundingBox", "X", "Y", "Lat", "Lon"}; - for ii=1:numel (tags) - idx = find (strcmp (tags{ii}, atts(1, :))); - if (! isempty (idx)) - atts(1, idx) = ["_" atts{1, idx} "_"]; - endif - endfor - ## Matlab style map-/geostruct. Divide attribute values over struct elems - if (nargout < 2) - ## Attributes appended to outs struct - for ii=1:size (atts, 2) - [outs.(atts{1, ii})] = deal (atts(2:end, ii){:}); + if (! isempty (atts)) + if (outopts < 2) + ## First check if any attributes match fieldnames; if so, pre-/append "_" + tags = {"Geometry", "BoundingBox", "X", "Y", "Lat", "Lon"}; + for ii=1:numel (tags) + idx = find (strcmp (tags{ii}, atts(1, :))); + if (! isempty (idx)) + atts(1, idx) = ["_" atts{1, idx} "_"]; + endif endfor - oatt = []; + ## Matlab style map-/geostruct. Divide attribute values over struct elems + if (nargout < 2) + ## Attributes appended to outs struct + for ii=1:size (atts, 2) + [outs.(atts{1, ii})] = deal (atts(2:end, ii){:}); + endfor + oatt = []; + else + ## Attributes separately in oatt struct + oatt(size (atts, 1) - 1).(atts{1, 1}) = []; + for ii=1:size (atts, 2) + [oatt.(atts{1, ii})] = deal (atts(2:end, ii){:}); + endfor + oatt = oatt'; + endif else - ## Attributes separately in oatt struct - oatt(size (atts, 1) - 1).(atts{1, 1}) = []; + ## Octave output struct. Add attributes columns as struct fields + ## Check if any attributes match fieldnames; if so, pre-/append "_" + tags = {"shpbox", "vals", "bbox", "npt", "npr", "idx", "Geometry", ... + "BoundingBox", "X", "Y", "Lat", "Lon"}; + for ii=1:numel (tags) + idx = find (strcmp (tags{ii}, atts(1, :))); + if (! isempty (idx)) + atts(1, idx) = ["_" atts{1, idx} "_"]; + endif + endfor for ii=1:size (atts, 2) - [oatt.(atts{1, ii})] = deal (atts(2:end, ii){:}); + outs.fields(end+1) = atts{1, ii}; + if (islogical (atts{2, ii}) || isnumeric (atts{2, ii})) + outs = setfield (outs, atts{1, ii}, cell2mat (atts(2:end, ii))); + else + outs = setfield (outs, atts{1, ii}, atts(2:end, ii)); + endif endfor - oatt = oatt'; + atts = []; endif - else - ## Octave output struct. Add attributes columns as struct fields - ## Check if any attributes match fieldnames; if so, pre-/append "_" - tags = {"shpbox", "vals", "bbox", "npt", "npr", "idx", "Geometry", ... - "BoundingBox", "X", "Y", "Lat", "Lon"}; - for ii=1:size (atts, 2) - outs.fields(end+1) = atts{1, ii}; - if (islogical (atts{2, ii}) || isnumeric (atts{2, ii})) - outs = setfield (outs, atts{1, ii}, cell2mat (atts(2:end, ii))); - else - outs = setfield (outs, atts{1, ii}, atts(2:end, ii)); - endif - endfor - atts = []; endif catch printf ("shaperead: file %s couldn't be read;\n", [ bname ".dbf" ]); diff --git a/inst/shapewrite.m b/inst/shapewrite.m index 02b1574..e673b6b 100644 --- a/inst/shapewrite.m +++ b/inst/shapewrite.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Philip Nienhuis +## Copyright (C) 2014-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 @@ -103,8 +103,8 @@ function [status] = shapewrite (shp, fname="", atts=[]) ## Check for dbfwrite function if (isempty (which ("dbfwrite"))) - error ("shapewrite: dbfwrite function not found. (io package installed \ -and loaded?)"); + error (["shapewrite: dbfwrite function not found. (io package ", ... + "installed and loaded?)"]); return; endif @@ -128,8 +128,8 @@ and loaded?)"); if (nargin > 2) if (isstruct (atts)) if (! n_dbfspec) - warning ("shapewrite: DbfSpec not implemented; including requested \ -attributes\n"); + warning (["shapewrite: DbfSpec not implemented; including ", ... + "requested attributes\n"]); n_dbfspec = 1; endif ## Get attribute names from field "FieldName"; allow lowercase and camelcase @@ -139,11 +139,12 @@ attributes\n"); ## Get fieldnames out of first match atts = {atts.(fieldnames (atts (fnidx)){1})}; else - warning ("shapewrite: no field 'fieldname' (case-insensitive) found \ -in struct\n=> input arg. #3 ignored"); + warning (["shapewrite: no field 'fieldname' (case-insensitive) ", ... + "found in struct\n=> input arg. #3 ignored"]); endif elseif (! iscellstr (atts) && ! ischar (atts)) - error ("shapewrite: arg.#3: attribute name or cellstr array of attribute names expected"); + error (["shapewrite: arg.#3: attribute name or cellstr array of ", ... + "attribute names expected"]); endif ## Check if requested attributes exist at all in shapestruct atts = unique (atts); @@ -257,7 +258,7 @@ in struct\n=> input arg. #3 ignored"); ## Nr. of vertices npt = numel (shp(ishp).X); ## Pointers to parts - nptr = idx(1:end-1) .- [0:numel(idx)-2]; + nptr = idx(1:end-1) - [0:numel(idx)-2]; ## Write record contents switch (stype) diff --git a/inst/sm2deg.m b/inst/sm2deg.m index 829dc65..cec5c34 100644 --- a/inst/sm2deg.m +++ b/inst/sm2deg.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/sm2km.m b/inst/sm2km.m index 6555536..9f1774d 100644 --- a/inst/sm2km.m +++ b/inst/sm2km.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/sm2nm.m b/inst/sm2nm.m index ed0db61..0f52d79 100644 --- a/inst/sm2nm.m +++ b/inst/sm2nm.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> +## Copyright (C) 2014-2022 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> ## ## 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 diff --git a/inst/sm2rad.m b/inst/sm2rad.m index e30525e..2b85da5 100644 --- a/inst/sm2rad.m +++ b/inst/sm2rad.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Pooja Rao <poojarao12@gmail.com> +## Copyright (C) 2014-2022 Pooja Rao <poojarao12@gmail.com> ## ## 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 diff --git a/inst/sph_chk.m b/inst/sph_chk.m new file mode 100644 index 0000000..dac177d --- /dev/null +++ b/inst/sph_chk.m @@ -0,0 +1,73 @@ +## Copyright (C) 2021-2022 The Octave Project Developers +## +## 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 {} {@var{ref_ell} =} sph_chk (@var{spheroid}) +## @deftypefnx {} {@var{ref_ell} =} sph_chk (@var{spheroid}, @var{req_fields}) +## Check validity of iput ellipsoids / spheroids. +## +## Inputs: +## @itemize +## @item +## @var{spheroid}: input ellipsoid or spheroid. Can be an ellipsoid or +## spheroid name, or ditto code (character string or scalar numeric value), or +## a numeric 2-element vector containing SemimajorAxis and Eccentricity +## values, or an ellipsoid / spheroid struct. If omitted a WGS84 ellipsoid +## struct is returned. +## +## @item +## @var{req_fields}: cellstr array of required fields of input ellipsoid / +## spheroid struct; only used/useful when an ellipsoid / spheroid struct +## value is given. The default required fields are "SemimajorAxis", +## "Flattening" and "LengthUnit". @* +## Note: the validity of the fields input this way isn't checked! +## @end itemize +## +## Output: @* +## Ellipsoid / spheroid struct. +## +## @seealso{} +## @end deftypefn + +function [E] = sph_chk (spheroid, req_flds = {"SemimajorAxis", "Flattening", "LengthUnit"}) + + if (isempty (spheroid)) + E = wgs84Ellipsoid; + elseif (isstruct (spheroid)) + E = spheroid; + ## Check fields + flds = isfield (E, req_flds); + if (! all (flds)) + error ("Vital fields missing from ellipsoid: %s", ... + sprintf ("%s ", req_flds(! flds))); + endif + elseif (isnumeric (spheroid) && isreal (spheroid) && isvector (spheroid)) + if (numel (spheroid) != 2) + error("sph_chk: 2-element vector expected"); + elseif (spheroid(1) < 1) + E.Eccentricity = spheroid(1); + E.SemimajorAxis = spheroid(2); + elseif (spheroid(2) < 1) + E.Eccentricity = spheroid(2); + E.SemimajorAxis = spheroid(1); + else + error("sph_chk: eccentricity expected to be between 0 and 1"); + endif + elseif (ischar (spheroid) || isnumeric (spheroid) && isscalar (spheroid)) + E = referenceEllipsoid (spheroid); + endif + +endfunction + diff --git a/inst/str2angle.m b/inst/str2angle.m index 7b51d0b..79a2120 100644 --- a/inst/str2angle.m +++ b/inst/str2angle.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 Philip Nienhuis +## Copyright (C) 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 @@ -82,7 +82,7 @@ function deg = str2angle (txt, verbose = 0) - fmt = [ '([-+]?[0123456789]{1,3})([^+-]\s?)([+-]?[0123456789]{2})' ... + fmt = [ '([-+]?[0123456789]{1,3})([^+-]\s?|°\s?)([+-]?[0123456789]{2})' ... '[''mM]\s?([+-]?[0123456789\.].*?)((?:["sS]|'''')[eEnNsSwW]?)' ]; if (iscellstr (txt)) @@ -103,7 +103,7 @@ function deg = str2angle (txt, verbose = 0) endif deg(1, ierr) = deg(3, ierr) = deg(4, ierr) = 0; endif - deg(1, :) = (abs (deg(2, :)) .+ (deg(3, :) .+ deg(4, :) / 60) / 60) .* sign (deg(2, :)); + deg(1, :) = (abs (deg(2, :)) + (deg(3, :) + deg(4, :) / 60) / 60) .* sign (deg(2, :)); deg(2:end, :) = []; deg(1, ierr) = NaN; ## Set coordinates in western or southern hemisphere to negative. diff --git a/inst/toDegrees.m b/inst/toDegrees.m index ecc9f06..5577990 100644 --- a/inst/toDegrees.m +++ b/inst/toDegrees.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/toRadians.m b/inst/toRadians.m index 96590d1..b20813b 100644 --- a/inst/toRadians.m +++ b/inst/toRadians.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/unitsratio.m b/inst/unitsratio.m index ac26500..c19aea9 100644 --- a/inst/unitsratio.m +++ b/inst/unitsratio.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 diff --git a/inst/utmzone.m b/inst/utmzone.m index 61437c7..22d252a 100644 --- a/inst/utmzone.m +++ b/inst/utmzone.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 diff --git a/inst/validateLengthUnit.m b/inst/validateLengthUnit.m index 7d1c616..4bf892c 100644 --- a/inst/validateLengthUnit.m +++ b/inst/validateLengthUnit.m @@ -1,4 +1,4 @@ -## Copyright (C) 2014-2020 Carnë Draug <carandraug@octave.org> +## Copyright (C) 2014-2022 Carnë Draug <carandraug@octave.org> ## ## 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 @@ -176,7 +176,9 @@ function std_unit = validateLengthUnit (unit, varargin) ## Built start of error message from the extra optional arguments func_name = ""; arg_id = "input"; - if (nargin > 1) + if (nargin < 1) + print_usage () + elseif (nargin > 1) second = varargin{1}; if (ischar (second)) func_name = [second ": "]; diff --git a/inst/vincenty.m b/inst/vincenty.m new file mode 100644 index 0000000..d9423f2 --- /dev/null +++ b/inst/vincenty.m @@ -0,0 +1,169 @@ +## Copyright (C) 2014-2022 Alfredo Foltran <alfoltran@gmail.com> +## +## 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 meters) between two (sets of) locations on +## an ellipsoid. +## +## The formula devised by Thaddeus Vincenty is used with an accurate +## ellipsoidal model of the earth (@var{ellipsoid}). The default ellipsoidal +## model is 'WGS84', 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. +## Optional argument @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 <alfoltran@gmail.com> +## +## Octave style fixes and some re-coding to avoid sneaky errors by +## Philip Nienhuis <prnienhuis@users.sf.net> + +function [dist, az] = vincenty (pt1, pt2, ellipsoid) + + if (nargin < 3) + ellipsoid = referenceEllipsoid (7030); + endif + + major = ellipsoid.SemimajorAxis; + minor = ellipsoid.SemiminorAxis; + f = ellipsoid.Flattening; + ## Avoid confusion of length units impsed by ellipsoid, standardize on meters + if (isfield (ellipsoid, "LengthUnit") && ! isempty (ellipsoid.LengthUnit)) + major *= unitsratio ("meters", ellipsoid.LengthUnit); + minor *= unitsratio ("meters", ellipsoid.LengthUnit); + endif + + 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++; + [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 (abs (sin_sigma < eps)) + 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 (abs (cos_sq_alpha > eps)) + 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 alpha2]); + endif + +endfunction + + +%!test +%! assert (vincenty ([37, -76], [37, -9]), 5830081.063, 1e-2); + +%!test +%! assert (vincenty ([37, -76], [67, -76], referenceEllipsoid (7019)), ... +%! 3337842.871, 1e-2) diff --git a/inst/vincentyDirect.m b/inst/vincentyDirect.m new file mode 100644 index 0000000..9b910f5 --- /dev/null +++ b/inst/vincentyDirect.m @@ -0,0 +1,117 @@ +## Copyright (C) 2014-2022 Alfredo Foltran <alfoltran@gmail.com> +## Copyright (C) 2021-2022 Philip Nienhuis (<prnienhuis@users.sf.net> +## +## 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{lato}, @var{lono}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim}, @var{ellipsoid})} +## @deftypefnx {Function File} {[@var{lato}, @var{lono}, @var{azo}] = } {vincentyDirect(@var{lat}, @var{lon}, @var{range}, @var{azi}, @var{dim}, @var{ellipsoid})} +## +## Compute the coordinates of the end-point of a displacement on a geodesic. +## @var{lat}, @var{lon} are the coordinates of the starting point, @var{range} +## is the covered distance of the displacements along a specified geodesic and +## @var{azi} is the direction of the displacement relative to the North. +## The units of all input and output parameters must be 'radians' and/or the +## length unit of @var{range} should match that of the used ellipsoid. +## +## The possible values for @var{dim} are 'angle' (default) or 'length'. +## +## 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 +## +## @seealso{geodeticfwd, meridianfwd, reckon, referenceEllipsoid, vincenty} +## @end deftypefn + +## Author: Alfredo Foltran <alfoltran@gmail.com> +## Created: 2014-04-13 + +function [lato, lono, azo] = vincentyDirect (lat, lon, rng, azi, dim = "angle", ellipsoid) + + if (nargin < 6) + ellipsoid = referenceEllipsoid (7030); + endif + + major = ellipsoid.SemimajorAxis; + minor = ellipsoid.SemiminorAxis; + f = ellipsoid.Flattening; + + iter_limit = 20; + + tanU1 = (1 - f) * tan (lat); + U1 = atan (tanU1); + sigma1 = atan2 (tanU1, cos (azi)); + cosU1 = cos (U1); + sinAlpha = cosU1 * sin (azi); + cos2Alpha = (1 - sinAlpha) * (1 + sinAlpha); + u2 = cos2Alpha * (major ^ 2 - minor ^ 2) / minor ^ 2; + + A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))); + B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2))); + + if (strcmpi (dim, "length")) + sigma = rng / (minor * A); + lastSigma = sigma + 1; + i = 0; + while (abs (lastSigma - sigma) > 10e-12 && i <= iter_limit) + i++; + lastSigma = sigma; + doubleSigmaM = 2 * sigma1 + sigma; + deltaSigma = B * sin (sigma) * (cos (doubleSigmaM) + ... + 0.25 * B * (cos (sigma) * (-1 + 2 * cos (doubleSigmaM) ^ 2) ... + - 1/6 * B * cos (doubleSigmaM) * (-3 + 4 * sin (sigma) ^ 2) ... + * (-3 * 4 * cos (doubleSigmaM) ^ 2))); + sigma = rng / (minor * A) + deltaSigma; + endwhile + + if (i > iter_limit) + error ("Direct Vincenty's formulae failed to converge!"); + endif + elseif (strcmpi (dim, "angle")) + sigma = rng; + else + error ("Parameter \"dim\" must be \"angle\" or \"length\"!"); + endif + + doubleSigmaM = 2 * sigma1 + sigma; + sinU1 = sin (U1); + lato = atan2 (sinU1 * cos (sigma) + cosU1 * sin (sigma) * cos (azi), ... + (1 - f) * sqrt (sinAlpha ^ 2 + (sinU1 * sin (sigma) - ... + cosU1 * cos (sigma) * cos (azi)) ^ 2)); + lambda = atan2 (sin (sigma) * sin (azi), ... + cosU1 * cos (sigma) - sinU1 * sin (sigma) * cos(azi)); + C = f/16 * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha)); + L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sin (sigma) * ... + (cos (doubleSigmaM) + C * cos (sigma) * (-1 + 2 * cos (doubleSigmaM) ^ 2))); + lono = L + lon; + lono = wrapToPi (lono); + + if (nargout() > 2) + azo = atan2 (sinAlpha, -sinU1 * sin (sigma) + cosU1 * cos (sigma) * cos (azi)); + endif + +endfunction diff --git a/inst/wgs84Ellipsoid.m b/inst/wgs84Ellipsoid.m index 9235449..a7bd613 100644 --- a/inst/wgs84Ellipsoid.m +++ b/inst/wgs84Ellipsoid.m @@ -1,4 +1,4 @@ -## Copyright (C) 2018-2020 Philip Nienhuis +## Copyright (C) 2018-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 diff --git a/inst/wrapTo180.m b/inst/wrapTo180.m index 7b2798e..e7b5ea9 100644 --- a/inst/wrapTo180.m +++ b/inst/wrapTo180.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Oscar Monerris Belda +## Copyright (C) 2015-2022 Oscar Monerris Belda ## ## 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 diff --git a/inst/wrapTo2Pi.m b/inst/wrapTo2Pi.m index f6901be..8e01162 100644 --- a/inst/wrapTo2Pi.m +++ b/inst/wrapTo2Pi.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Oscar Monerris Belda +## Copyright (C) 2015-2022 Oscar Monerris Belda ## ## 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 diff --git a/inst/wrapTo360.m b/inst/wrapTo360.m index 41f2c72..c50c75a 100644 --- a/inst/wrapTo360.m +++ b/inst/wrapTo360.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Oscar Monerris Belda +## Copyright (C) 2015-2022 Oscar Monerris Belda ## ## 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 diff --git a/inst/wrapToPi.m b/inst/wrapToPi.m index c1886d8..6aedf92 100644 --- a/inst/wrapToPi.m +++ b/inst/wrapToPi.m @@ -1,4 +1,4 @@ -## Copyright (C) 2015-2020 Oscar Monerris Belda +## Copyright (C) 2015-2022 Oscar Monerris Belda ## ## 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 |