
n = 20;   % plots portion of a 2n x n x nz grid in spherical shell
nz = 4;


clf
colormap([1 1 1; 1 .8 .8; .8 .8 1]);   
daspect([1 1 1])
axis off
hold on

% part of outer sphere:

[xc,yc] = meshgrid(linspace(-3,0,1.5*n+1), linspace(-1,1,n+1));
zc = ones(size(xc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,ones(size(zp)));   

% part of inner sphere:

%[xc,yc] = meshgrid(linspace(-3,1,2*n+1), linspace(-1,1,n+1));
[xc,yc] = meshgrid(linspace(-3,0,1.5*n+1), linspace(-1,1,n+1));
zc = zeros(size(xc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,2*ones(size(zp)));   

% mesh in z-direction:

[yc,zc] = meshgrid(linspace(-1,1,n+1),linspace(0,1,nz+1));
xc = ones(size(yc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,zeros(size(zp)));   

xc = zeros(size(yc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,zeros(size(zp)));   

[xc,zc] = meshgrid(linspace(0,1,.5*n+1),linspace(0,1,nz+1));
yc = ones(size(xc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,zeros(size(zp)));   

yc = -ones(size(xc));
[xp,yp,zp] = mapc2p(xc,yc,zc);
surf(xp,yp,zp,zeros(size(zp)));   

hold off
view([40 30])

return

