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.
[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);
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.
