How to make a surface in polar coordinates using polar3d?

Question

Hexe 2023-6-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977014-how-to-make-a-surface-in-polar-coordinates-using-polar3d

评论： Star Strider 2023-6-9

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 个评论
显示无隐藏无

Dyuman Joshi 2023-6-2

编辑：Dyuman Joshi 2023-6-2

在 MATLAB Online 中打开

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);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

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 2023-6-2

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;

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('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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-6-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977014-how-to-make-a-surface-in-polar-coordinates-using-polar3d#answer_1249154

在 MATLAB Online 中打开

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 个评论
显示 22更早的评论隐藏 22更早的评论

Star Strider 2023-6-5

在 MATLAB Online 中打开

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))

.

Hexe 2023-6-6

编辑：Hexe 2023-6-6

Unfortunately, with my code, where polar r is R, which is calculated through the integral, it does not work, and the same plots are obtained as before. I did not think that in matlab this is such a problem. Apparently, I will have to try to build such a graph in another package. I think, that the code for a plot I want to obtain may be something like this, but it does not work "Matrix dimension must agree in X=cos(a).*cos(tv).*r;".

clear all

s = 3;

n = 1;

p = 0.1; % this is time

tv = 0:10:360; %3; % this is angle theta for polar rotation

ra = 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+ra^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).*sin(c))).*(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));

toc % Total Time To This Point

LIT = toc/k % Mean Loop Iteration Time

end

[a,tv]=meshgrid(linspace(-pi/2,pi/2),linspace(0,2*pi));

r=R;

X=cos(a).*cos(tv).*r;

Y=sin(a).*sin(tv).*r;

Z=cos(a).*r;

surf(X,Y,Z)

Hexe 2023-6-8

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');

Hexe 2023-6-9

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 2023-6-9

Interesting!

Thank you for the follow-up.

请先登录，再进行评论。

How to make a surface in polar coordinates using polar3d?

2 个评论
显示无隐藏无

采纳的回答

24 个评论
显示 22更早的评论隐藏 22更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to make a surface in polar coordinates using polar3d?

2 个评论 显示 无隐藏 无

采纳的回答

24 个评论 显示 22更早的评论隐藏 22更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

24 个评论
显示 22更早的评论隐藏 22更早的评论