Home > atmlab > geodetic > circle_plane_intersect.m

circle_plane_intersect

PURPOSE ^

CIRCLE_PLANE returns the intesection points of a circle and plane in 3 dimensions

SYNOPSIS ^

function [y]=circle_plane_intersect(p1,n1,r1,p2,n2)

DESCRIPTION ^

CIRCLE_PLANE returns the intesection points of a circle and plane in 3 dimensions

 FORMAT y=circle_plane_intersect(p1,n1,r1,p2,n2)

 OUT      
         y     cartesian coordinates of the intersection points
                y is nan if there are no intersections
        
 IN     p1  the circle center in cartesian coordinates
          n1  the normal to the circle
          r1   the radius of the circle
          p2  a point in the plane 
          n2  the normal to the plane

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

circle_plane_intersect.m

SOURCE CODE ^

0001 %CIRCLE_PLANE returns the intesection points of a circle and plane in 3 dimensions
0002 %
0003 % FORMAT y=circle_plane_intersect(p1,n1,r1,p2,n2)
0004 %
0005 % OUT
0006 %         y     cartesian coordinates of the intersection points
0007 %                y is nan if there are no intersections
0008 %
0009 % IN     p1  the circle center in cartesian coordinates
0010 %          n1  the normal to the circle
0011 %          r1   the radius of the circle
0012 %          p2  a point in the plane
0013 %          n2  the normal to the plane
0014 
0015 % HISTORY: Created by Bengt Rydberg 2011-11-22
0016 function [y]=circle_plane_intersect(p1,n1,r1,p2,n2)
0017 
0018 %check input
0019 rqre_datatype( p1, @istensor1 );                        %&%
0020 rqre_datatype( n1, @istensor1 );                        %&%
0021 rqre_datatype( r1, @istensor1 );                        %&%
0022 rqre_datatype( p2, @istensor1 );                        %&%
0023 rqre_datatype( n2, @istensor1 );                        %&%
0024 if length(p1)~=3                            %&%
0025     error('p1 must have length 3')                        %&%
0026 end                                                %&%
0027 if length(n1)~=3                            %&%
0028     error('n1 must have length 3')                    %&%
0029 end                                                %&%
0030 if length(r1)~=1                            %&%
0031     error('r1 must have length 1')                        %&%
0032 end                                                %&%
0033 if length(p2)~=3                            %&%
0034     error('p2 must have length 3')                        %&%
0035 end                                                %&%
0036 if length(n2)~=3                            %&%
0037     error('n2 must have length 3')                    %&%
0038 end                                                %&%
0039 
0040 p1=p1';
0041 n1=n1';
0042 p2=p2';
0043 n2=n2';
0044 
0045 
0046 if sum(cross(n1,n2))==0
0047  %the circle plane and the other plane is parallel
0048 %and there are no crossings
0049  y=nan;
0050  return
0051 end
0052 
0053 %rotate the coordinate system to make the code simpler
0054 ind1=find(abs(n2)==max(abs(n2)));
0055 
0056 if ind1==1
0057    indv=[1,2,3];
0058 elseif ind1==2
0059     indv=[2,1,3];
0060 else
0061    indv=[3,1,2];
0062 end
0063 
0064 %equation of the plane is then
0065 %a2x+b2y+c2z=n2(1)*x+n2(2)*y+n2(3)*z=d2=n2*p2'
0066 d2=n2*p2';
0067 a2=n2(indv(1));
0068 b2=n2(indv(2));
0069 c2=n2(indv(3));
0070 
0071 %the circle plane is
0072 %a1x+b1x+c1x=d1=n1(1)*x+n1(2)*y+n1(3)*z=d1=n1*p1'
0073 d1=n1*p1';
0074 a1=n1(indv(1));
0075 b1=n1(indv(2));
0076 c1=n1(indv(3));
0077 
0078 %
0079 xa=p1(indv(1));
0080 ya=p1(indv(2));
0081 za=p1(indv(3));
0082 
0083 %x=(d2-b2*y-c2*z)/a2;
0084 b0=b1-b2*a1/a2;
0085 c0=c1-c2*a1/a2;
0086 d0=d1-d2*a1/a2;
0087 
0088 
0089 
0090 if b0~=0
0091         y1=d0/b0;
0092         y2=c0/b0;
0093         %y=y1-y2*z;
0094         %x=d2/a2*(1-b2/a2*y1)-z*(b2/a2*y2-c2/a2)
0095         %x=x1-x2*z
0096         x1=d2/a2-b2/a2*d0/b0;
0097         x2=c2/a2-b2/a2*y2;
0098         %z^2*z0+z*z1+z2=0
0099         z0=x2^2+y2^2+1;
0100         z1=2*(x2*(xa-x1)+y2*(ya-y1)-za);
0101         z2=-r1^2+(x1-xa)^2+(y1-ya)^2+za^2;
0102         %%%%%%%
0103         zs1=(-z1+sqrt(z1^2-4*z0*z2))/(2*z0);
0104         zs2=(-z1-sqrt(z1^2-4*z0*z2))/(2*z0);
0105         xs1=x1-x2*zs1;
0106         xs2= x1-x2*zs2;
0107         ys1=y1-y2*zs1;
0108         ys2= y1-y2*zs2;
0109         y=[xs1,ys1,zs1;xs2,ys2,zs2];
0110    else
0111         z1=d0/c0;
0112         z2=b0/c0;
0113         x1=d2/a2-c2/a2*d0/c0;
0114         x2=b2/a2-c2/a2*z2;
0115         y0=x2^2+z2^2+1;
0116         y1=2*(x2*(xa-x1)+z2*(za-z1)-ya);
0117         y2=-r1^2+(x1-xa)^2+(z1-za)^2+ya^2;
0118         ys1=(-y1+sqrt(y1^2-4*y0*y2))/(2*y0);
0119         ys2=(-y1-sqrt(y1^2-4*y0*y2))/(2*y0);
0120         xs1=x1-x2*ys1;
0121         xs2= x1-x2*ys2;
0122         zs1=z1-z2*ys1;
0123         zs2= z1-z2*ys2;
0124         y=[xs1,ys1,zs1;xs2,ys2,zs2];
0125 end
0126 
0127 y=y(:,indv);
0128 if ~isreal(y)
0129   y=nan;
0130 end
0131 
0132

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005