How to make a surface in polar coordinates using polar3d?

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]

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]

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]

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

.

[Hexe]

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]

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]

It is interesting, because I obtain the same figure as you got after changing the colormap from turbo to hsv.

[Star Strider]

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]

You are right. I work in R2018b and do not have turbo colormap. Thank you for advice with square axis - now the plot looks much better. It looks like, I will just reduce the axis z and cut the interval to show the most informative part of the plot.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Hexe]

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]

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]

Fig. 2 looks like the plot I want to obtain. In my code a is a polar angle, tv is a rotation angle and R obtained by integratio is a radius. The fig. 2 is created by this code:

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

surf(X,Y,Z)

[Star Strider]

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]

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)

[Star Strider]

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]

Honestly, I do not understand, how it works for you, because I obtain the same plot as before. Could you, please, show me now your plot with the last code?

[Star Strider]

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]

Thank you much for you help. I'll be waiting when you have time)

[Star Strider]

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]

Thank you. These codes are separately different and do not relate to each other. I just added the code (first Comment in your message) to show what kind of plot I want to obtain from my code (second comment in your message).

[Star Strider]

Thank you!

I am going to experiment with ‘r’ later, and see if I can make the code work.

[Hexe]

I just have got an idea, that it should be not polar3d, but 3d spherical coordinates. Is it possible to plot this code this way?

[Star Strider]

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]

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]

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]

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]

Interesting!

Thank you for the follow-up.