%《MATLAB及其在理工課程中的應用指南》第五章例5－2－1程序
% 二次曲面參數對形狀的影響
% 西安電子科技大學出版社出版  陳懷琛編著 1999年10月初版，2004年11月第二版
%
echo off
for m=1:3
   a=input('a= （書上為a(1)=5,a(2)=5i,a(3)=5i) ');
   b=input('b= （書上為b(1)=4,b(2)=4,b(3)=4i) ');
   c=input('c= （書上為c(1)=3,c(2)=3,c(3)=3) ');	%輸入參數，N為網格線的數目
   d=input('d= （書上為d(1)=1,d(2)=1,d(3)=0.1) ');
   N=input('N= （書上為N(1)=20,N(2)=15,N(3)=10) ');
   xgrid=linspace(-abs(a),abs(a),N);		% 建立x網格坐標
   ygrid=linspace(-abs(b),abs(b),N);		% 建立y網格坐標
   [x,y]=meshgrid(xgrid,ygrid);		        % 建立N×N個點的xy網格坐標
   z=c*sqrt(d-y.*y/b/b-x.*x/a/a);u=1;		% u=1,表示z要取正負值
   z1=real(z);								% 取z的實部z1
   for k=2:N-1							    % 以下七行程序為取消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
   subplot(1,3,m),surf(x,y,z1),hold on		% 畫空間曲面
   set(gcf,'color','w')                     % 置圖形背景色為白色
   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
end