function KoordNBRio

% http://geometrie.hs-nb.de/MATLAB-Programme/KoordNBRio.m   Bearbeitungsstand: 6.3.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);
lambdaB=-dms2rad(43,13,22);
phiB   =-dms2rad(22,53,42);
[xA,yA,zA]=sph2cart(lambdaA,phiA,1);
[xB,yB,zB]=sph2cart(lambdaB,phiB,1);
grau={'color' [1 1 1]*0.3};
plot3([xA 0 xB],[yA 0 yB],[zA 0 zB],grau{:},'erase','none');
plot3([xA xB],[yA yB],[zA zB],'.',grau{:},'MarkerSize',20,'erase','none');
[xK,yK,zK]=sph2cart(lambdaA,[-90:90]*deg,1);
plot3(xK,yK,zK,grau{:},'erase','none')
[xK,yK,zK]=sph2cart([-180:180]*deg,phiA,1);
plot3(xK,yK,repmat(zK,size(xK)),grau{:},'erase','none')
[xK,yK,zK]=sph2cart(lambdaB,[-90:90]*deg,1);
plot3(xK,yK,zK,grau{:},'erase','none')
[xK,yK,zK]=sph2cart([-180:180]*deg,phiB,1);
plot3(xK,yK,repmat(zK,size(xK)),grau{:},'erase','none')
t=[-180:180]*deg;
plot3(cos(t),sin(t),zeros(size(t)),'k','erase','none')
t=[-90:90]*deg;
plot3(cos(t),zeros(size(t)),sin(t),'k','erase','none')
t=[0:35]*deg;
d=0.5;
plot3(cos(t),sin(t),zeros(size(t)),'k','erase','none','linewidth',1.5)
Spitze([cos(t(end));sin(t(end));0],[-sin(t(end));cos(t(end));0],d/4,45*deg);
Schrift([cos(t(end));sin(t(end));0],'^{ }\lambda','erase','none','HorizontalAlignment','center','VerticalAlignment','top')
plot3(cos(t),zeros(size(t)),sin(t),'k','erase','none','linewidth',1.5)
Spitze([cos(t(end));0;sin(t(end))],[-sin(t(end));0;cos(t(end))],d/4,45*deg);
Schrift([cos(t(end));0;sin(t(end))],'\phi_{ }','erase','none','HorizontalAlignment','right','VerticalAlignment','middle')
plot3([1+d;-1],[0;0],[0;0],'k','erase','none')
plot3([0;0],[1+d;-1],[0;0],'k','erase','none')
plot3([0;0],[0;0],[1+d;-1],'k','erase','none')
Spitze([1+d;0;0],[1;0;0],d/4,45*deg);
Spitze([0;1+d;0],[0;1;0],d/4,45*deg);
Spitze([0;0;1+d],[0;0;1],d/4,45*deg);
Schrift([cos(45*deg);sin(45*deg);0],'^{ }Äquator ','erase','none','HorizontalAlignment','left','VerticalAlignment','middle')
Schrift([0;0;-1],'^{  }Greenwichmeridian','erase','none','HorizontalAlignment','left','VerticalAlignment','top')
Schrift([1+d;0;0],'\itx','erase','none','HorizontalAlignment','left','VerticalAlignment','middle')
Schrift([0;1+d;0],'\ity','erase','none','HorizontalAlignment','left','VerticalAlignment','middle')
Schrift([0;0;1+d],'\itz  ','erase','none','HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([xA;yA;zA],'NB',grau{:},'erase','none','HorizontalAlignment','left','VerticalAlignment','bottom')
Schrift([xB;yB;zB],'Rio',grau{:},'erase','none','HorizontalAlignment','center','VerticalAlignment','top')
hold('off')
view(65,24);
axis('off') 
print('-depsc','KoordNBRio','-opengl')

% 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',24,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)],...
      'k','erase','none',varargin{:});
    end
  end

function t=degGrid(t1,t2,d);
if nargin<3
  d=1;
end
deg=pi/180;
delta=d*deg;
t=linspace(t1,t2,round((t2-t1)/delta));
