function Kugelbereich

% http://geometrie.hs-nb.de/MATLAB-Programme/Kugelbereich.m   Bearbeitungsstand: 20.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); 
axis('equal') 
colormap(topocolors(min(min(H)),max(max(H)))) 
shading('interp')
title('') 
xlabel('x') 
ylabel('y') 
zlabel('z') 
caxis([min(min(H)) max(max(H))])
view(135,28); 
axis([-0.9 0.9 -0.6 0.6 0.4 0.6])
hold('on') 
phi0=30*deg;
phi1=50*deg;
lambda0=10*deg;
lambda1=70*deg;
t=[lambda0/deg:d:lambda1/deg]*deg;
plot3(cos(phi0)*cos(t),cos(phi0)*sin(t),sin(phi0)*ones(size(t)),'k','linewidth',1.5);
plot3(cos(phi1)*cos(t),cos(phi1)*sin(t),sin(phi1)*ones(size(t)),'k','linewidth',1.5);
t=[lambda1/deg:d:(360+lambda0)/deg]*deg;
plot3(cos(phi0)*cos(t),cos(phi0)*sin(t),sin(phi0)*ones(size(t)),'k','linewidth',0.5);
plot3(cos(phi1)*cos(t),cos(phi1)*sin(t),sin(phi1)*ones(size(t)),'k','linewidth',0.5);
t=[phi0/deg:d:phi1/deg]*deg;
plot3(cos(t)*cos(lambda0),cos(t)*sin(lambda0),sin(t),'k','linewidth',1.5);
plot3(cos(t)*cos(lambda1),cos(t)*sin(lambda1),sin(t),'k','linewidth',1.5);
t=[-90:d:phi0/deg]*deg;
plot3(cos(t)*cos(lambda0),cos(t)*sin(lambda0),sin(t),'k','linewidth',0.5);
plot3(cos(t)*cos(lambda1),cos(t)*sin(lambda1),sin(t),'k','linewidth',0.5);
t=[phi1/deg:d:90]*deg;
plot3(cos(t)*cos(lambda0),cos(t)*sin(lambda0),sin(t),'k','linewidth',0.5);
plot3(cos(t)*cos(lambda1),cos(t)*sin(lambda1),sin(t),'k','linewidth',0.5);

hold('off')
h=15*deg;
[x,y,z]=sph2cart(lambda0-h,phi0,1);
text(x,y,z,'\phi_0','EraseMode','none','FontName','TimesNewRoman',...
  'FontSize',48,'HorizontalAlignment','center','VerticalAlignment','baseline');
[x,y,z]=sph2cart(lambda0-h,phi1,1);
text(x,y,z,'\phi_{1  }','EraseMode','none','FontName','TimesNewRoman',...
  'FontSize',48,'HorizontalAlignment','center','VerticalAlignment','baseline');
[x,y,z]=sph2cart(lambda0,phi0-h,1);
text(x,y,z,'\lambda_0','EraseMode','none','FontName','TimesNewRoman',...
  'FontSize',48,'HorizontalAlignment','left','VerticalAlignment','middle');
[x,y,z]=sph2cart(lambda1,phi0-h,1);
text(x,y,z,'\lambda_1','EraseMode','none','FontName','TimesNewRoman',...
  'FontSize',48,'HorizontalAlignment','left','VerticalAlignment','middle');
axis('off') 
print('-deps','Kugelbereich','-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
