## 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} =} gcxgc (@var{lat1}, @var{lon1}, @var{az1}, @var{lat2}, @var{lon2}, @var{az2})
## @deftypefnx {Functional File} {@var{lat}, @var{lon} =} 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.
##
## 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] = gcxgc (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
## 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
endif
if (! ischar (angleUnit))
error ("gcxgc: character value expected for 'angleUnit'");
elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7)))
vect = deg2rad ([varargin{:}]);
if (vect(:, 1) > pi / 2 || vect(:, 1) < - pi / 2 ||...
vect(:, 4) > pi / 2 || vect(:, 4) < - pi / 2 ||...
vect(:, 2) > pi || vect(:, 2) < -pi || ...
vect(:, 5) > pi || vect(:, 5) < -pi)
error("gcxgc: values out of expected range lat [-90, 90] and long [-180, 180]")
endif
elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7)))
vect = [varargin{:}];
if (vect(:, 1) > pi / 2 || vect(:, 1) < - pi / 2 ||...
vect(:, 4) > pi / 2 || vect(:, 4) < - pi / 2 ||...
vect(:, 2) > pi || vect(:, 2) < -pi || ...
vect(:, 5) > pi || vect(:, 5) < -pi)
error("gcxgc: values out of expected range lat [-pi/2, pi/2] and long [-pi, pi]")
endif
else
error ("gcxgc: illegal input for 'angleUnit'");
endif
#Algorithm from Michael Geyer
#EARTH-REFERENCED AIRCRAFT NAVIGATION AND SURVEILLANCE ANALYSIS" pg 111 & 123
c1 = lla2xyz (vect(:, 1), vect(:, 2), vect(:, 3));
c2 = lla2xyz (vect(:, 4), vect(:, 5), vect(:, 6));
N = cross (c1, c2);
lat3 = atan2 (N(3), hypot (N(1), N(2)));
del_su = vect(:,5) - vect(:,2);
phi_su = atan2 ((cos (vect(:,4)) .* sin (del_su)), ...
(sin (vect(:, 4)) .* cos (vect(:, 1)) .- ...
cos (vect(:, 4)) .* sin (vect(:, 1)) .* cos(del_su)));
phi_us = atan2 ((cos (vect(:,1)) .* sin (-del_su)), ...
(sin (vect(:, 1)) .* cos (vect(:, 4)) .- ...
cos (vect(:, 1)) .* sin (vect(:, 4)) .* cos(-del_su)));
beta_u = min ([abs(phi_su .- vect(:,3)),
abs(phi_su .- vect(:,3) .+ 2 * pi),
abs(phi_su .- vect(:,3) .- 2*pi)]);
beta_s = min ([abs(phi_us .- vect(:,6)),
abs(phi_us .- vect(:,6) .+ 2 * pi),
abs(phi_us .- vect(:,6) .- 2*pi)]);
if (sind (rad2deg (vect(:,3))) == 0 && sind (rad2deg (vect(:,6))) == 0)
#Note: use sind because sin (pi) != 0
lon3 =0;
elseif ( beta_u == beta_s || beta_u + beta_s >= pi || beta_u + beta_s <= 0)
warning("gcxgc: Non-unique intersection.\n")
lat3 = NaN;
lon3 = NaN;
elseif ((phi_su + beta_u == vect(:, 3)) ||
(phi_su - beta_u == vect(:, 3)) ||
(phi_us - beta_s == vect(:, 6)) ||
(phi_us + beta_s == vect(:, 6)))
lon3 = atan2(N(2), N(1));
endif
[alat3 alon3] = antipode (lat3, lon3, "r");
lat = [lat3 alat3];
lon = [lon3 alon3];
if (strncmpi (angleUnit, "degrees", length (angleUnit)))
lat = rad2deg (lat);
lon = rad2deg (lon);
endif
endfunction
function [C] = lla2xyz (lat, long, az)
C(:,1) = cos (lat) ^ 2 .* sin (long) .* cos (az) .- sin (lat) .* ...
(sin (az) .* cos (long) .- sin (lat) .* cos (az) .* sin (long));
C(:,2) = -cos (lat) ^ 2 .* cos (long) .* cos (az) .- sin (lat) .* ...
(sin (az) .* sin (long) .+ sin (lat) .* cos (az) .* cos (long));
C(:,3) = cos (lat) .* sin (az);
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] = gcxgc (20, -5, 45, 30, 5, 15);
%! assert (lat3(1), 28.0620, 10-3)
%! assert (lon3(1), 4.4121, 10-3)
%!test
%! [lat3, lon3] = gcxgc (0, 0, 90, 0, 180, 80);
%! assert (lat3(1), 0, 10-3)
%! assert (lon3(1), 180, 10-3)
%!warning [lat3, lon3] = gcxgc (0, 0, 45, 0, 180, -45);
%!warning [lat3, lon3] = gcxgc (45, 45, 90, -45, -135, 90);
%!error gcxgc ("s", 0, 100, 10, 30, 0);
%!error gcxgc (3i, 0, 100, 10, 30, 0);
%!error gcxgc (50, "s", 100, 10, 30, 0);
%!error gcxgc (50, 2i, 10, 10, 30, 0);
%!error gcxgc (50, 0, "s", 10, 30, 0);
%!error gcxgc (50, 0, 100i, 10, 30, 0);
%!error gcxgc (50, 0, 100, "s", 30, 0);
%!error gcxgc (50, 0, 100, 10i, 30, 0);
%!error gcxgc (50, 0, 100, 10, "s", 0);
%!error gcxgc (50, 0, 100, 10, 30i, 0);
%!error gcxgc (50, 0, 100, 10, 30, "s");
%!error gcxgc (50, 0, 100, 10, 30, 2i);
%!error gcxgc (50, 0, 100, 10, 30, 0, "f");
%!error gcxgc (50, 0, 100, 10, 30, 0, "degreef");
%!error gcxgc (190, 0, 90, -90.000001, 180, 80);