Rotate Spherical Cap On Sphere
11 次查看(过去 30 天)
显示 更早的评论
Dear All,
I am trying to produce a plot where I have a sphere with another spherical section representation a region of interest. I also want to be able to rotate this spherical section over any theta,phi. I will walk though what I have done and explain where I am stuck.
The first think I do is produce a sphere that is transparent. Here is the code for that.
%%%%%%%%%%% Generate a sphere consisting of 200 by 200 faces %%%%%%%%%%%%%
[x,y,z]=sphere(200);
rad=1;
hSurface1=surf(rad*x,rad*y,rad*z);hold on
set(hSurface1,'FaceColor',[0 0 1],'FaceAlpha',0.25,'FaceLighting','gouraud','EdgeColor','none')
hSurface2=surf(rad*x*.7,rad*y*.7,rad*z*.7);hold on
set(hSurface2,'FaceColor',[0 1 0],'FaceAlpha',0.25,'FaceLighting','gouraud','EdgeColor','none')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This is one sphere with unit radius =1. The inner sphere represents the shell thickness. It is for visual purposes only. I next define a range for the spherical section. As expected the phi range is from 0-2pi. The theta range is determined by a region of interest for the problem. For example:
numPoints = 4;
thetaMin = 0
thetaMax = pi/2-0.68;
thetaRange = linspace(thetaMin, thetaMax, numPoints);
phiMin = 0;
phiMax = 2*pi;
phiRange = linspace(phiMin, phiMax, numPoints);
[theta,phi] = meshgrid(thetaRange,phiRange,numPoints);
[x1,y1,z1] = sph2cart(phi, theta, rad);
x1 = round(x1,4);
y1 = round(y1,4);
z1 = round(z1,4);
%surf(x1,y1,z1)
plot3(x1,y1,z1,'*w');
This will generate a spherical cap (or points) that sits on top of the sphere. My question is how do I rotate the cap (or indvidual points) using any theta and phi. Is there a unique way to transform the mesh points? I have tried to convert the individual points in the following approach:
theta_rot = 0;phi_rot = 0;
These are the new theta and phi I want to rotate mesh points to. I tried to find the correct difference in theta and phi to calculate a new vector b.
for i=1:length(x1)
for j=1:length(x1)
a=[x1(i,j) y1(i,j) z1(i,j)]
[azimuth,elevation,r] = cart2sph(x1(i,j),y1(i,j),z1(i,j))
b=[sind(theta_rot+elevation)*cosd(phi_rot+azimuth) sind(theta_rot+elevation)*sind(phi_rot+azimuth) cosd(theta_rot+elevation)]
w=cross(a,b);
d=dot(a,b);
Id = [1 0 0;0 1 0;0 0 1];
vx = [0 -w(3) w(2);w(3) 0 w(1);-w(2) w(1) 0];
R = Id + vx + vx^2*((1-d)/norm(w)^2);
rot_vals = mtimes(R,[x1(i,j) y1(i,j) z1(i,j)]');
rot_x1(i,j) = rot_vals(1);
rot_y1(i,j) = rot_vals(2);
rot_z1(i,j) = rot_vals(3);
end
end
This approach does not work since I am unable to find a better way to construct vector (b) by finding the difference in the points. I hope this is not too confusing and any help is appreaciated.
0 个评论
回答(1 个)
darova
2019-9-11
I rotated data about X axis first (theta angle). Then i am rotating about Z (phi angle)
clc,clear
p = linspace(0,2*pi,50);
t = linspace(pi/2,pi/6,10);
[phi,theta] = ndgrid(p,t);
[x1,y1,z1] = sph2cart(phi, theta, 1);
[X,Y,Z] = sphere(20);
h = surf(X,Y,Z);
set(h,'EdgeColor','none')
% set(h,'FaceCOlor', 'b');
alpha(h,0.5)
t = 45;
st1 = sind( t*sind(t) );
ct1 = cosd( t*sind(t) );
hold on
for t = linspace(0,180,100)
% ---------- PLAY WITH THIS
% st1 = sind( t*sind(t) );
% ct1 = cosd( t*sind(t) );
st3 = sind( t );
ct3 = cosd( t );
% ------------
Rx = [1 0 0; 0 ct1 -st1; 0 st1 ct1]; % rotation matrix around X axis
Rz = [ct3 -st3 0; st3 ct3 0; 0 0 1]; % rotation matrix around Z axis
V = Rx*[x1(:) y1(:) z1(:)]'; % rotate data
V = Rz*V;
x2 = reshape(V(1,:),size(x1));
y2 = reshape(V(2,:),size(x1));
z2 = reshape(V(3,:),size(x1));
h = plot3(x2,y2,z2,'g');
pause(0.05)
delete(h)
end
hold off
axis equal
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!