3D plot in polar coordinates

6 ビュー (過去 30 日間)
Hexe
Hexe 2023 年 5 月 21 日
コメント済み: Star Strider 2023 年 5 月 23 日
Hello! I made a code for solving the integral and it looks realistic in polar coordinates. But, how to present it in 3D as a figure in volume? Perhaps I should add an aditional rotation angle for this. Maybe it is easy for those who know ow to do it. If somebody knows, please, help me. Thank you.
s = 3;
n = 1;
t = 0.1;
r = 1;
a = 0:1:360;
a = a*pi/180;
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(s.*k.*sqrt(1-u.^2).*sin(a).*(cos(c)+sin(c))/(L*sqrt(2)))));
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;
polar(a,R);

採用された回答

Star Strider
Star Strider 2023 年 5 月 22 日
This took a while to get working correctly, and takes about 500 seconds to run, so I will post the code here slthough not the plot —
s = 3;
n = 1;
tv = 0:0.1:2;%3;
r = 1;
a = 0:2:360;
a = deg2rad(a);
tic
for k = 1:numel(tv)
t = tv(k)
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(s.*k.*sqrt(1-u.^2).*sin(a).*(cos(c)+sin(c))/(L*sqrt(2)))));
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(a));
[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
surf(X, Y, Z)
colormap(turbo)
axis('equal')
axis('off')
hold on
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
figure
surf(X, Y, Z, 'EdgeColor','none')
colormap(turbo)
axis('equal')
axis('off')
hold on
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
% END
I only run it here from 0 to 2, because there is not much detail beyond that. The time dimension is the ‘Z’ dimension. Most of the detail is between 0 and 0.5.
.
  6 件のコメント
Hexe
Hexe 2023 年 5 月 23 日
編集済み: Hexe 2023 年 5 月 23 日
Dear Star Strider.
Your previous answer is good and acceptable. I just wanted to clarify about another way of presentation the result, but probably it is not possible.
Thank you a lot.
Star Strider
Star Strider 2023 年 5 月 23 日
As always, my pleasure!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeGraphics についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by