1## Copyright (C) 2004-2020 Andrew Collier <abcollier@users.sourceforge.net> 2## Copyright (C) 2006-2020 Alexander Barth <abarth93@users.sourceforge.net> 3## 4## This program is free software; you can redistribute it and/or modify it under 5## the terms of the GNU General Public License as published by the Free Software 6## Foundation; either version 3 of the License, or (at your option) any later 7## version. 8## 9## This program is distributed in the hope that it will be useful, but WITHOUT 10## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 11## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 12## details. 13## 14## You should have received a copy of the GNU General Public License along with 15## this program; if not, see <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {} @var{az} = azimuth(@var{lat1},@var{lon1},@var{lat2},@var{lon2}) 19## @deftypefnx {Function File} {} @var{az} = azimuth(@var{lat1},@var{lon1},@var{lat2},@var{lon2},@var{units}) 20## @deftypefnx {Function File} {} @var{az} = azimuth(@var{pt1}, @var{pt2}) 21## @deftypefnx {Function File} {} @var{az} = azimuth(@var{pt1}, @var{pt2},@var{units}) 22## 23## Calculates the great circle azimuth from a point 1 to a point 2. 24## The latitude and longitude of these two points can either be given 25## independently or as columns of the matrices @var{pt1} and @var{pt2} in 26## the form [latitude longitude]. 27## 28## The units for the input coordinates and output angles can be 29## "degrees" (the default) or "radians". 30## 31## @example 32## >> azimuth([10,10], [10,40]) 33## ans = 87.336 34## >> azimuth([0,10], [0,40]) 35## ans = 90 36## >> azimuth(pi/4,0,pi/4,-pi/2,"radians") 37## ans = 5.3279 38## @end example 39## 40## @seealso{elevation,distance} 41## @end deftypefn 42 43## Author: Andrew Collier <abcollier@users.sourceforge.net> 44## Adapted-by: Alexander Barth <abarth93@users.sourceforge.net> 45 46## Uses "four-parts" formula. 47 48function az = azimuth(varargin) 49 ## default units are degrees 50 51 units = "degrees"; 52 53 [reg,prop] = parseparams(varargin); 54 55 if length(reg) == 2 56 pt1 = reg{1}; 57 pt2 = reg{2}; 58 59 a = pt1(:,1); 60 b = pt2(:,1); 61 C = pt2(:,2) - pt1(:,2); 62 elseif length(reg) == 4 63 a = reg{1}; 64 b = reg{3}; 65 C = reg{4} - reg{2}; 66 else 67 error("Wrong number of type of arguments"); 68 end 69 70 if length(prop) == 1 71 units = prop{1}; 72 73 if (~strcmp(units,"degrees") && ~strcmp(units,"radians")) 74 error("Only degrees and radians are allowed as units"); 75 end 76 elseif length(prop) > 1 77 error("Wrong number of type of arguments"); 78 end 79 80 if (strcmp(units,"degrees")) 81 a = deg2rad(a); 82 b = deg2rad(b); 83 C = deg2rad(C); 84 end 85 86 az = atan2(sin(C) , cos(a) .* tan(b) - sin(a) .* cos(C) ); 87 88 ## bring the angle in the interval [0 2*pi[ 89 90 az = mod(az,2*pi); 91 92 ## convert to degrees if desired 93 94 if (strcmp(units,"degrees")) 95 az = rad2deg(az); 96 end 97endfunction 98 99%!test 100%! assert(azimuth([10,10], [10,40]), 87.336, 1e-3) 101%! assert(azimuth([0,10], [0,40]), 90, 1e-3) 102%! assert(azimuth(pi/4,0,pi/4,-pi/2,"radians"), 5.3279, 1e-4) 103