intersecdraw.m

matlab绘制平面与圆柱交面

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