octave_packages/mapping-1.0.7/azimuth.m

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