function [distkm,bear] = dist(lat1,lon1,lat2,lon2)
%
%function [distkm,bear] = dist(lat1,lon1,lat2,lon2)
%
% Function to calculate cumulative distance along a ship track
%    relative to the first point.
%
% Does not Assume Lat & Lon are in decimal degrees.
% If your positions are in decimal minutes it will call a subroutine
% called "ddmm2dec.m" to convert to decimal degrees.
% 
%  Limited  to less than 100 or so miles.  If over that limit
%  distorsions may exist due to sphere to plane projection.
%  For a more accurate calculation see "distcalc.m".
%
% Jim Manning 12/92 .. conversion of Fortran Code "Distan.for"
%
% next correct for longitude sign
if lon1>0;
  lon1=-1*lon1;
  lon2=-1*lon2;
end;

% first correct for position form
if lat1>90;
 [lat1,lon1]=ddmm2dec(lat1,lon1);
 [lat2,lon2]=ddmm2dec(lat2,lon2);
end;

pi=3.1415926;
pid180=pi/180.;
alat=lat1*pid180;
alon=lon1*pid180;
blat=lat2*pid180;
blon=lon2*pid180;
cix=cos(alat)*cos(alat)*cos(alon-blon)+sin(alat)*sin(alat);
% Note the problem with this calulation is we are assuming
% that the earth is a perfect sphere when the radius
% is 6378 at equator vs 6356 at poles
distkmx=6371*tan(sqrt(abs(1-cix*cix)/cix));
ciy=cos(alat)*cos(blat)*cos(alon-alon)+sin(alat)*sin(blat);
distkmy=6371*tan(sqrt(abs(1-ciy*ciy)/ciy));
if distkmx>0.001;
  bear=abs(atan(distkmy/distkmx));
  if blon>alon&blat>alat;
     bear=90-bear*180/pi;
  elseif  blon>alon&blat<=alat;
     bear=90+bear*180/pi;
  elseif  blon<alon&blat<=alat;
     bear=270-bear*180/pi;
  elseif  blon<alon&blat>alat;
     bear=270+bear*180/pi;
  end;
else;
  if blat>=alat
    bear=0;
  else;
    bear=180;
  end;
end;

ci=cos(alat)*cos(blat)*cos(alon-blon)+sin(alat)*sin(blat);
distkm=6371*tan(sqrt(abs(1-ci*ci)/ci));

%disp([num2str(distkm),' kilometers'])
%disp([num2str(distkm/1.8532),'nautical miles'])