function unitsphere(p)

if nargin < 1,
p = 2;
end

[X,Y] = meshgrid(0:.01:1,0:.01:1);
Z = (1 - X.^p - Y.^p).^(1/p);
ind = find(imag(Z)>eps | real(Z)<0);
Z(ind) = NaN;
surf([X X X X -X -X -X -X ],[Y Y -Y -Y Y Y -Y -Y],[Z -Z Z -Z Z -Z Z -Z]);
hold on
surf([X X X X -X -X -X -X ],[Z -Z Z -Z Z -Z Z -Z],[Y Y -Y -Y Y Y -Y -Y]);
surf([Z -Z Z -Z Z -Z Z -Z],[X X X X -X -X -X -X ],[Y Y -Y -Y Y Y -Y -Y]);
axis equal
shading interp
light('Position',5*[4,-9,8],'Style','local')
view(-16,18);
colormap spring

%%
set(gca,'XTick',-1:0.5:1)
set(gca,'YTick',-1:0.5:1)
set(gca,'ZTick',-1:0.5:1)
set(gcf,'Color',[1 1 1]);
set(gca,'FontSize',13)
title(['p = ' num2str(p) ])
hold off
shg