## Copyright (C) 2020 ## ## This program is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program. If not, see . ## -*- texinfo -*- ## @deftypefn {Functional File} {@var{lat}, @var{lon}, @var{idl} =} gcxgc (@var{lat1}, @var{lon1}, @var{az1}, @var{lat2}, @var{lon2}, @var{az2}) ## @deftypefnx {Functional File} {@var{lat}, @var{lon}, @var{idl} =} gcxgc (@var{lat1}, @var{lon1}, @var{az1}, @var{lat2}, @var{lon2}, @var{az2}, @var{angleUnit}) ## Determines the intersection points between two great circles. ## ## Input: ## @itemize ## @item ## @var{lat1}, @var{lon1}, @var{az1}: latitude, longitude, and azimuth of ## great circle #1. These must be scalar values or vectors of equal length. ## @end item ## ## @item ## @var{lat2}, @var{lon2}, @var{az2}: latitude, longitude, and azimuth of ## great circle #2. These must be scalar values or vectors of equal length. ## @end item ## ## @item ## @var{angleUnit}: string for angular units ('degrees' or 'radians', ## case-insensitive, just the first character will do). Default is 'degrees'. ## @var{angleUnit} applies to all inputs and outputs. ## @end item ## @end itemize ## ## Output: ## If scalar values have been input, @var{lat} and @var{lon} are both 1x2 ## vectors. If vectors have been input @var{lat} and @var{lon} are Nx2 arrays ## where N is the number of great circle pairs. The results for multiple ## great circle pairs are concatenated vertically no matter the orientation of ## input vectors. ## ## Optional ouput @var{idl} lists pairs of coinciding great circles, if any. ## ## Example: ## @example ## lat1 = 51.8853; ## lon1 = 0.2545; ## az1 = 108.55; ## lat2 = 49.0034; ## lon2 = 2.5735; ## az2 = 32.44; ## [newlat, newlon] = gcxgc (lat1, lon1, az1, lat2, lon2, az2) ## newlat = ## 50.908 -50.908 ## newlon = ## 4.5086 -175.4914 ## @end example ## @end deftypefn function [lat, lon, idl] = gcxgc3 (varargin) if (nargin < 6) print_usage(); elseif (nargin == 6) angleUnit = "degrees"; else angleUnit = varargin{7}; endif if (! (all (cellfun ("isnumeric", varargin(1:6)) && ... all (cellfun ("isreal", varargin(1:6)))))) error ("gcxgc: numeric values expected for first six inputs"); endif isv = ! cellfun ("isscalar", varargin(1:6)); if (any (isv)) ## At least one of the location inputs is a vector. Check sizes numval = cellfun ("numel", varargin(isv)); if (any (diff (numval))) error ("gcxgc: all vector inputs must have same lengths"); endif nv = numval(1); ## Make sure all inputs are column vectors of same length for ii=1:6 if (isv(ii)) varargin(ii) = {varargin{ii}(:)}; else varargin(ii) = {(repmat (varargin{ii}, numval(1), 1))}; endif endfor else nv = 1; endif if (! ischar (angleUnit)) error ("gcxgc: character value expected for 'angleUnit'"); elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7))) vect = deg2rad ([varargin{:}]); elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7))) vect = [varargin{:}]; else error ("gcxgc: illegal input for 'angleUnit'"); endif [lat, lon] = get_intscs (vect); ## Check for coinciding great circles. Done by comparing (Lat,Lon) where ## azimuths are exactly EW / 90 degrees ## 1. Intersection with equator ## FIXME (speedup): check for lat = azimuth == 0 (as those ARE on equator) [~, loni1] = ... get_intscs ([(zeros (nv, 2)), (pi / 2 * ones (nv, 1)), vect(:, 1:3)]); [~, loni2] = ... get_intscs ([(zeros (nv, 2)), (pi / 2 * ones (nv, 1)), vect(:, 4:6)]); ## 2. Intersection of equator + great circles tru poles midway through loni* [lati1, loni1] = ... get_intscs ([(zeros (nv, 1)), (mean (loni1')'), (zeros (nv, 1)), vect(:, 1:3)]); [lati2, loni2] = ... get_intscs ([(zeros (nv, 1)), (mean (loni2')'), (zeros (nv, 1)), vect(:, 4:6)]); ## We don't need to compare polar great circle axis; just comparing lati* ## and antipodes on N hemisphere will do. First select those coordinates id1 = lati1(:, 1) < 0; lati1(id1, 1) = lati1(id1, 2); loni1(id1, 1) = loni1(id1, 2); lati1 = lati1(:, 1); loni1 = loni1(:, 1); id2 = lati2(:, 1) < 0; lati2(id2, 1) = lati2(id2, 2); loni2(id2, 1) = loni2(id2, 2); lati2 = lati2(:, 1); loni2 = loni2(:, 1); ## Find out which lati & loni coincide. ("and()" accepts multidimensions) idl = and (abs (lati1 .- lati2) < 2 * eps, abs (loni1 .- loni2) < 2 * eps); ## Set output relating to coinciding great circles op NaN, NaN lat(idl, :) = NaN; lon(idl, :) = NaN; idl = find (idl); if (! isempty (idl)) warning ("gcxgc: non-unique intersection(s).\n") endif if (strncmpi (angleUnit, "degrees", length (angleUnit))) lat = rad2deg (lat); lon = rad2deg (lon); endif 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)) .* ... cos (vect(:, 2)) .* sin (vect(:, 3)); 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)) .* ... cos (vect(:, 5)) .* sin (vect(:, 6)); 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)); N = cross (c1, c2, 2); lat3 = atan2 (N(:, 3), hypot (N(:, 1), N(:, 2))); if (sind (rad2deg (vect(:, 3))) == 0 && sind (rad2deg (vect(:, 6))) == 0) #Note: use sind because sin (pi) != 0 lon3 = zeros (size (vect, 1)); else lon3 = atan2 (N(:, 2), N(:, 1)); endif [alat3 alon3] = antipode (lat3, lon3, "r"); lat = [lat3 alat3]; lon = [lon3 alon3]; endfunction %!test %! [lat3, lon3] = gcxgc ( 51.8853, 0.2545, 108.55, 49.0034, 2.5735, 32.44); %! assert (degrees2dms (lat3(1)), [50 54 27], 10-3) %! assert (degrees2dms (lon3(1)), [04 30 31], 10-3) %!test %! [lat3, lon3] = gcxgc3 (20, -5, 45, 30, 5, 15); %! assert (lat3(1), 28.0620, 10-3) %! assert (lon3(1), 4.4121, 10-3) %!warning gcxgc3 (0, 0, 45, 0, 180, -45); %!error gcxgc3 ("s", 0, 100, 10, 30, 0) %!error gcxgc3 (3i, 0, 100, 10, 30, 0) %!error gcxgc3 (50, "s", 100, 10, 30, 0) %!error gcxgc3 (50, 2i, 10, 10, 30, 0) %!error gcxgc3 (50, 0, "s", 10, 30, 0) %!error gcxgc3 (50, 0, 100i, 10, 30, 0) %!error gcxgc3 (50, 0, 100, "s", 30, 0) %!error gcxgc3 (50, 0, 100, 10i, 30, 0) %!error gcxgc3 (50, 0, 100, 10, "s", 0) %!error gcxgc3 (50, 0, 100, 10, 30i, 0) %!error gcxgc3 (50, 0, 100, 10, 30, "s") %!error gcxgc3 (50, 0, 100, 10, 30, 2i) %!error gcxgc3 (50, 0, 100, 10, 30, 0, "f") %!error gcxgc3 (50, 0, 100, 10, 30, 0, "degreef") %!error gcxgc3 ([50 0], 0, 0, 0, 0, [1 2 3])