how to plot a function on a 3D sphere?

9 views (last 30 days)
Wiqas Ahmad
Wiqas Ahmad on 3 Jan 2022
Commented: Wiqas Ahmad on 4 Jan 2022
N0=100;
r_med=[0.445 0.889 1.445 2];
sigma_g=7;
N_ang=91;
theta = (0:1:360)/180*pi;
phi = (0:1:360)/180*pi;
[x,y]=meshgrid(theta,phi);
Num_r = 50e3;
r = linspace(1,50,Num_r)./2;
I0=1;Q0=0;U0=0;V0=0;
for i = 1:length(r_med)
% [P11(i,:),P12(i,:),P33(i,:),P34(i,:),Qsca_c(i,:),~,~] = ZK_W_Cloud_PhaseFunc(N0,r_med(i),sigma_g,N_ang);
P11_t(i,:) = [P11(i,:) fliplr(P11(i,2:end))];
P12_t(i,:) = [P12(i,:) fliplr(P11(i,2:end))];
P33_t(i,:) = [P33(i,:) fliplr(P11(i,2:end))];
P34_t(i,:) = [P34(i,:) fliplr(P11(i,2:end))];
P1=repmat(P11_t(i,:),length(phi),1);
P2=repmat(P12_t(i,:),length(phi),1);
P3=repmat(P33_t(i,:),length(phi),1);
P4=repmat(P34_t(i,:),length(phi),1);
z = P1.*I0+((P2.*Q0.*cos(2*y))+(P2.*U0.*sin(2*y)));
end
[v,u,w]=sph2cart(x,y,z);
v=z.*cos(y).*cos(x);
u=z.*sin(y).*cos(x);
w=z.*sin(y);
figure,
surf(v,u,w),xlabel('x'),ylabel('y'),zlabel('z')
%set(gca,'zscale','log');
I want to obtain the figure as shown
However, the output of my code is just a circle
Please correct me.
  2 Comments
Wiqas Ahmad
Wiqas Ahmad on 3 Jan 2022
I had some corrction with the help of my colleague and now this code almost approaches my desire program. The range of theta and phi have been changed so that they run now for the whole values making a complete circle instead of half as shown in the figure. Similalry, the elements P11 and P12 elements have been flipped so that they also became vectors of 1*361 instead of 1*181. (the elements P33 and P34 can be ignored for instance). Now all of my data are fulfilled to plot them to obtain the desire figures. I have tried using the codes used for sphere but got the circle.

Sign in to comment.

Answers (1)

Matt J
Matt J on 3 Jan 2022
Edited: Matt J on 3 Jan 2022
My recommendation would be to build these surfaces using the cylinder function, e.g.,
R=1; %radius of circular cross-section
d=1.5; %radius of torus
fn=@(z) 2*R*(z-0.5);
x=fn( linspace(0,1) );
r1=d-sqrt(R^2-x.^2); %inner half surface
r2=d+sqrt(R^2-x.^2); %outer half surface
[X1,Y1,Z1]=cylinder(r1,40);
[X2,Y2,Z2]=cylinder(r2,40);
h(1)=surf(X1, Y1, fn(Z1)); hold on
h(2)=surf(X2, Y2, fn(Z2)); hold off
rotate(h,[1,0,0],90)
axis equal
xlabel X; ylabel Y; zlabel Z;
view(35,30)
  5 Comments
Matt J
Matt J on 4 Jan 2022
Yes, but that doesn't explain why my posted solution is insufficient.
Wiqas Ahmad
Wiqas Ahmad on 4 Jan 2022
I used it as follows and gives me the same output as yours
N0=100;
r_med=[0.445];
sigma_g=7;
N_ang=91;
theta0 = (0:1:360)/180*pi;
phi0 = (0:1:360)/180*pi;
[theta,phi]=meshgrid(theta0,phi0);
Num_r = 50e3;
r = linspace(1,50,Num_r)./2;
I0=1;Q0=1;U0=0;V0=0;
for i = 1:length(r_med)
%[P11(i,:),P12(i,:),P33(i,:),P34(i,:),Qsca_c(i,:),~,~] = ZK_W_Cloud_PhaseFunc(N0,r_med(i),sigma_g,N_ang);
P11_t(i,:) = [P11(i,:) fliplr(P11(i,2:end))];
P12_t(i,:) = [P12(i,:) fliplr(P11(i,2:end))];
P33_t(i,:) = [P33(i,:) fliplr(P11(i,2:end))];
P34_t(i,:) = [P34(i,:) fliplr(P11(i,2:end))];
P1=repmat(P11_t(i,:),length(phi0),1);
P2=repmat(P12_t(i,:),length(phi0),1);
P3=repmat(P33_t(i,:),length(phi0),1);
P4=repmat(P34_t(i,:),length(phi0),1);
phf = P1.*I0+((P2.*Q0.*cos(2*phi))+(P2.*U0.*sin(2*phi)));
end
[x,y,z]=pol2cart(theta,phi,phf);
x=z.*cos(theta);
y=z.*sin(theta);
z=z;
R=1; %radius of circular cross-section
d=1.5; %radius of torus
fn=@(z) 2.*R.*(z-0.5);
x=fn( linspace(0,1) );
r1=d-sqrt(R^2-x.^2); %inner half surface
r2=d+sqrt(R^2-x.^2); %outer half surface
[X1,Y1,Z1]=cylinder(r1,30);
[X2,Y2,Z2]=cylinder(r2,30);
figure,
h(1)=surf(abs(X1), abs(Y1), fn(abs(Z1))); hold on
h(2)=surf(abs(X2), abs(Y2), fn(abs(Z2))); hold off
rotate(h,[1,0,0],90)
axis equal
xlabel X; ylabel Y; zlabel Z;
view(35,30)
%mesh(x,y,z,abs(z)),xlabel('x'),ylabel('y'),zlabel('z')
%set(gca,'zscale','log');
The problem is that the output figure is independt of the values of I0 and Q0. By changing the light polarization from linear to circular, the figure doesn't change. I'm not clear how can I use it for my case?

Sign in to comment.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!