N = 50; thetavec = linspace(0,pi,N); phivec = linspace(0,2*pi,2*N); [th, ph] = meshgrid(thetavec,phivec); R = ones(size(th)); % should be your R(theta,phi) surface in general R = rand(size(th))*0.1 + 1; % should be your R(theta,phi) surface in general R = ones(size(th)) .* (sin(th*3)+1); % should be your R(theta,phi) surface in general x = R.*sin(th).*cos(ph); y = R.*sin(th).*sin(ph); z = R.*cos..