xvec        = linspace(0,1,50);
yvec        = linspace(-1,1,50);
zvec        = yvec;
[X,Y,Z]     = meshgrid(xvec,yvec,zvec);
R           = (X.*X+Y.*Y+Z.*Z).^(0.5);
cosmat      = cos(2*pi*R*4)/4./R;
indx        = find(R>1);
cosmat(indx)= NaN*ones(1,length(indx));
indx        = find(R<0.01);
cosmat(indx)= NaN*ones(1,length(indx));
figure;

colormap(hot(100))
mapdata = colormap;
caxis([-1 1])

D  = cosmat;
p1 = patch(isosurface(X,Y,Z,R, 1),'FaceColor','red',...
    'EdgeColor','none');
p2 = patch(isocaps(X,Y,Z,D,-1000,'whichplane','xmin'),'FaceColor','interp',...
    'EdgeColor','none');
view(3);  daspect([1,1,1])
%axis([0 1 -1 1 -1 1])

camlight headlight;
camlight;lighting gouraud
camlight headlight;
colorbar
%isonormals(D,p1)



%p=patch(isosurface(X,Y,Z,cosmat,0));
%isonormals(X,Y,Z,cosmat,p)
%set(p,'FaceColor','red','EdgeColor','none');
%daspect([1 1 1])
%view(3)
%camlight
%lighting phong