R=10;
z=real(zeros(201,201));
m=0;
n=0;
step = 0.1;
for x=-R:step:R
m = m + 1;
%x
for y=-sqrt(R*R - x*x):step:sqrt(R*R - x*x)
%y
n = int32(y / step) + R / step + 1;
%n
z(n, m) = real(sqrt(R*R - x*x - y*y));
endfor
endfor
%z
mesh(z);
另一种方法(from octave):
function [xx, yy, zz] = sphere (varargin)
[hax, varargin, nargin] = __plt_get_axis_arg__ ("sphere", varargin{:});
if (nargin > 1)
print_usage ();
elseif (nargin == 1)
n = varargin{1};
else
n = 20;
endif
theta = linspace (0, 2*pi, n+1);
phi = linspace (-pi/2, pi/2, n+1);
[theta,phi] = meshgrid (theta, phi);
x = cos (phi) .* cos (theta);
y = cos (phi) .* sin (theta);
z = sin (phi);
if (nargout > 0)
xx = x;
yy = y;
zz = z;
else
oldfig = [];
if (! isempty (hax))
oldfig = get (0, "currentfigure");
endif
unwind_protect
hax = newplot (hax);
surf (x, y, z);
unwind_protect_cleanup
if (! isempty (oldfig))
set (0, "currentfigure", oldfig);
endif
end_unwind_protect
endif
endfunction