function Astro

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

clc;
close('all')
global ax
deg=pi/180;
lamA=-(55+22/60+47/3600)*deg; phiA=(29+49/60+31/3600)*deg; alph=54*deg;
lamB= (19+46/60+26/3600)*deg; phiB=-(8+11/60+51/3600)*deg; beta=62*deg;
[xa,ya,za]=sph2cart(lamA,phiA,1);
[xb,yb,zb]=sph2cart(lamB,phiB,1);
a=[xa;ya;za]
b=[xb;yb;zb]
v=cross(a,b)
P0=[a';b';0 0 1]\[cos(alph);cos(beta);0]
p=[dot(v,v) 2*dot(v,P0) dot(P0,P0)-1]
u=roots(p)
X1=P0+u(1)*v
X2=P0+u(2)*v
[lam1,phi1,r1]=cart2sph(X1(1),X1(2),X1(3));
[lam1d,lam1m,lam1s]=rad2dms(lam1)
[phi1d,phi1m,phi1s]=rad2dms(phi1)
[lam2,phi2,r2]=cart2sph(X2(1),X2(2),X2(3));
[lam2d,lam2m,lam2s]=rad2dms(lam2)
[phi2d,phi2m,phi2s]=rad2dms(phi2)

% 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); 
ax=[-1 1 -1 1 -1 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))])
O=[0;0;0];
drawPolygon(@plot3,[O a],'--k','linewidth',2.0,'erase','none')
drawPolygon(@plot3,[O b],'--k','linewidth',2.0,'erase','none')
drawPolygon(@plot3,[O X1],'--k','linewidth',2.0,'erase','none')
Spitze(2*a,a,0.2,30*deg,'k');
Spitze(a,a,0.2,30*deg,'k','erase','none');
Spitze(2*b,b,0.2,30*deg,'k');
Spitze(b,b,0.2,30*deg,'k','erase','none');
Spitze(X1,X1,0.2,30*deg,'k','erase','none');
Spitze(X1+a,a,0.2,30*deg,'k','erase','none');
Spitze(X1+b,b,0.2,30*deg,'k','erase','none');
drawPolygon(@plot3,[a 2*a],'k','linewidth',2.0)
drawPolygon(@plot3,[b 2*b],'k','linewidth',2.0)
drawPolygon(@plot3,[X1 X1+a],'k','linewidth',2.0)
drawPolygon(@plot3,[X1 X1+b],'k','linewidth',2.0)
drawPolygon(@plot3,[X1 2*X1],'k','linewidth',2.0)
view(69,-14)
r=view*[X1-a;1];
Schrift((X1+a)/2+a,'   Fixstern \itA','rotation',90+atan2(r(2),r(1))/deg,...
  'HorizontalAlignment','center','VerticalAlignment','baseline');
r=view*[X1-b;1];
Schrift((X1+b)/2+b,'Fixstern \itB','rotation',90-atan2(r(2),r(1))/deg,...
  'HorizontalAlignment','center','VerticalAlignment','bottom');
Schrift([0;0;0],'\itO','HorizontalAlignment','center','VerticalAlignment','top','erase','none');
Schrift(0.5*a,'\bf\ita   ','HorizontalAlignment','center','VerticalAlignment','middle','erase','none');
Schrift(1.5*a,'\bf\ita   ','HorizontalAlignment','center','VerticalAlignment','middle','erase','none');
Schrift(X1+0.5*a,'\bf\ita   ','HorizontalAlignment','center','VerticalAlignment','middle','erase','none');
Schrift(0.5*b,'\bf\itb','HorizontalAlignment','center','VerticalAlignment','top','erase','none');
Schrift(1.5*b,'\bf\itb','HorizontalAlignment','center','VerticalAlignment','top','erase','none');
Schrift(X1+0.6*b,'\bf\itb','HorizontalAlignment','center','VerticalAlignment','top','erase','none');
Schrift(0.5*X1,'\bf\itx   ','HorizontalAlignment','center','VerticalAlignment','middle','erase','none');
Schrift(X1,'\itX','HorizontalAlignment','left','VerticalAlignment','top','erase','none');
t=[0:360]*deg;
drawPolygon(@plot3,[cos(t);sin(t);zeros(size(t))],'k')
t=[-90:90]*deg;
drawPolygon(@plot3,[cos(lamA)*cos(t);sin(lamA)*cos(t);sin(t)],'k')
drawPolygon(@plot3,[cos(lamB)*cos(t);sin(lamB)*cos(t);sin(t)],'k')
drawPolygon(@plot3,[cos(lam1)*cos(t);sin(lam1)*cos(t);sin(t)],'k')
r=0.4;
c=cross(cross(a,X1)/norm(cross(a,X1)),a);
c=c/norm(c);
t=degGrid(0,alph);
drawPolygon(@plot3,repmat(X1,size(t))+r*repmat(cos(t),3,1).*repmat(a,size(t)) ...
                                     +r*repmat(sin(t),3,1).*repmat(c,size(t)),'k','linewidth',2.0)
Schrift(X1+r*(a*cos(alph/2)+c*sin(alph/2)),'\alpha','HorizontalAlignment','center','VerticalAlignment','top','erase','none');
c=cross(cross(b,X1)/norm(cross(b,X1)),b);
c=c/norm(c);
t=degGrid(0,beta);
drawPolygon(@plot3,repmat(X1,size(t))+r*repmat(cos(t),3,1).*repmat(b,size(t)) ...
                                     +r*repmat(sin(t),3,1).*repmat(c,size(t)),'k','linewidth',2.0)
Schrift(X1+0.9*r*(b*cos(beta/2)+c*sin(beta/2)),'\beta^{   }','HorizontalAlignment','center','VerticalAlignment','middle','erase','none');
axis([-0.6 1.8 -1 0.9 0 1.1])
axis('off')
hold('off')
print('-deps','Astro');

% 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 drawPolygon(draw,P,varargin);
global ax
feval(draw,P(1,:),P(2,:),P(3,:),varargin{:})
if min(P(1,:))<ax(1)
  ax(1)=min(P(1,:));
end
if max(P(1,:))>ax(2)
  ax(2)=max(P(1,:));
end
if min(P(2,:))<ax(3)
  ax(3)=min(P(2,:));
end
if max(P(2,:))>ax(4)
  ax(4)=max(P(2,:));
end
if min(P(3,:))<ax(5)
  ax(5)=min(P(3,:));
end
if max(P(3,:))>ax(6)
  ax(6)=max(P(3,:));
end

function [d,m,s]=rad2dms(rad)
D=rad/pi*180;
d=fix(D);
M=(D-d)*60;
m=fix(M);
s=(M-m)*60;

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));
function Schrift(P,txt,varargin)
text(P(1),P(2),P(3),txt,'FontName','TimesNewRoman','FontSize',36,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
