## 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 point between two great circles ## ## Input: ## @itemize ## @item ## @var{lat1}, @var{lon1}, @var{az1}: latitude, longitude, and azimuth of great circle #1 ## @end item ## ## @item ## @var{lat2}, @var{lon2}, @var{az2}: latitude, longitude, and azimuth of great circle #2 ## @end item ## ## @item ## @var{angleUnit}: string for angular units ('degrees' or 'radians', ## case-insensitive, just the first character will do). Default is 'degrees'. ## @end item ## @end itemize ## ## 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 (lat1, lon1, az1, lat2, lon2, az2, angleUnit = "degrees") if (nargin < 6) print_usage(); endif if (! isnumeric (lat1) || ! isreal (lat1) || ... ! isnumeric (lon1) || ! isreal (lon1) || ... ! isnumeric (az1) || ! isreal (az1) || ... ! isnumeric (lat2) || ! isreal (lat2) || ... ! isnumeric (lon2) || ! isreal (lon2) || ... ! isnumeric (az2) || ! isreal (az2) ) error ("gcxgc : numeric values expected for first five inputs"); endif if (! ischar (angleUnit)) error ("gcxgc: character value expected for 'angleUnit'"); elseif (strncmpi (angleUnit, "degrees", length (angleUnit))) vect=deg2rad([lat1 lon1 az1 lat2 lon2 az2]); elseif (! strncmpi (angleUnit, "radians", length (angleUnit))) error ("gcxgc: illegal input for 'angleUnit'"); else vect=[lat1 lon1 az1 lat2 lon2 az2]; endif ## Algorithm from http://www.movable-type.co.uk/scripts/latlong.html del_lat=abs(vect(1)-vect(4)); del_lon=abs(vect(2)-vect(5)); del_12=2*asin(sqrt(sin(del_lat/2)^2+cos(vect(1))*cos(vect(4))*sin(del_lon/2)^2)); theta_a=acos((sin(vect(4))-sin(vect(1))*cos(del_12))/(sin(del_12)*cos(vect(1)))); theta_b=acos((sin(vect(1))-sin(vect(4))*cos(del_12))/(sin(del_12)*cos(vect(4)))); if sin (vect(5) - vect(2)) > 0 theta_12 = theta_a; theta_21 = 2 * pi - theta_b; else theta_12 = 2 * pi - theta_a; theta_21 = theta_b; endif alpha1 = vect(3) - theta_12; alpha2 = theta_21 - vect(6); alpha3 = acos(-cos (alpha1) * cos(alpha2) + ... sin (alpha1) * sin (alpha2) * cos (del_12)); del_13 = atan2 (sin (del_12) * sin (alpha1) * sin (alpha2),... cos(alpha2)+cos(alpha1)*cos(alpha3)); lat3 = asin (sin (vect(1)) * cos (del_13) + ... cos (vect(1)) * sin (del_13) * cos (vect(3))); del_lon13 = atan2 (sin (vect(3)) * sin (del_13) * cos (vect(1)), ... cos (del_13) - sin (vect(1)) * sin (lat3)); lon3 = vect(2) + del_lon13; [alat3 alon3] = antipode (lat3, lon3, 'r'); lat = [lat3 alat3]; lon = [lon3 alon3]; if (sin (alpha1) * sin(alpha2)) < 0 ## Causes ambiguity lon = [0 pi]; endif if (strncmpi (angleUnit, "degrees", length (angleUnit))) lat = rad2deg (lat); lon = rad2deg (lon); endif 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) %!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")