Restructured calculateEllipk2 and calculateGeodist.

lib/+artoa/+data/calculateEllipk2.m

-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

lib/+artoa/+data/calculateGeodist.m

-% 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