function KongruenteWelten

% http://geometrie.hs-nb.de/MATLAB-Programme/KongruenteWelten.m   Bearbeitungsstand: 4.2.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/

close('all')
clc
deg=pi/180;
dreh1=inline('[1 0 0;0 cos(t) -sin(t);0 sin(t) cos(t)]*P','P','t');
dreh3=inline('[cos(t) -sin(t) 0;sin(t) cos(t) 0;0 0 1]*P','P','t');

% 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=6;
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
[L,B]=meshgrid([-180:d:180]*deg,[-90:d:90]*deg);
[X,Y,Z]=sph2cart(L,B,1);
[X,Y,Z]=transformMesh(dreh3,X,Y,Z,-120*deg);
[X,Y,Z]=transformMesh(dreh1,X,Y,Z,15*deg);
surf(X,Y,Z,H);
hold('on')
surf(X,Y,-2.2-Z,H);

[X,Y,Z]=transformMesh(dreh3,X,Y,Z,-45*deg);
surf(X+2.1,Y,Z,H);
surf(X+2.1,Y,-2.2-Z,H);

[X,Y,Z]=transformMesh(dreh1,X,Y,Z,-60*deg);
surf(X+4.2,Y,Z,H);
surf(X+4.2,Y,-2.2-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(0,15);
zlim([-3.8 0.8])
hold('off')
axis('off')
print('-deps','KongruenteWelten','-opengl')

function [X1,Y1,Z1]=transformMesh(transform,X0,Y0,Z0,varargin)
XYZ=transform([X0(:)';Y0(:)';Z0(:)'],varargin{:});
X1=reshape(XYZ(1,:),size(X0));
Y1=reshape(XYZ(2,:),size(Y0));
Z1=reshape(XYZ(3,:),size(Z0));

% 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

