How to make a surface in polar coordinates using polar3d?
Show older comments
Hello! I have an integral for a function which I plot in polar coordinates at a fixed polar angle theta (th). How to write a 3D plot for R in polar coordinates (angles a and th changes)? I made it using pol2cart, but this is not the best way to present the result. If somebody can help me to plot using polar3D or something like this, it would be great. Thank you.
clear all
s = 3;
n = 1;
r = 1;
t = 0.1;
th = 0:10:360; % angle theta
a = 0:1:360; % angle alpha
b = sqrt(2*n*t);
L = sqrt((4*t+r^2)/3);
fun = @(k,u,c,a) ((k.^2).*exp(-1.5*k.^2)).*((u.^2).*(1-u.^2).*exp(-(b.*u).^2).*(cos(s.*k.*u.*cos(a)/L))).*(((cos(c)).^2).*(cos(k.*sqrt(1-u.^2).*(s.*sin(a).*cos(th).*cos(c)+s.*sin(a).*sin(th).*sin(c))/(L))));
f3 = arrayfun(@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi),a);
B = ((6*sqrt(6)*(b^3))/(erf(b)*(pi^2)))*(1-(3/(2*b^2))*(1-((2*b*exp(-b^2))/(erf(b)*sqrt(pi)))))^(-1);
R = B*f3;
figure(3)
polar(a,R);
This is how it look in Cartesian coordinates (theta in radians), but it should look much better in polar3D.
2 Comments
Your code does not seem to be working here.
Note that there is no inbuilt command/function in MATLAB that produces a 3D polar plot.
As a workaround, you can convert polar to cartesian via pol2cart and use surf. However, if you need the plot in Polar 3D, check out this FileEx submission 3D Polar Plot
s = 3;
n = 1;
r = 1;
t = 0.1;
th = 0:10:360; % angle theta
a = 0:1:360; % angle alpha
b = sqrt(2*n*t);
L = sqrt((4*t+r^2)/3);
fun = @(k,u,c,a) ((k.^2).*exp(-1.5*k.^2)).*((u.^2).*(1-u.^2).*exp(-(b.*u).^2).*(cos(s.*k.*u.*cos(a)/L))).*(((cos(c)).^2).*(cos(k.*sqrt(1-u.^2).*(s.*sin(a).*cos(th).*cos(c)+s.*sin(a).*sin(th).*sin(c))/(L))));
f3 = arrayfun(@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi),a);
B = ((6*sqrt(6)*(b^3))/(erf(b)*(pi^2)))*(1-(3/(2*b^2))*(1-((2*b*exp(-b^2))/(erf(b)*sqrt(pi)))))^(-1);
R = B*f3;
figure(3)
polar(a,R);
Hexe
on 2 Jun 2023
Accepted Answer
Star Strider
on 2 Jun 2023
This takes too long to run here (it took 336.315050 seconds — 00:05:36.315049 — just now on MATLAB Online) however it plots the surface on a ‘synthetic’ polar coordinate axis in 3D. It would be relatively straightforward to change the Z axis coordinates (or ellimiinate them entirely). The original problem was that the code for the polar axes was not complete. It now works, as it worked in my earlier code.
Try this —
s = 3;
n = 1;
p = 0.1; % this is time
tv = 0:10:360; %3; % this is angle theta for polar rotation
r = 1;
a = 0:10:360;
a = deg2rad(a);
tv = deg2rad(tv);
tic
for k = 1:numel(tv)
t = tv(k)
b = sqrt(2*n*p);
L = sqrt((4*p+r^2)/3);
fun = @(k,u,c,a) ((k.^2).*exp(-1.5*k.^2)).*((u.^2).*(1-u.^2).*exp(-(b.*u).^2).*(cos(s.*k.*u.*cos(a)/L))).*(((cos(c)).^2).*(cos(k.*sqrt(1-u.^2).*(s.*sin(a).*cos(t).*cos(c)+s.*sin(a).*sin(t).*sin(c))/(L))));
f3 = arrayfun(@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi),a);
B = ((6*sqrt(6)*(b^3))/(erf(b)*(pi^2)))*(1-(3/(2*b^2))*(1-((2*b*exp(-b^2))/(erf(b)*sqrt(pi)))))^(-1);
R = B*f3;
ta = t*ones(size(tv));
[x,y,z] = pol2cart(a, R, ta);
X(k,:) = x;
Y(k,:) = y;
Z(k,:) = z;
toc % Total Time To This Point
LIT = toc/k % Mean Loop Iteration Time
end
toc
figure (2)
surfc(X, Y, Z)
colormap(turbo)
axis('equal')
axis('off')
hold on
a = deg2rad(0:2:360); % Create & Plot Polar Cylindrical Coordinates
r = 1;
xc = [0:0.25*r:r].'*cos(a);
yc = [0:0.25*r:r].'*sin(a);
xr = [0;r]*cos(0:pi/4:2*pi);
yr = [0;r]*sin(0:pi/4:2*pi);
plot3(xc.', yc.', zeros(size(xc)), '-k')
plot3(xr, yr, zeros(size(xr)), '-k')
zt = (0:1:max(tv)).';
zt1 = ones(size(zt));
surf(zt1*xc(end,:), zt1*yc(end,:), zt*ones(size(a)), 'FaceAlpha',0, 'MeshStyle','row')
hold off
text(xr(2,1:end-1)*1.1,yr(2,1:end-1)*1.1, zeros(1,size(xr,2)-1), compose('%3d°',(0:45:315)), 'Horiz','center','Vert','middle')
text(ones(size(zt))*xc(end,end-10), ones(size(zt))*yc(end,end-10), zt, compose('%.1f',zt))
toc
TTmin = toc/60
.
24 Comments
Hexe
on 2 Jun 2023
The long process is not a problem (in my MatLab it was 417 sec). But as a result, I got the same plot as before.
Star Strider
on 2 Jun 2023
This is what I get with the same code I posted —
I can’t figure out how to get it from MATLAB Drive where it’s saved (there appear to be no other options with MATLAB Online) to here, so I did a quick screencap instead.
.
Hexe
on 2 Jun 2023
Star Strider
on 2 Jun 2023
The turbo colormap was introduced in R2020b. If you have an earlier release, you will not have it, so the code stops there and does not draw the radial coordinates. Use a colormap available to you that produces the desired result.
Also, setting the last value of ‘r’ in the code to either 0.25 or 0.4 or 0.5 and using axis('square') instead of axis('equal') produces a better plot and keeps the radial coordinates circular rather than plotting them as ellipses.
Hexe
on 2 Jun 2023
Star Strider
on 2 Jun 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
Hexe
on 5 Jun 2023
Dear Star Strider!
I still have a question about presentation. The current code gives me the next plot (fig. 1). But actually this is not what I want. Is it possible someway to present this plot in polar3d, like in fig. 2?
Sincerely.
Star Strider
on 5 Jun 2023
I am not certain how figure 2 is created.
One option could be isosurface since figure 2 appears to be an isosurface plot. That sort of plot connects equal-valued ponts in a 3D volume to forma 3D surface. I am not certain how to get there from your code, since I do not understand what you are doing.
The approach in creating an isosurface plot is to create ‘X’, ‘Y’, and ‘Z’ volume matrices, and then define the isosurface using a function of all three matrices, and choosing the constant value that defines the surface.
Hexe
on 5 Jun 2023
That certainly works.
I added the polar coordinate axis to it. The polar coordiate position in the ‘z’ plane is defined by the zeros vector in the plot3 calls, so it is currently at the 0 level. To change its position, add any value (preferably between -2 and +2) to those vectors.
[theta,phi]=meshgrid(linspace(-pi/2,pi/2),linspace(0,2*pi));
r=1+sin(theta*3).*cos(phi*2);
X=cos(theta).*cos(phi).*r;
Y=cos(theta).*sin(phi).*r;
Z=sin(theta).*r;
figure
surf(X,Y,Z)
colormap(turbo)
axis('square')
axis('off')
hold on
tv = 0:10:360; %3; % this is angle theta for polar rotation
a = deg2rad(0:2:360); % Create & Plot Polar Cylindrical Coordinates
r = 2;
xc = [0:0.25*r:r].'*cos(a);
yc = [0:0.25*r:r].'*sin(a);
xr = [0;r]*cos(0:pi/4:2*pi);
yr = [0;r]*sin(0:pi/4:2*pi);
plot3(xc.', yc.', zeros(size(xc)), '-k')
plot3(xr, yr, zeros(size(xr)), '-k')
zt = (0:1:max(tv)).';
zt1 = ones(size(zt));
% surf(zt1*xc(end,:), zt1*yc(end,:), zt*ones(size(a)), 'FaceAlpha',0, 'MeshStyle','row')
% hold off
text(xr(2,1:end-1)*1.1,yr(2,1:end-1)*1.1, zeros(1,size(xr,2)-1), compose('%3d°',(0:45:315)), 'Horiz','center','Vert','middle')
% text(ones(size(zt))*xc(end,end-10), ones(size(zt))*yc(end,end-10), zt, compose('%.1f',zt))
.
Star Strider
on 6 Jun 2023
I am not certain what the problem is, however it works here.
My only suggestion is that it may be necessary to define a different ‘r’ variable, for example ‘ra’ for ‘r axis’ or ‘r2’ or something similar for the polar coordinate code creation, then change all other references to it where appropriate. Providing that the variable references in my code are consistent so that it works to create and plot the polar coordinates, the variable names can be anything that work as MATLAB variable names. There is nothing particularly special about them.
Hexe
on 6 Jun 2023
Star Strider
on 6 Jun 2023
The last code is throwing this error in line 24:
Arrays have incompatible sizes for this operation.
X=cos(a).*cos(tv).*r;
I have some errands to run, so I will deal with this later to understand what the problem is and how to correct it. (I have to run it offline because it requires more than the allotted 55 seconds to run here.)
.
Hexe
on 6 Jun 2023
Star Strider
on 6 Jun 2023
I added this assignment to the code just before the line that threw the error:
QS = [size(a); size(tv); size(r)]
and got this result:
QS =
100 100
100 100
1 37
I will spend some time on that later to see if I can determine what the problem is and how to fix it. It may just require making ‘r’ to be (1x100). If so, that simply requires a linspace call for ‘r’.
I would like to understand something that I seem to be missing. Does the code in this Comment relate to the code in this Comment, or are they entirely separate?
.
Hexe
on 7 Jun 2023
Star Strider
on 7 Jun 2023
Thank you!
I am going to experiment with ‘r’ later, and see if I can make the code work.
Hexe
on 7 Jun 2023
Star Strider
on 8 Jun 2023
I am not exactly certain what you want to do, however you can calculate them in spehrical coordinates, then use the sph2cart function to convert them to Cartesian coordinates that you can then plot, probably with surf. I have never done that (at least that I remember) because I have never needed to, however I will help with it as I can to get it to work.
I subscripted ‘R’ as ‘R(k,:)’ in my latest experiment, since since the original ‘R’ was (1x37), however that produces a (37x37) matrix, with continuing dimension problems for the ‘X’, ‘Y’, and ‘Z’ calculations. I do not have enough insight into your code to figure out how to make ‘R’ have compatible dimensions with ‘a’ and ‘tv’, so those calculations will work, and not also change ‘a’ and ‘tv’ at the same time, since ‘R’ is essentially a function of ‘a’. I thought that might be a straightforward fix. It is not straightforward at all.
Hexe
on 8 Jun 2023
Well, today I obtained the thing I wanted. If it is interesting for you, here is the code and the plot I got.
clear all
a_range = linspace(0, 2*pi, 100); % Adjust the number of points as needed
t_range = linspace(0, 2*pi, 100); % Adjust the number of points as needed
R = zeros(length(a_range), length(t_range));
progress_bar = waitbar(0, 'Calculating...');
s = 3;
n = 1;
p = 1;
r = 1;
b = sqrt(2*n*p);
L = sqrt((4*p+r^2)/3);
for i = 1:length(a_range)
for j = 1:length(t_range)
a = a_range(i);
t = t_range(j);
integrand = @(k, u, c) ((k.^2).*exp(-1.5*k.^2)).*((u.^2).*(1-u.^2).*exp(-(b.*u).^2).*(cos(s.*k.*u.*cos(a)/L))).*(((cos(c)).^2).*(cos(k.*sqrt(1-u.^2).*(s.*sin(a).*cos(t).*cos(c)+s.*sin(a).*sin(t).*sin(c))/(L)))); % Replace with your actual integrand function
result = integral3(integrand, 0, Inf, -1, 1, 0, 2*pi);
B = ((6*sqrt(6)*(b^3))/(erf(b)*(pi^2)))*(1-(3/(2*b^2))*(1-((2*b*exp(-b^2))/(erf(b)*sqrt(pi)))))^(-1);
R(i, j) = result*B;
end
waitbar(i / length(a_range), progress_bar);
end
close(progress_bar);
[A, T] = meshgrid(a_range, t_range);
Rho = R;
X = Rho .* cos(T) .* sin(A);
Y = Rho .* sin(T) .* sin(A);
Z = Rho .* cos(A);
figure;
surf(X, Y, Z);
xlabel('X');
ylabel('Y');
zlabel('Z');
title('3D Plot in Spherical Coordinates');
Star Strider
on 8 Jun 2023
It is interesting, because I never would have guessed what it was supposed to look like. I see that you used your own version of sph2cart to produce the ‘X’, ‘Y’, and ‘Z’ matrices.
Out of curiosity, what is it and does it represent?
Hexe
on 9 Jun 2023
Star Strider
on 9 Jun 2023
Interesting!
Thank you for the follow-up.
More Answers (0)
Categories
Find more on Polar Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
