function [phi,s,aik,bik,cik,dik]=join(ei,ni,ek,nk);

twopi  = 2*pi;
rhodeg = 180/pi;
rhosec = rhodeg*3600;

de  = ek-ei;           % east coordinate  difference
dn  = nk-ni;           % north coordinate difference
s2  = de*de + dn*dn;   % distance squared
phi = atan2(de,dn);    % bearing in radians
if phi<0
  phi = twopi +phi;
end  
phi = phi*rhodeg;           % bearing in seconds of arc
s   = sqrt(s2);             % distance in metres
aik = -de/(s2*100)*rhosec;  % seconds/cm
bik =  dn/(s2*100)*rhosec;  % seconds/cm
cik = dn/s;
dik = de/s;
