function dSNBRio

% http://geometrie.hs-nb.de/MATLAB-Programme/dSNBRio.m   Bearbeitungsstand: 9.8.2004
% (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);
lambdaB=-dms2rad(43,13,22);
phiB   =-dms2rad(22,53,42);
[sZ,lambdaZ,phiZ,KursZ,sU,lambdaU,phiU,KursU]=...
  RouteBerechnen(lambdaA,phiA,lambdaB,phiB,250/6371);
[P1,P2,P3]=RouteDarstellen3D(lambdaA,phiA,lambdaB,phiB,lambdaZ,phiZ,lambdaU,phiU);
hold('off')
view(100,0);
axis([-0.9 0.9 -0.9 0.9 -1 1])
axis('off') 
print('-deps','dSNBRio')

% Berechnung der Route

function [sZ,lambdaZ,phiZ,KursZ,sU,lambdaU,phiU,KursU]=...
  RouteBerechnen(lambdaA,phiA,lambdaB,phiB,ds)

% Zwischenpunkte

dS=acos(cos(phiA)*cos(phiB)*cos(lambdaA-lambdaB)+sin(phiA)*sin(phiB));
[x,y,z]=sph2cart([lambdaA lambdaB], [phiA phiB], 1);
n=cross([x(1);y(1);z(1)],[x(2);y(2);z(2)]);
n0=n/norm(n);
v=[x(1);y(1);z(1)];
w=cross(n0,v);
sZ=unique([0:ds:dS dS]);
co=cos(sZ);
si=sin(sZ);
Z=v*co+w*si;
[lambdaZ phiZ rZ]=cart2sph(Z(1,:),Z(2,:),Z(3,:));

% Kurs in den Zwischenpunkten

coLam=cos(lambdaZ);
siLam=sin(lambdaZ);
coPhi=cos(phiZ);
siPhi=sin(phiZ);
Nord=[-siPhi.*coLam;-siPhi.*siLam;coPhi];
Ost=[-siLam;coLam;zeros(size(siLam))];
dZ=-v*si+w*co;
KursZ=atan2(dot(Ost,dZ),dot(Nord,dZ));
KursZ=KursZ+(KursZ<0)*2*pi;

% Umkehrpunkt

sU=atan(w(3)/v(3));
if sU<0
  sU=sU+2*pi;
end
if sU>dS
  sU=[];
  lambdaU=[];
  phiU=[];
  KursU=[];
else
  ZU=v*cos(sU)+w*sin(sU);
  [lambdaU phiU rU]=cart2sph(ZU(1,:),ZU(2,:),ZU(3,:));
  
  % Kurs im Umkehrpunkt
  
  coLam=cos(lambdaU);
  siLam=sin(lambdaU);
  coPhi=cos(phiU);
  siPhi=sin(phiU);
  Nord=[-siPhi.*coLam;-siPhi.*siLam;coPhi];
  Ost=[-siLam;coLam;zeros(size(siLam))];
  dU=-v*sin(sU)+w*cos(sU);
  KursU=atan2(dot(Ost,dU),dot(Nord,dU));
end

% Darstellung der Route

function [P1,P2,P3]=RouteDarstellen3D(lambdaA,phiA,lambdaB,phiB,...
   lambdaZ,phiZ,lambdaU,phiU)
[R1,R2,R3]=sph2cart([lambdaA lambdaB],[phiA phiB],1);
[P1,P2,P3]=sph2cart(lambdaZ,phiZ,1);
[Ux,Uy,Uz]=sph2cart(lambdaU,phiU,1);
plot3(P1,P2,P3,'k','erase','none','linewidth',1.5);
plot3([R1(1) 0 R1(2)],[R2(1) 0 R2(2)],[R3(1) 0 R3(2)],'--k','erase','none','linewidth',1.0);
plot3(Ux,Uy,Uz,'ok','erase','none');

% 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
