# I want to visualize a circle of intersection of a sphere by a plane

14 views (last 30 days)
Atom on 22 Sep 2023
Answered: David Lovell on 5 Dec 2023
I want to visualize a circle of intersection of a sphere x^2+y^2+z^2=1, by a plane x+z=0.
N = 20;
thetavec = linspace(0,pi,N);
phivec = linspace(0,2*pi,2*N);
[th, ph] = meshgrid(thetavec,phivec);
R = ones(size(th));
x = R.*sin(th).*cos(ph);
y = R.*sin(th).*sin(ph);
z = R.*cos(th);
figure;
surf(x,y,z);
hold on;
[x z] = meshgrid(-2:0.1:2);
y = zeros(size(x, 1));
z=-x;
surf(x, y, z)
axis vis3d
% %
Please help me to plot the Fig such that one can visualize the axes x,y,z passes through (0,0,0) and circle of intersection is clearly visible. My plane plot process is not working. I am unable to make the sphere transparent, so that the plane along with the circle is clearly visible. Please help me.

Angelo Yeo on 22 Sep 2023
N = 20;
my_sphere = struct('x', [], 'y', [], 'z', []);
my_plane = struct('x', [], 'y', [], 'z', []);
[my_sphere.x, my_sphere.y, my_sphere.z] = sphere(N);
%{
intersection
x.^2 + y.^2 + z.^2 = 1
x + z = 0
x.^2 + y.^2 + (-x).^2 = 1
2x^2 + y^2 = 1
y^2 = 1-2x.^2
y = \pm sqrt(1-2x.^2)
%}
x1 = linspace(-1/sqrt(2), 1/sqrt(2), 100);
y1 = sqrt(1-2*x1.^2);
z1 = -x1;
x2 = linspace(-1/sqrt(2), 1/sqrt(2), 100);
y2 = -sqrt(1-2*x2.^2);
z2 = -x2;
[my_plane.x, my_plane.y] = meshgrid(-2:0.1:2);
my_plane.z = -1 * my_plane.x;
figure;
mesh(my_sphere.x, my_sphere.y, my_sphere.z,'FaceAlpha',0.5)
axis equal;
hold on;
mesh(my_plane.x, my_plane.y, my_plane.z,'FaceAlpha',0.5)
plot3(x1, y1, z1,'r','linewidth',2)
plot3(x2, y2, z2,'r','linewidth',2)
plot3(0,0,0,'ro','linewidth',2)
view(30, 10)
xlabel('x')
ylabel('y')
zlabel('z')
##### 3 CommentsShow 1 older commentHide 1 older comment
Angelo Yeo on 22 Sep 2023
plot3(linspace(-2,2,100), zeros(1, 100), zeros(1,100), 'k')
plot3(zeros(1, 100), linspace(-2,2,100), zeros(1,100), 'k')
plot3(zeros(1, 100), zeros(1,100), linspace(-2,2,100), 'k')
xlim([-2, 2])
ylim([-2, 2])
zlim([-2, 2])
Atom on 22 Sep 2023
Thank you very much.

David Lovell on 5 Dec 2023
Your plane exhibits a very simple relationship between x and z, so this one is easy. For a more general plane ax+by+cz=d, you can exploit the stereographic projection. Under this projection, the circle that represents the intersection between the sphere and the plane is mapped onto a circle on the 2D projection plane with center (-a/(c-d), -b/(c-d)) and radius sqrt(c^2-d^2+a^2+b^2)/(c-d). You can generate a set of points that constitute this circle, then invert these across the stereographic projection, and you get the circle you're looking for. Here's an example that can show this intersection for any plane that intersects the sphere:
% define inverse of stereographic projection (u,v) --> (x,y,z)
syms x(u,v) y(u,v) z(u,v)
x(u,v) = 2*u./(u.^2 + v.^2 + 1);
y(u,v) = 2*v./(u.^2 + v.^2 + 1);
z(u,v) = (u.^2 + v.^2 - 1)/(u.^2 + v.^2 + 1);
% plot sphere
N = 20;
my_sphere = struct('x', [], 'y', [], 'z', []);
my_plane = struct('x', [], 'y', [], 'z', []);
[my_sphere.x, my_sphere.y, my_sphere.z] = sphere(N);
mesh(my_sphere.x, my_sphere.y, my_sphere.z,'FaceAlpha',0.5)
axis equal
hold on
% plot plane
a = 1;
b = 0;
c = 1;
d = 0;
[X Y] = meshgrid(-1:0.1:1); % Generate x and y data
Z = -1/c*(a*X + b*Y - d); % Solve for z data
surf(X,Y,Z)
% plot intersection
theta = linspace(0,2*pi,100); % 100 equally spaced angles on a circle
uproj = -a/(c-d)+radius*cos(theta); % coordinates of circle on projection plane
xs = x(uproj,vproj); % invert the projection back into 3-D space
ys = y(uproj,vproj);
zs = z(uproj,vproj);
plot3(xs,ys,zs,'r-','Linewidth',5) % this is the circle in 3-D space
view(30,10)