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.
%     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,
     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;
           end
         end
       ang = 2 * pi - ang;
       ang = ang / pi * 180;
       end
     end
    else
          
     z=NaN;
     ang = NaN;
%        fprintf('error because distance <= 0')
     end
abstand=z;
%the end