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

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.

Answer 1

[Angelo Yeo]

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')

[Atom]

Excellent... Thanks a lot. If I want to pass the axes through the centre (0,0,0), what I have to do.

[Angelo Yeo]

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])

[Atom]

Thank you very much.

Answer 2

[David Lovell]

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)