MATLAB绘制二维等高线和二次曲面

clear;x=-3:0.1:3;y=-3:0.1:3;
[X,Y]=meshgrid(x,y);Z=sqrt(X.^2+Y.^2);
contour (X,Y,Z,10);
xlabel('X-axis'),ylabel('Y-axis');
title('Contour3 of Surface')
grid on

a=input('a='); b=input('b='); c=input('c=');
d=input('d='); N=input('N='); %输入参数，N为网格线数目
xgrid=linspace(-abs(a),abs(a),N); %建立x网格坐标
ygrid=linspace(-abs(b),abs(b),N); %建立y网格坐标
[x,y]=meshgrid(xgrid,ygrid); %确定N×N个点的x,y网格坐标
z=c*sqrt(d-y.*y/b^2-x.*x/a^2);u=1; %u=1,表示z要取正值
z1=real(z); %取z的实部z1
for k=2:N-1; %以下7行程序的作用是取消z中含虚数的点
for j=2:N-1
if imag(z(k,j))~=0 z1(k,j)=0;end
if all(imag(z([k-1:k+1],[j-1:j+1])))~=0 z1(k,j)=NaN;end
end
end
surf(x,y,z1),hold on %画空间曲面
if u==1 z2=-z1;surf (x,y,z2);%u=1时加画负半面
axis([-abs(a),abs(a),-abs(b),abs(b),-abs(c),abs(c)]);
end
xlabel('x'),ylabel('y'),zlabel('z')
hold off

运行程序，当a=5,b=4,c=3,d=1,N=50时结果为