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
%     other details make this routine different from an earlier version 
%     ellips used at yale.
%     6377. is the equatorial radius of the earth minus the depth of the
%     sound channel(=1200meters).
%     6378 is the equatorial radius of the earth, approximately.
%
%     changed to matlab by c.schmid march 94
%     Version 2.0          14 Nov 1996          R.Goldsmith
%             2.0  fixed angle computation as is was not done correctly
%                  if the angle approached 270 degrees.  270.1 flipped
%                  over to 90.1.   Might still be wrong for other cases
%                  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));
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);
                                         %   compute bearing...
   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
           angle = 2 * pi - angle;
           angle = angle / pi * 180;
       end
   end
else
    z = NaN;
    angle = NaN;
%        fprintf('error because distance <= 0')
end

abstand = z;

end
