Clear Filters
Clear Filters

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

27 views (last 30 days)
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.

Accepted Answer

Angelo Yeo
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 Comments
Angelo Yeo
Angelo Yeo on 22 Sep 2023
Maybe adding such lines can be helpful.
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])

Sign in to comment.

More Answers (1)

David Lovell
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
radius = sqrt(c^2-d^2+a^2+b^2)/(c-d); % radius of circle on projection plane
uproj = -a/(c-d)+radius*cos(theta); % coordinates of circle on projection plane
vproj = -b/(c-d)+radius*sin(theta);
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)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!