I want to visualize a circle of intersection of a sphere by a plane
18 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Atom
le 22 Sep 2023
Réponse apportée : David Lovell
le 5 Déc 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.
0 commentaires
Réponse acceptée
Angelo Yeo
le 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 commentaires
Angelo Yeo
le 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])
Plus de réponses (1)
David Lovell
le 5 Déc 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)
0 commentaires
Voir également
Catégories
En savoir plus sur Surface and Mesh Plots dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!