how to plot a function on a 3D sphere?

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.

[Walter Roberson]

How does this differ from https://www.mathworks.com/matlabcentral/answers/1620195-how-to-plot-a-3d-function?s_tid=srchtitle ?

[Wiqas Ahmad]

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.

Answer 1

[Matt J]

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)

[Wiqas Ahmad]

Thank you for your comment. But my scenario is quite different and I think it doesn't suit that. I have tried few ones before like this

[Matt J]

If your scenario is different, you will need to edit your question, so that people can understand how. The cylinder command works for both th surfaces you say you want ot plot.

[Wiqas Ahmad]

Let me explain to you in full detail. But the description involves physics terminologies, sorry for any inconvenience. This is the program for plotting the scattering phase function of a polarized light from a spherical water cloud droplet on a sphere. We used Mie-scattering theory to obtain the vectors P11,P12,P33 and P34 from a file "ZK_W_Cloud_PhaseFunc" which interpret the scattering properties of a spherical droplet. "N0,r_med(i),sigma_g,N_ang" these parameters represent the properties of droplet sizes. We use the following command

[P11(i,:),P12(i,:),P33(i,:),P34(i,:),Qsca_c(i,:),~,~] = ZK_W_Cloud_PhaseFunc(N0,r_med(i),sigma_g,N_ang);

For instance, consider i=1. Actually, P11 is 1*181 vector representing the properties of a half sphere, so we flipped to get the vector for the whole sphere (1*361) and similary for all others P12,P33,P34 as well. Note that all P's are the function of only theta (x-axis), so we did some manipulation to obtain all the vector P's for phi as well (y-axis). Now, the scattering phase function of a spherical water cloud droplet for polarized light can be obtained by using the command

z = P1.*I0+((P2.*Q0.*cos(2*y))+(P2.*U0.*sin(2*y)));

where

I0=1;Q0=1;U0=0;V0=0; 

is used for linearly polarized light, and

I0=1;Q0=0;U0=0;V0=1;

is used for circular polarized light. The function z is indepedent of the value of V0 and thus V0 can be ignored.

The first figure is obtained using linearly polarized light and the second one by using circularly polarized light. I hope everything will be clear now. Now I'm trying to obtain the same figures with the following code

N0=100;
r_med=[0.445];
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=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(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.*sin(y).*cos(x);
u=z.*sin(y).*sin(x);
w=z.*cos(y);
figure,
mesh(v,u,w),xlabel('x'),ylabel('y'),zlabel('z')
%set(gca,'zscale','log');

The output is

I have rotated the figure so that every view is clearly visible. I don't know how can I improve it further to get the desired figures.

[Matt J]

Yes, but that doesn't explain why my posted solution is insufficient.

[Wiqas Ahmad]

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?