function NBTokio

% http://geometrie.hs-nb.de/MATLAB-Programme/NBTokio.m   Bearbeitungsstand: 25.3.2005
% (C) Martin Nitschke, nitschke@hs-nb.de

clc
deg=pi/180;
R=6371;        % Erdradius
dms2rad=inline('(d+m/60+s/3600)*pi/180','d','m','s');  
               % (deg,min,sec) -> rad
lambdaA= dms2rad(13,15,43); phiA=dms2rad(53,33,22);  
               % Neubrandenburg
lambdaB=dms2rad(139,32,22); phiB=dms2rad(35,40,21);  
               % Tokio
delta=500/R;   % Abstand der Zwischenpunkte
[xA,yA,zA]=sph2cart(lambdaA,phiA,1);
[xB,yB,zB]=sph2cart(lambdaB,phiB,1);
e1Strich=[xA;yA;zA];
ab=cross(e1Strich,[xB;yB;zB]);
q=norm(ab);
e3Strich=ab/q;
e2Strich=cross(e3Strich,e1Strich);
dS=acos(dot(e1Strich,[xB;yB;zB]));
tU=atan(e2Strich(3)/e1Strich(3));
tU=tU+2*pi*(tU<0);
for t=unique([tU 0:delta:dS dS])
  s=R*t;
  Z=  cos(t)*e1Strich+sin(t)*e2Strich;
  dZ=-sin(t)*e1Strich+cos(t)*e2Strich;
  [lambda,phi]=cart2sph(Z(1),Z(2),Z(3));
  [dL,mL,sL]=rad2dms(abs(lambda));
  [dB,mB,sB]=rad2dms(abs(phi));
  eN=[-sin(phi)*cos(lambda)
      -sin(phi)*sin(lambda)
       cos(phi)];
  eO=[-sin(lambda)
       cos(lambda)
       0];
  alpha=atan2(dot(eO,dZ),dot(eN,dZ));
  alpha=alpha+2*pi*(alpha<0);
  [dA,mA,sA]=rad2dms(alpha);
  fprintf(['s=%5.0fkm lambda=%4.0fº%02.0f''%02.0f" '...
           'phi=%3.0fº%02.0f''%02.0f" ' ...
           'alpha=%3.0fº%02.0f''%02.0f"'],...
   s,dL*sign(dL),mL,sL,dB*sign(dB),mB,sB,dA,mA,sA)
  if abs(t-tU) < sqrt(eps)
    fprintf(' U');
  end
  fprintf('\n')
end

function [d,m,s]=rad2dms(rad)
  D=rad/pi*180;
  d=fix(D);
  M=(D-d)*60;
  m=fix(M);
  s=(M-m)*60;
