function IntersecDraw1(A,B,C,D)
%平面方程AX+BY+CZ+D=0
%圓柱面方程　X＊X＋Y＊Y＝16　0＜＝Z＜＝12
%本函數用來求平面與圓柱面的關系。

Costhita=abs(C)/sqrt(A^2+B^2+C^2);  %曲面與xOy面的夾角的余弦
a=4/Costhita;       %橢圓長半軸，短半軸固定為4。
theta=0:0.01:2*pi;%角度循環
xx=a*cos(theta);%x
yy=4*sin(theta);%y
[m,n]=size(xx);
if C==0 %輸入平面與xoy面垂直
    if abs(D)/sqrt(A^2+B^2)>4       %外
        disp('無相交部分')
    elseif abs(D)/sqrt(A^2+B^2)==4  %切
        disp('相交部分為直線')        
    else
        disp('相交部分為矩形') %內
      
        b=sqrt(16-D^2/(A^2+B^2));%矩形的半寬(長固定為12)

        x=[b,b,-b,-b];  %繪制矩形x四個點
        y=[0,12,12,0];  %繪制矩形y四個點
        
        strTitle = '矩形,其';
        strTitle=strcat(strTitle,'x范圍為(');
        strTitle=strcat(strTitle,num2str(-b));
        strTitle=strcat(strTitle,',');
        strTitle=strcat(strTitle,num2str(b));
        strTitle=strcat(strTitle,')');        
        strTitle=strcat(strTitle,';y范圍為(0,12)');
        disp(strTitle);
        
        fill(x,y,'b');  %繪制矩形
        
               
        
%         x=[b,b,-b,-b b];
%         y=[0,12,12,0 0];
%         plot(x,y,'b');  %繪制矩形
    end

else % 輸入平面不與xoy面垂直
%過原點作平面的垂線,考慮垂線與z軸構成的平面與輸入平面的交線
    K=sqrt(A^2+B^2)/abs(C);%交線的斜率   
        if K==0%斜面與xoy面平行
            if -D/C<0 || -D/C>12 %平面與圓柱面相離　
            disp('無相交部分')
            else%平面切或割圓柱面
            disp('相交部分為圓')
            ezplot('x^2+y^2-16',[-5,5]);    %圓，x^2+y^2-16=0
            disp('圓的方程為：x^2+y^2=16');
            axis equal;	%坐標軸xy尺度相等
            end
        elseif K>0 && K<=1.5%斜率位于此區間時包含完整橢圓
            if -D/C==4*K+12 || -D/C==-4*K%與圓柱面相切時
                disp('相交部分為點')
            elseif -D/C>=-4*K+12 && -D/C<4*K+12%平面割圓柱面的上平面
                disp('相交部分為半橢圓')                
                x0=(12+D/C)/K/Costhita;%平面與上平面的割線在以橢圓中心為原點長短軸為坐標軸的坐標系下的方程。
                
                strTitle = '部分橢圓,長軸為';
                strTitle=strcat(strTitle,num2str(a));
                strTitle=strcat(strTitle,'，短軸為4,x范圍為(');
                strTitle=strcat(strTitle,num2str(-a));
                strTitle=strcat(strTitle,',');
                strTitle=strcat(strTitle,num2str(x0));
                strTitle=strcat(strTitle,')');
                disp(strTitle);
                
                count = sum(xx<=x0)+1;%半橢圓部分
              
                xxx=zeros(count);%以0初始化值
                yyy=zeros(count);
                j=1;
                for i=1:n
                    if(xx(i)<=x0)	%繪圖區間-a，x0
                        xxx(j)=xx(i);%半橢圓部分
                        yyy(j)=yy(i);
                        j=j+1;
                    end
                end
                xxx(count)=xxx(1);	%x首末點相連
                yyy(count)=yyy(1); 
                plot(xxx,yyy);  %繪橢圓，x^2/a^2+y^2/16=1            
                axis equal;%坐標系xy坐標軸尺度相等
                
                
                                
            %elseif -D/C>＝4*K && -D/C<＝-4*K+12%平面與圓柱面上下平面無交線
            elseif -D/C>=4*K && -D/C<=-4*K+12
                disp('相交部分為完整橢圓') 
                strTitle = '完整橢圓,長軸為';
                strTitle=strcat(strTitle,num2str(a));
                strTitle=strcat(strTitle,'，短軸為4');
                disp(strTitle);
                
                plot(xx,yy);
                axis equal;       %以上四行為繪制完整橢圓 
                
                
                
               % ezplot('x^2/a^2+y^2/16-1')%[-a:.001:a]);
               
            elseif -D/C>-4*K && -D/C<4*K%平面與圓柱面下平面相交
                disp('相交部分為半橢圓')               
                x0=D/C/K;%平面與下平面的割線在以橢圓中心為原點長短軸為坐標軸的坐標系下的方程。
                
                strTitle = '部分橢圓,長軸為';%標題，標示信息
                strTitle=strcat(strTitle,num2str(a));%字符串相連
                strTitle=strcat(strTitle,'，短軸為4,x范圍為(');
                strTitle=strcat(strTitle,num2str(x0));
                strTitle=strcat(strTitle,',');
                strTitle=strcat(strTitle,num2str(a));
                strTitle=strcat(strTitle,')');
                disp(strTitle);%顯示標題
                
                count = sum(xx>=x0)+1;
               
                xxx=zeros(count);%初始化
                yyy=zeros(count);
                j=1;
                for i=1:n
                    if(xx(i)>=x0)%在區間[x0,a]上繪制橢圓
                        xxx(j)=xx(i);%半橢圓部分
                        yyy(j)=yy(i);
                        j=j+1;
                    end
                end
                xxx(count)=xxx(1);%首末點相連
                yyy(count)=yyy(1);               
                plot(xxx,yyy); 
                
                          
             
            else %平面與圓柱面相離
                disp('無相交部分')
            end
        
        elseif K>1.5%斜率位于此部分時，相交部分不會出現完整橢圓
            if -D/C==4*K+12 || -D/C==-4*K%平面與圓柱面相切
                disp('相交部分為點')
            elseif -D/C>=4*K && -D/C<4*K+12%平面與圓柱面上平面相割
                disp('相交部分為半橢圓')
                %橢圓范圍為-a~x0
                x0=(12+D/C)/K/Costhita;%平面與上平面的割線在以橢圓中心為原點長短軸為坐標軸的坐標系下的方程。
                
                strTitle = '部分橢圓,長軸為';%標題，標示信息
                strTitle=strcat(strTitle,num2str(a));%字符串相連
                strTitle=strcat(strTitle,'，短軸為4,x范圍為(');
                strTitle=strcat(strTitle,num2str(-a));
                strTitle=strcat(strTitle,',');
                strTitle=strcat(strTitle,num2str(x0));
                strTitle=strcat(strTitle,')');
                disp(strTitle);%顯示標題
                
                count=sum(xx<=x0)+1;%半橢圓部分點計數
                xxx=zeros(count);%初始化值
                yyy=zeros(count);             
                j=1;
                for i=1:n
                    if (xx(i)<=x0)
                        xxx(j)=xx(i);%半橢圓部分
                        yyy(j)=yy(i);
                        j=j+1;
                    end
                end
                xxx(count)=xxx(1);%首末點相連
                yyy(count)=yyy(1);
                plot(xxx,yyy);%畫圖
                
            elseif -D/C>-4*K+12 && -D/C<4*K%平面與上下兩平面均相割　
                 disp('相交部分為半橢圓')  
                 
                x1=D/C/K/Costhita;%橢圓范圍為x1~x2    
                x2=(12+D/C)/K/Costhita;
                
                strTitle = '部分橢圓,長軸為';%標題，標示信息
                strTitle=strcat(strTitle,num2str(a));%字符串相連
                strTitle=strcat(strTitle,'，短軸為4,x范圍為(');
                strTitle=strcat(strTitle,num2str(x1));
                strTitle=strcat(strTitle,',');
                strTitle=strcat(strTitle,num2str(x2));
                strTitle=strcat(strTitle,')');
                disp(strTitle);%顯示標題
                
                [m,n]=size(xx);
                count = 0;%計數器初始化
                for i=1:n
                    if xx(i)>=x1&&xx(i)<=x2%在[x1,x2]上的點
                        count=count+1;
                    end
                end                      
                xxx=zeros(count);%初始化
                yyy=zeros(count);
                j=1;
                for i=1:n
                    if xx(i)>=x1&&xx(i)<=x2%橢圓范圍[x1,x2]
                        xxx(j)=xx(i);
                        yyy(j)=yy(i);
                        j=j+1;
                    end
                end
                plot(xxx,yyy);%畫圖
                                
            elseif -D/C>-4*K && -D/C<-4*K+12%平面只與下平面相割　
                disp('相交部分為半橢圓')
                
                strTitle = '部分橢圓,長軸為';%標題，標示信息
                strTitle=strcat(strTitle,num2str(a));%字符串相連
                strTitle=strcat(strTitle,'，短軸為4,x范圍為(');
                strTitle=strcat(strTitle,num2str(x0));
                strTitle=strcat(strTitle,',');
                strTitle=strcat(strTitle,num2str(a));
                strTitle=strcat(strTitle,')');
                disp(strTitle);%顯示標題
                
                x0=D/K/C/Costhita;%橢圓范圍為x0~a  
                count=sum(xx>=x0)+1;%在[x0,,a]上的點計數
                xxx=zeros(count);%初始化
                yyy=zeros(count);
                j=1;
                for i=1:n
                    if xx(i)>=x0;%橢圓范圍為x0~a 
                        xxx(j)=xx(i);
                        yyy(j)=yy(i);
                        j=j+1;
                    end
                end
                xxx(count)=xxx(1);%首末點相連
                yyy(count)=yyy(1);
                plot(xxx,yyy);%畫圖
                
            else %平面與圓柱面相離
                disp('無相交部分');
            end
        end
end