1 ## Copyright (C) 2008 Alexander Barth <barth.alexander@gmail.com>
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.
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.
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
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)
27 ## This function can also be used to define a spherical coordinate system
28 ## with rotated poles.
31 ## Author: Alexander Barth <barth.alexander@gmail.com>
33 function [lato,lono] = reckon(varargin);
37 [reg,prop] = parseparams(varargin);
55 reg{i} = repmat(reg{i},sz);
56 elseif !isequal(size(reg{i}),sz)
63 elseif length(prop) > 1
64 error("reckon: wrong number of type of arguments");
72 if strcmp(units,"degrees")
74 elseif strcmp(units,"radians")
77 error(["reckon: unknown units: " units]);
87 lato = pi/2 - acos(sin(lat).*cos(range) + cos(lat).*sin(range).*cos(azimuth));
89 cos_gamma = (cos(range) - sin(lato).*sin(lat))./(cos(lato).*cos(lat));
90 sin_gamma = sin(azimuth).*sin(range)./cos(lato);
92 gamma = atan2(sin_gamma,cos_gamma);
96 ## bring the lono in the interval [-pi pi[
98 lono = mod(lono+pi,2*pi)-pi;
100 ## convert to degrees
109 %! [lato,lono] = reckon(30,-80,20,40);
110 %! assert(lato,44.16661401448592,1e-10)
111 %! assert(lono,-62.15251496909770,1e-10)
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)
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)