Home > atmlab > geodetic > circle_sphere_intersect.m

circle_sphere_intersect

PURPOSE ^

CIRCLE_SPHERE_INTERSECT returns the intersection points of a circle and sphere,

SYNOPSIS ^

function [y]=circle_sphere_intersect(p1,n1,r1,r2)

DESCRIPTION ^

CIRCLE_SPHERE_INTERSECT returns the intersection points of a circle and sphere,
        The sphere is assumed to be centered at 
             cartesian origin [0,0,0].

 FORMAT  y=circle_sphere_intersect(p1,n1,r1,r2)

  OUT        y   the intersection points in cartesian coordinates
                      y is nan if there are no intersections

  IN            p1 the circle center in cartesian coordinates
                  n1 the normal of the cricle plane
                  r1  the radius of the circle
                  r2 the radius of the sphere (the sphere is placed 
                  at 0,0,0)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

circle_sphere_intersect.m

SOURCE CODE ^

0001 %CIRCLE_SPHERE_INTERSECT returns the intersection points of a circle and sphere,
0002 %        The sphere is assumed to be centered at
0003 %             cartesian origin [0,0,0].
0004 %
0005 % FORMAT  y=circle_sphere_intersect(p1,n1,r1,r2)
0006 %
0007 %  OUT        y   the intersection points in cartesian coordinates
0008 %                      y is nan if there are no intersections
0009 %
0010 %  IN            p1 the circle center in cartesian coordinates
0011 %                  n1 the normal of the cricle plane
0012 %                  r1  the radius of the circle
0013 %                  r2 the radius of the sphere (the sphere is placed
0014 %                  at 0,0,0)
0015 
0016 %HISTORY: created be Bengt Rydberg 2011-11-22
0017 function [y]=circle_sphere_intersect(p1,n1,r1,r2)
0018 
0019 %check input
0020 rqre_datatype( p1, @istensor1 );                        %&%
0021 rqre_datatype( n1, @istensor1 );                        %&%
0022 rqre_datatype( r1, @istensor1 );                        %&%
0023 rqre_datatype( r2, @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(r2)~=1                            %&%
0034     error('r2 must have length 1')                        %&%
0035 end                                                %&%
0036 
0037 p1=p1';
0038 n1=n1';
0039 
0040 %rotate coordinates if necessary to simplify code
0041 ind1=min(find(p1));
0042 if ind1==1
0043    indv=[1,2,3];
0044 elseif ind1==2
0045     indv=[2,1,3];
0046 else
0047    indv=[3,1,2];
0048 end
0049 
0050 
0051 indv=[1,2,3];
0052 
0053 %equation of the plane is then
0054 %ax+by+cz=n1(1)*x+n1(2)*y+n1(3)*z=d2=n1*p1'
0055 d=n1*p1';
0056 a=n1(indv(1));
0057 b=n1(indv(2));
0058 c=n1(indv(3));
0059 xa=p1(indv(1));
0060 ya=p1(indv(2));
0061 za=p1(indv(3));
0062 
0063 %we have three equations
0064 %x^2+y^2+z^2=r2^2                     (1)
0065 %(x-xa)^2+(y-ya)^2+(z-za)^2=r1^2          (2)
0066 %a*x+b*y+c*z=d                            (3)
0067 
0068 %(1) in (2) =>
0069 %-2*x*xa-2*y*ya-2*z*za=r1^2-r2^2-ra^2=r0^2
0070 ra2=p1*p1';
0071 r02=r1^2-r2^2-ra2;
0072 %=> x=-(r0^2+2*y*ya+2*z*za)/(2*xa)
0073 %=>x=-r0^2/(2*xa)-y*ya/xa-z*za/xa=x0+x1*y+x2*z     (4)
0074 x0=-r02/(2*xa);
0075 x1=-ya/xa;
0076 x2=-za/xa;
0077 
0078 if (a*x1+b)~=0
0079     %(4) in (3) =>
0080     %y=(d-a*x0-z*(c+a*x2))/(a*x1+b)=y1+y2*z      (5)
0081     y1=(d-a*x0)/(a*x1+b);
0082     y2=-(c+a*x2)/(a*x1+b);
0083     %(5) in (4) =>
0084     %x=x01+x02*z                               (6)
0085     x01=x0+x1*y1;
0086     x02=x1*y2+x2;
0087     %(5) & (6) in 1=>
0088     %z1*z^2+z2*z+z3=0
0089     z0=1+x02^2+y2^2;
0090     z1=2*x01*x02+2*y1*y2;
0091     z2=x01^2+y1^2-r2^2;
0092     zs1=(-z1+sqrt(z1^2-4*z0*z2))/(2*z0);
0093     zs2=(-z1-sqrt(z1^2-4*z0*z2))/(2*z0);
0094     xs1=x01+x02*zs1; 
0095     xs2=x01+x02*zs2; 
0096     ys1=y1+y2*zs1;
0097     ys2=y1+y2*zs2;
0098 else
0099     %(4) in (3) =>
0100     %z=(d-a*x0-y*(b+a*x1))/(a*x2+c)=z1+z2*y     (5)
0101     z1=(d-a*x0)/(a*x2+c);
0102     z2=-(b+a*x1)/(a*x2+c);
0103     %(5) in (4) =>
0104     %x=x01+x02*y                              (6)
0105     x01=x0+x2*z1;
0106     x02=x2*z2+x1;
0107     %(5) & (6) in 1=>
0108     %y1*y^2+y2*y+y3=0
0109     y0=1+x02^2+z2^2;
0110     y1=2*x01*x02+2*z1*z2;
0111     y2=x01^2+z1^2-r2^2;
0112     ys1=(-y1+sqrt(y1^2-4*y0*y2))/(2*y0);
0113     ys2=(-y1-sqrt(y1^2-4*y0*y2))/(2*y0);
0114     xs1=x01+x02*ys1; 
0115     xs2=x01+x02*ys2; 
0116     zs1=z1+z2*ys1;
0117     zs2=z1+z2*ys2;
0118 end
0119 
0120 y=[xs1,ys1,zs1;xs2,ys2,zs2];
0121 y=y(:,indv);
0122 if ~isreal(y)
0123   y=nan;
0124 end
0125 
0126     
0127 
0128 
0129 
0130 
0131 
0132

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