function show(lat,lng,dst)

%
% function show(lat,lng,dst)
%
% display the current computational grid
% on a spherical topographic map
%
% stereographic projection to eye at lat,lng
% dst sets magnification, default = 4.5
%


global N topo pow cmap x y z;

if nargin<2 help show; return; end
if nargin<3 dst=4.5; end                  % fills window

title=get(gcf,'name');
aks=get(gcf,'child');                     % handle

if length(title)>3 & title(1:4)=='show'   % is it there 
   set(aks,'cameratarget',[0 0 0]);       % do these
   set(aks,'cameraviewangle',20);         % to gain control
   lat=pi*lat/180; lng=pi*lng/180;        % radians
   cam=dst*[cos(lat)*cos(lng) cos(lat)*sin(lng) sin(lat)];
   set(aks,'cameraposition',cam);         % easy - after all
   return;
end

disp('Generate 2-D contour map');         % draw contours on the sphere

atopo=topo(:,[360 1:360]);                     % don't ask
lg=(linspace(0,360,361)-0.5)*pi/180;; 
lt=(0+linspace(-89.5,89.5,180))*pi/180;
[C,H]=contour(lg,lt,atopo,[-4000 10 1000 3000],'k');  % pretty slow



view(3); disp('Move it to the sphere');

for k=1:length(H)                              % the contour lines
   lg=get(H(k),'xdata');                       % longitude
   lt=get(H(k),'ydata');                       % latitude
   zz=cos(lt);
   set(H(k),'xdata',zz.*cos(lg));
   set(H(k),'ydata',zz.*sin(lg));              % on unit sphere
   set(H(k),'zdata',sin(lt));
end
    
disp('Connect the nodes');             % connect nodes with black lines

r=1+2/5/N/N;
rx=r*x; ry=r*y; rz=r*z;
surface(rx,ry,rz,'facecolor','none');

M=N+1; nn=diag(1:M)*ones(M); nn=nn'-nn;     % constant on diags
x1=-flipud(rx(1:M,1:M));
y1=flipud(ry(1:M,1:M));
z1=flipud(rz(1:M,1:M));

for k=2-M:M-2                               % diagonal lines
    n=find(nn==k); 
    line(x1(n),y1(n),z1(n),'color',[0 0 0]);
    line(-y1(n),x1(n),z1(n),'color',[0 0 0]);
    line(-x1(n),-y1(n),z1(n),'color',[0 0 0]);

    line(y1(n),-x1(n),z1(n),'color',[0 0 0]);
end

disp('Fill in the colors');                   % fill in the colors

[xx,yy,zz]=sphere(90);
surface(-xx,-yy,zz,'facecolor','texture','cdata',topo,'edgecolor','none');

colormap(cmap); caxis([-7450 5750]);          % slight correction

axis equal; axis off;
set(gcf,'numbertitle','off');
set(gcf,'menu','none');
% set(gcf,'position',[775 500 500 500])
set(gcf,'name',sprintf('show  N=%d  and  pow = %g',N,pow));
axis([-1 1 -1 1 -1 1]);

show(lat,lng,dst);