diff --git a/lib/+artoa/+data/calculateEllipk2.m b/lib/+artoa/+data/calculateEllipk2.m index 27ec762fc69c56153b35977b7aeaeff8fda62033..981d0cdffb6707b8a7d0a89b917ee5639cbe79bb 100644 --- a/lib/+artoa/+data/calculateEllipk2.m +++ b/lib/+artoa/+data/calculateEllipk2.m @@ -1,8 +1,6 @@ -function [abstand,ang] = calculateEllipk2(st,sc) -%function [abstand,ang] = ellipk2(st,sc) -% -% routine computes great arc distance z in kilometers on a spheroid -% with flattening =1/295=.00339 using the andoyer-lambert approximation. +function [abstand, angle] = calculateEllipk2(st, sc) +%CALCULATEELLIPK2 Computes great arc distance z in kilometers on a spheroid. +% Does flattening =1/295=.00339 using the andoyer-lambert approximation. % z is the distance from sc to st and ang is the bearing from sc to % st in true radians rel. to north. % the coordinates for sc and st(lat,long) must be in radians. this and @@ -20,41 +18,42 @@ function [abstand,ang] = calculateEllipk2(st,sc) % but works for my case. % -snst=sin(st(:,1)); -snsc=sin(sc(:,1)); -cfz = snst .* snsc + cos(st(:,1)) .* cos(sc(:,1)) .* cos(st(:,2)-sc(:,2)); +snst = sin(st(:, 1)); +snsc = sin(sc(:, 1)); +cfz = snst .* snsc + cos(st(:, 1)) .* cos(sc(:, 1)) .* cos(st(:, 2) - sc(:, 2)); z = 1 - cfz .* cfz; if z > 0 - sfz=sqrt(z); - z=atan2(sfz,cfz); - a=snst-snsc; - b=snst+snsc; - delz = (z+3*sfz) ./ (1.-cfz) .*a .*a; - delz = .84752e-3 * (delz + (z-3.*sfz) ./ (1.+cfz) .*b .*b); - z=6378.0*(z-delz); + sfz = sqrt(z); + z = atan2(sfz, cfz); + a = snst - snsc; + b = snst + snsc; + delz = (z + 3 * sfz) ./ (1. - cfz) .* a .* a; + delz = .84752e-3 * (delz + (z-3. * sfz) ./ (1. + cfz) .* b .* b); + z = 6378.0 * (z - delz); % compute bearing... - ang1 = cos(sc(:,1)) .* tan(st(:,1)) - sin(sc(:,1)) .* cos(sc(:,2)-st(:,2)); - if ang1 ~= 0, - ang = atan(sin(sc(:,2)-st(:,2))./ang1); - for ii=1:length(st(:,1)), - if ang1 < 0, - if sc(ii,1) > st(ii,1), - ang = ang + pi; - end - if sc(ii,1) < st(ii,1), - ang = ang - pi; + ang1 = cos(sc(:, 1)) .* tan(st(:, 1)) - sin(sc(:, 1)) .* cos(sc(:, 2) - st(:, 2)); + if ang1 ~= 0 + angle = atan(sin(sc(:,2)-st(:,2))./ang1); + for ii=1:length(st(:,1)) + if ang1 < 0 + if sc(ii,1) > st(ii,1) + angle = angle + pi; + end + if sc(ii,1) < st(ii,1) + angle = angle - pi; + end end - end - ang = 2 * pi - ang; - ang = ang / pi * 180; + angle = 2 * pi - angle; + angle = angle / pi * 180; end - end - else - - z=NaN; - ang = NaN; + end +else + z = NaN; + angle = NaN; % fprintf('error because distance <= 0') - end -abstand=z; -%the end +end + +abstand = z; + +end diff --git a/lib/+artoa/+data/calculateGeodist.m b/lib/+artoa/+data/calculateGeodist.m index b73edd57d8ab8e7fe0170ea4e9342089fc7925e5..87be99ad42b00a0301bc4f5ac1184677c781e912 100644 --- a/lib/+artoa/+data/calculateGeodist.m +++ b/lib/+artoa/+data/calculateGeodist.m @@ -1,49 +1,52 @@ -% GEODIST calculate the distance and bearing between two positions +function [distance, angle] = calculateGeodist(pPosition1, pPosition2, pUnit) +%CALCULATEGEODIST Calculate the distance and bearing between the two positions. +% Calculates the distance and bearing between pos1 and pos2. The optional +% parameter 'pUnit' specifies if the positions are in degrees ('degree') or +% in radiant ('rad'). If unit is not given, 'degree' would be the default. % -% SYNOPSIS: -% function [dist,ang] = geodist(pos1,pos2,unit) +% Goedist returns distance from pPosition1 to pPosition2 (dist) and the bearing +% from pPosition1 to pPosition2 (ang) in degrees rel. to north. % -% DESCRIPTION: -% Calculates the distance and bearing between pos1 and pos2. The optional -% parameter 'unit' specifies if the positions are in degrees ('degree') or -% in radiant ('rad'). If unit is not given, 'degree' would be the default. +% Parameters: +% pPosition1 (double) 1D-Vector with two values, latitude and +% longitude. +% pPosition2 (double) 1D-Vector with two values, latitude and +% longitude. +% pUnit (char) Specifies if positions are in degrees or +% radiant. % -% Goedist returns distance from pos1 to pos2 (dist) and the bearing from -% pos1 to pos2 (ang) in degrees rel. to north. -% -% -% rlat=[pos1(1),pos2(1)] -% rlon=[pos1(2),pos2(2)] -% -% calls function ellipk2 to calculate distance and bearing - -function [dist,ang] = calculateGeodist(pos1,pos2,unit) +% Returns: +% distance (double) The calculated distance. +% angle (double) The calculated bearing in degrees rel. to +% north. +%% Parameter check if nargin < 3 - unit='degree'; + pUnit='degree'; end -rlat=[pos1(1),pos2(1)]; -rlon=[pos1(2),pos2(2)]; +%% Convert input +rlat = [pPosition1(1), pPosition2(1)]; +rlon = [pPosition1(2), pPosition2(2)]; -if strcmp(unit,'degree') +if strcmp(pUnit,'degree') % rad = degree/180*pi ! rlat = rlat / 180 * pi; rlon = rlon / 180 * pi; end -s=size(rlat); -if s(1) > 1 +rlat_size = size(rlat); +if rlat_size(1) > 1 % make sure that rlat and rlon have more rows than columns - if s(1) < s(2) - rlat=rlat'; - rlon=rlon'; + if rlat_size(1) < rlat_size(2) + rlat = rlat'; + rlon = rlon'; end end -st=[rlat(:,1) rlon(:,1)]; -sc=[rlat(:,2) rlon(:,2)]; -[dist,ang] = artoa.data.calculateEllipk2(st,sc); -% dist is the distance from sc to st and ang is the bearing from sc to -% st in degrees rel. to north. +st = [rlat(:,1) rlon(:,1)]; +sc = [rlat(:,2) rlon(:,2)]; + +[distance, angle] = artoa.data.calculateEllipk2(st, sc); +end