16 views (last 30 days)

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.

darova
on 11 Sep 2019

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

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.