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