Code covered by the BSD License

# geom3d

### David Legland (view profile)

19 Jun 2009 (Updated )

Library to handle 3D geometric primitives: create, intersect, display, and make basic computations

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

intersectPlaneSphere(plane, sphere)
```function circle = intersectPlaneSphere(plane, sphere)
%INTERSECTPLANESPHERE Return intersection circle between a plane and a sphere
%
%   CIRC = intersectPlaneSphere(PLANE, SPHERE)
%   Returns the circle which is the intersection of the given plane
%   and sphere.
%   PLANE  : [x0 y0 z0  dx1 dy1 dz1  dx2 dy2 dz2]
%   SPHERE : [XS YS ZS  RS]
%   CIRC   : [XC YC ZC  RC  THETA PHI PSI]
%   [x0 y0 z0] is the origin of the plane, [dx1 dy1 dz1] and [dx2 dy2 dz2]
%   are two direction vectors,
%   [XS YS ZS] are coordinates of the sphere center, RS is the sphere
%   [XC YC ZC] are coordinates of the circle center, RC is the radius of
%   the circle, [THETA PHI] is the normal of the plane containing the
%   circle (THETA being the colatitude, and PHI the azimut), and PSI is a
%   rotation angle around the normal (equal to zero in this function, but
%   kept for compatibility with other functions). All angles are given in
%   degrees.
%
%   planes3d, spheres, circles3d, intersectLinePlane, intersectLineSphere
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY
%   27/06/2007: change output format of circle, add support for multiple
%       data
%   2011-06-21 use degrees for angles

% number of inputs of each type
Ns = size(sphere, 1);
Np = size(plane, 1);

% unify data dimension
if Ns ~= Np
if Ns == 1
sphere = sphere(ones(Np, 1), :);
elseif Np == 1
plane = plane(ones(Ns, 1), :);
else
error('data should have same length, or one data should have length 1');
end
end
% center of the spheres
center  = sphere(:,1:3);

if size(sphere, 2) == 4
Rs  = sphere(:,4);
else
% assume default radius equal to 1
Rs  = ones(size(sphere, 1), 1);
end

% projection of sphere center on plane -> gives circle center
circle0 = projPointOnPlane(center, plane);

d   = distancePoints3d(center, circle0);
Rc  = sqrt(Rs.*Rs - d.*d);

% normal of planes = normal of circles
nor = planeNormal(plane);

% convert to angles
[theta, phi] = cart2sph2(nor(:,1), nor(:,2), nor(:,3));
psi = zeros(Np, 1);

% create structure for circle
k = 180 / pi;
circle = [circle0 Rc [theta phi psi]*k];
```