**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

19 views (last 30 days)

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

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

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;

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

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

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

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

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

.

Hexe
on 6 Jun 2023

Edited: Hexe
on 6 Jun 2023

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

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

.

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

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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)