I want to visualize a circle of intersection of a sphere by a plane
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Atom
am 22 Sep. 2023
Beantwortet: David Lovell
am 5 Dez. 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 Kommentare
Akzeptierte Antwort
Angelo Yeo
am 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 Kommentare
Angelo Yeo
am 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])
Weitere Antworten (1)
David Lovell
am 5 Dez. 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 Kommentare
Siehe auch
Kategorien
Mehr zu Surface and Mesh Plots finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!