# How to make a surface in polar coordinates using polar3d?

56 views (last 30 days)
Hexe on 1 Jun 2023
Commented: Star Strider on 9 Jun 2023
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.
Dyuman Joshi on 2 Jun 2023
Edited: Dyuman Joshi on 2 Jun 2023
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);
Arrays have incompatible sizes for this operation.

Error in solution>@(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)))) (line 9)
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))));

Error in solution>@(k,u,c)fun(k,u,c,a) (line 10)
f3 = arrayfun(@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi),a);

Error in integral3>@(y,z)fun(x(1)*ones(size(z)),y,z) (line 129)
@(y,z)fun(x(1)*ones(size(z)),y,z), ...

Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3>innerintegral (line 128)
Q1 = integral2Calc( ...

Error in integral3>@(x)innerintegral(x,fun,yminx,ymaxx,zminxy,zmaxxy,integral2options) (line 111)
f = @(x)innerintegral(x, fun, yminx, ymaxx, ...

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)

Error in integral3 (line 113)
Q = integralCalc(f,xmin,xmax,integralOptions);

Error in solution>@(a)integral3(@(k,u,c)fun(k,u,c,a),0,Inf,-1,1,0,2*pi) (line 10)
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
It works at fixed th as I wrote, for example, when th = 30. In this case I obtain the 2D plot. When I want to make th a variable to make surface, I can do it by pol2cart and surf. This way I obtained the plot attached (see the code below). But this is not what I want. I need a 3D plot exactly in polar coordinates, not in cartesian.
%Cut: sx = s*sin(a); sy = 0; sz = s*cos(a)
clear all
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;
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('on')
hold on
xc = [0:0.25*r:r];
yc = [0:0.25*r:r];
zc = [0:0.25*r:r];
xr = [0;r]*cos(0:pi/4:2*pi);
yr = [0;r]*sin(0:pi/4:2*pi);
zr = [0;r];
plot3(xc.', yc.', zc.', zeros(size(xc)), '-k')
plot3(xr, yr, zr, zeros(size(xr)), '-k')
zt = (0:1:max(tv)).';
zt1 = ones(size(zt));
surfc(zt1*xc(end,:), zt1*yc(end,:), zt*ones(size(tv)), '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

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;
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
.
Hexe on 9 Jun 2023
Yes, I had doubts on how it should be looked like, therefore made the same calculation in Mathematica and obtained there similar results, but the accuracy of the last package is not so good as in Matlab. Now I try to reproduce this figure also in Python but have a little bit problem with this task. About the sense there are many details but in general this is the correlation of internal electric fields in the material and how they interact.
Star Strider on 9 Jun 2023
Interesting!
Thank you for the follow-up.