function GrosskreisbogenT

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

% Ausserdem benötigt: Terrainbase 5-Minute Global Digital Elevation Model (tbase.bin)
%                     des National Geophysical Data Center (NGDC)
%                     http://www.ngdc.noaa.gov/mgg/gravity/1999/data/global/tbase/

deg=pi/180;

% Datengitter  1/12°
% Quelle: http://www.ngdc.noaa.gov/mgg/gravity/1999/data/global/tbase/

id=fopen('tbase.bin','r');
A=fread(id,[4320 2160],'int16=>int16')';
fclose(id);
A=A(end:-1:1,:);
hS=2810;                   % Höhe Südpol
n=4;
d=n/12;                    % Vergröberung des Gitters auf n/12°
A=A(n:n:end,1:n:end);
H=double([repmat(hS,1,size(A,2)+1); A A(:,1)]);
clear A 
close('all') 
[L,B]=meshgrid([-180:d:180]*deg,[-90:d:90]*deg);
[X,Y,Z]=sph2cart(L,B,1); 
surf(X,Y,Z,H); 
hold('on') 
axis('equal') 
colormap(topocolors(min(min(H)),max(max(H)))) 
shading('interp')
caxis([min(min(H)) max(max(H))])
dms2rad=inline('(d+m/60+s/3600)*pi/180','d','m','s');
lambdaA=dms2rad( 13,15,43); phiA=dms2rad(53,33,22); % A=NB
lambdaB=dms2rad(139,32,22); phiB=dms2rad(35,40,21); % B=Tokio
[xA,yA,zA]=sph2cart(lambdaA,phiA,1);
[xB,yB,zB]=sph2cart(lambdaB,phiB,1);
e1S=[xA;yA;zA];
e3S=cross([xA;yA;zA],[xB;yB;zB]);
e3S=e3S/norm(e3S);
e2S=cross(e3S,e1S);
R=[e1S e2S e3S];
dS=acos(dot(e1S,[xB;yB;zB]));
tRef=0.2*dS;
t=linspace(0,dS,round(dS/deg));
Z=R*[cos(t);sin(t);zeros(size(t))];
ZRef=R*[cos(tRef);sin(tRef);0];
[lambdaRef phiRef]=cart2sph(ZRef(1),ZRef(2),ZRef(3));
T=-sin(tRef)*e1S+cos(tRef)*e2S;
Schrift(e1S,'{\itA}','erase','none','HorizontalAlignment','center','VerticalAlignment','top')
Schrift(ZRef,'{\itP}({\itt})','erase','none','HorizontalAlignment','left','VerticalAlignment','top')
Schrift(ZRef+T,'{\itP}''({\itt})','erase','none','HorizontalAlignment','right','VerticalAlignment','bottom')
Schrift([xB;yB;zB],'{\itB}','erase','none','HorizontalAlignment','left','VerticalAlignment','bottom')
plot3(ZRef(1),ZRef(2),ZRef(3),'.k','MarkerSize',30,'erase','none');
plot3(ZRef(1)+[0 T(1)],ZRef(2)+[0 T(2)],ZRef(3)+[0 T(3)],'k','linewidth',2.0,'erase','none');
Spitze(ZRef+T,T,0.1,30*deg,'k');
plot3(Z(1,:),Z(2,:),Z(3,:),'k','linewidth',2.0,'erase','none')
t=[0:360]*deg;
plot3(cos(phiRef)*cos(t),cos(phiRef)*sin(t),repmat(sin(phiRef),size(t)),'k','erase','none')
t=[-90:90]*deg;
plot3(cos(t)*cos(lambdaRef),cos(t)*sin(lambdaRef),sin(t),'k','erase','none')
L=0.8*(pi/2-phiRef);
Nord=[-sin(phiRef)*cos(lambdaRef);-sin(phiRef)*sin(lambdaRef);cos(phiRef)];
Nord=L*Nord;
plot3(ZRef(1)+[0 Nord(1)],ZRef(2)+[0 Nord(2)],ZRef(3)+[0 Nord(3)],'k','erase','none');
Spitze(ZRef+Nord,Nord,0.1,30*deg,'k');
Schrift(ZRef+Nord,'Nord  ','erase','none',...
  'HorizontalAlignment','right','VerticalAlignment','middle')
Ost=[-sin(lambdaRef);cos(lambdaRef);0];
Ost=L*Ost;
plot3(ZRef(1)+[0 Ost(1)],ZRef(2)+[0 Ost(2)],ZRef(3)+[0 Ost(3)],'k','erase','none');
Spitze(ZRef+Ost,Ost,0.1,30*deg,'k');
Schrift(ZRef+Ost+[0;0;0.05],' Ost','erase','none',...
  'HorizontalAlignment','left','VerticalAlignment','top')
rW=0.4;
alpha=atan2(dot(T,Ost),dot(T,Nord));
t=linspace(0,alpha,round(alpha/deg));
W=repmat(ZRef,size(t))+rW*(repmat(cos(t),3,1).*repmat(Nord,size(t))+repmat(sin(t),3,1).*repmat(Ost,size(t)));
plot3(W(1,:),W(2,:),W(3,:),'k','erase','none')
Spitze(W(:,end),cross(T,ZRef),0.05,45*deg,'k');
Schrift(W(:,round(end/2)),'\alpha^{ }','erase','none',...
  'HorizontalAlignment','right','VerticalAlignment','middle')
view(157,24)
axis([-0.7 0.7 -0.7 0.7 0.85 0.851])
axis('off') 
print('-deps','GrosskreisbogenT')

% Hilfsfunktionen zur Farbgestaltung

function C=topocolors(a,b)
% schwarz-weiss
q=0.3;
C=[q*gray(-a)+(1-q)*ones(-a,3)
  q*gray(b)+(1-q)*ones(b,3)];
% % Farbe
% C=[topomaptrans(1:36,1:-a)
%   topomaptrans(37:64,1:b)];
% % aufhellen
% a=0.3;
% C=a*C+(1-a)*ones(size(C));

function T=topomaptrans(I,J)
if isempty(J)
  T=[];
else
  load('topo','topomap1')
  K=polyval( polyfit(J([1 end]),I([1 end]),1) , J);
  T=[topomap1(I(1),:)
    interp1(I,topomap1(I,:),K(2:end-1))
    topomap1(I(end),:)];
end

function Schrift(P,txt,varargin)
text(P(1),P(2),P(3),txt,'FontName','TimesNewRoman','FontSize',18,varargin{:});

function Spitze(P,v,L,lambda,varargin);
deg=pi/180;
if nargin<=2
  L=1;
end
if nargin<=3
  lambda=7.5*deg;
end
v=L*v./repmat(sqrt(dot(v,v)),size(P,1),1);
Q=P-v;
t=[0:6:360]*deg;
for k=1:size(P,2)
  N=tan(lambda/2)*null(v(:,k)');
  for w=2:length(t)
    plot3(...
      [P(1,k) Q(1,k)+cos(t(w))*L(:,k)*N(1,1)+sin(t(w))*L(:,k)*N(1,2) Q(1,k)+cos(t(w-1))*L(:,k)*N(1,1)+sin(t(w-1))*L(:,k)*N(1,2) P(1,k)],...
      [P(2,k) Q(2,k)+cos(t(w))*L(:,k)*N(2,1)+sin(t(w))*L(:,k)*N(2,2) Q(2,k)+cos(t(w-1))*L(:,k)*N(2,1)+sin(t(w-1))*L(:,k)*N(2,2) P(2,k)],...
      [P(3,k) Q(3,k)+cos(t(w))*L(:,k)*N(3,1)+sin(t(w))*L(:,k)*N(3,2) Q(3,k)+cos(t(w-1))*L(:,k)*N(3,1)+sin(t(w-1))*L(:,k)*N(3,2) P(3,k)],...
      varargin{:});
    end
  end
