octave_packages/mapping-1.0.7/reckon.m

   1 ## Copyright (C) 2008 Alexander Barth <barth.alexander@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, write to the Free Software
  15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} {[@var{lato},@var{lono}] = } reckon(@var{lat},@var{lon},@var{range},@var{azimuth})
  19 ## @deftypefnx {Function File} {[@var{lato},@var{lono}] = } reckon(@var{lat},@var{lon},@var{range},@var{azimuth},@var{units})
  20 ## Compute the coordinates of the end-point of a displacement on a
  21 ## sphere. @var{lat},@var{lon} are the coordinates of the starting point, @var{range}
  22 ## is the covered distance of the displacements along a great circle and
  23 ## @var{azimuth} is the direction of the displacement relative to the North.
  24 ## The units of all input and output parameters can be either 'degrees' (default)
  25 ## or 'radians'.
  26 ##
  27 ## This function can also be used to define a spherical coordinate system
  28 ## with rotated poles.
  29 ## @end deftypefn
  30
  31 ## Author: Alexander Barth <barth.alexander@gmail.com>
  32
  33 function [lato,lono] = reckon(varargin);
  34
  35   units = "degrees";
  36
  37   [reg,prop] = parseparams(varargin);
  38
  39   ## Input checking
  40   if length(reg) != 4
  41     print_usage ();
  42   endif
  43
  44   sz = [1 1];
  45
  46   for i=1:4
  47     if !isscalar(reg{i})
  48       sz = size(reg{i});
  49       break;
  50     endif
  51   endfor
  52
  53   for i=1:4
  54     if isscalar(reg{i})
  55       reg{i} = repmat(reg{i},sz);
  56     elseif !isequal(size(reg{i}),sz)
  57       print_usage();
  58     endif
  59   endfor
  60
  61   if length(prop) == 1
  62     units = prop{1};
  63   elseif length(prop) > 1
  64     error("reckon: wrong number of type of arguments");
  65   end
  66
  67   lat = reg{1};
  68   lon = reg{2};
  69   range = reg{3};
  70   azimuth = reg{4};
  71
  72   if strcmp(units,"degrees")
  73     d = pi/180;
  74   elseif strcmp(units,"radians")
  75     d = 1;
  76   else
  77     error(["reckon: unknown units: " units]);
  78   endif
  79
  80   ## convert to radians
  81
  82   lat = lat*d;
  83   lon = lon*d;
  84   range = range*d;
  85   azimuth = azimuth*d;
  86
  87   lato = pi/2 - acos(sin(lat).*cos(range) + cos(lat).*sin(range).*cos(azimuth));
  88
  89   cos_gamma = (cos(range) - sin(lato).*sin(lat))./(cos(lato).*cos(lat));
  90   sin_gamma = sin(azimuth).*sin(range)./cos(lato);
  91
  92   gamma = atan2(sin_gamma,cos_gamma);
  93
  94   lono = lon + gamma;
  95
  96   ## bring the lono in the interval [-pi pi[
  97
  98   lono = mod(lono+pi,2*pi)-pi;
  99
 100   ## convert to degrees
 101
 102   lono = lono/d;
 103   lato = lato/d;
 104
 105 endfunction
 106
 107
 108 %!test
 109 %! [lato,lono] = reckon(30,-80,20,40);
 110 %! assert(lato,44.16661401448592,1e-10)
 111 %! assert(lono,-62.15251496909770,1e-10)
 112
 113 %!test
 114 %! [lato,lono] = reckon(-30,80,[5 10],[40 45]);
 115 %! assert(lato,[-26.12155703039504 -22.70996703614572],1e-10)
 116 %! assert(lono,[83.57732793979254  87.64920016442251],1e-10)
 117
 118 %!test
 119 %! [lato,lono] = reckon([-30 31],[80 81],[5 10],[40 45]);
 120 %! assert(lato,[-26.12155703039504  37.76782079033356],1e-10)
 121 %! assert(lono,[83.57732793979254  89.93590456974810],1e-10)