How do I plot a toroid in MATLAB?

I would like to plot a toroid in MATLAB but MATLAB does not have a built in function to do this.

Answer 1

[MathWorks Support Team]

You will need to formulate the x, y, and z-coordinate matrices manually and then plot them using the SURF function.

The SURF and MESH functions accept only one set of x, y, and z-coordinates, but in a toroid, (x,y) ordered pairs can have two corresponding z-coordinates. Therefore, to plot a toroid in MATLAB, you will need to plot the top and bottom halves as two separate surfaces on the same plot. For example:

%%Create R and THETA data

theta = 0:pi/10:2*pi;

r = 2*pi:pi/20:3*pi;

[R,T] = meshgrid(r,theta);

%%Create top and bottom halves

Z_top = 2*sin(R);

Z_bottom = -2*sin(R);

%%Convert to Cartesian coordinates and plot

[X,Y,Z] = pol2cart(T,R,Z_top);

surf(X,Y,Z);

hold on;

[X,Y,Z] = pol2cart(T,R,Z_bottom);

surf(X,Y,Z);

axis equal

shading interp

[Alex Pedcenko]

This is not a torus, it made of two caps: top and bottom are made from sin curve, hence the vertical cross-section of this surface is not circular, so this is not torus!

[Stephen23]

"This is not a torus..."

That is correct: it is not a torus.

However it is a toroid, which is what the title and the answer state it is.

It might help to revise the difference between a toroid (what this question and answer are about) and a torus (which is what your comment is about), e.g. from Wikipedia:

• "a toroid is a surface of revolution with a hole in the middle, like a doughnut, forming a solid body... If the revolved figure is a circle, then the object is called a torus."
• "a torus (plural tori) is a surface of revolution generated by revolving a circle in three-dimensional space about an axis that is coplanar with the circle."

or from Wolfram Mathematics http://mathworld.wolfram.com/Toroid.html :

• "A surface of revolution obtained by rotating a closed plane curve about an axis parallel to the plane which does not intersect the curve. The simplest toroid is the torus."

The answer creates a toroid from sine curves, just as it states. It does not create a torus, nor does the answer state that it creates a torus.

[Alex Pedcenko]

you're right

Answer 2

[Alex Pedcenko]

R=5; % outer radius of torus
r=2; % inner tube radius

th=linspace(0,2*pi,36); % e.g. 36 partitions along perimeter of the tube 
phi=linspace(0,2*pi,18); % e.g. 18 partitions along azimuth of torus

% we convert our vectors phi and th to [n x n] matrices with meshgrid command:
[Phi,Th]=meshgrid(phi,th); 
% now we generate n x n matrices for x,y,z according to eqn of torus
x=(R+r.*cos(Th)).*cos(Phi);
y=(R+r.*cos(Th)).*sin(Phi);
z=r.*sin(Th);
surf(x,y,z); % plot surface

daspect([1 1 1]) % preserves the shape of torus
colormap('jet') % change color appearance 
title('Torus')
xlabel('X');ylabel('Y');zlabel('Z');

[Sonal Gupta]

How do I further add spheres (to represent nodes of a torus network) on the edges of the plotted torus?

[Stephen23]

Alex Pedcenko's "Answer" moved here:

you can just add:

hold on
plot3(x,y,z,'o')

by the end of my code, if this is what you meant.

[Stephen23]

Sonal Gupta 's comment moved here:

Thanks. This works as needed. How do I make these markers filled? 'MarkerFaceColor' and 'filled' options don't seem to work. Also, is there a way to plot these markers as 3D color balls like this?

[Stephen23]

Alex Pedcenko's "Answer" moved here:

Well... you can make them 3D if you generate a set of spheres (see: help sphere) with a small radius and position them at these coordinates, but it will be much slower (your PC will need to render huge number of 3D objects) than this method above.

e.g

R=5; % outer radius of torus
r=2; % inner tube radius
th=linspace(0,2*pi,18); % e.g. 18 partitions along perimeter of the tube 
phi=linspace(0,2*pi,36); % e.g. 36 partitions along azimuth of torus
% we convert our vectors phi and th to [n x n] matrices with meshgrid command:
[Phi,Th]=meshgrid(phi,th); 
% now we generate n x n matrices for x,y,z according to eqn of torus
x=(R+r.*cos(Th)).*cos(Phi);
y=(R+r.*cos(Th)).*sin(Phi);
z=r.*sin(Th);
mesh(x,y,z,'EdgeColor','k'); % plot surface
daspect([1 1 1]) % preserves the shape of torus
%colormap('jet') % change color appearance 
title('Torus')
xlabel('X');ylabel('Y');zlabel('Z');
%% spheres:
hold on
[X Y Z]=sphere(10);
R=0.15; % size of the spheres
lighting gouraud
for i = 1:numel(x)
    surf(R*X+x(i),R*Y+y(i),R*Z+z(i),'FaceColor',[0.2 0.2 0.2],'FaceAlpha',0.8,'FaceLighting','gouraud','EdgeColor','none');
end
camlight

[Sonal Gupta]

Works perfectly. Thanks for the quick resolution.

[Sonal Gupta]

Hi Alex. Could you help me a bit more? I have been trying to plot lesser lines on the surface plot using your code. I need to plot a 3x16 (Th x Phi resp) network on a Torus shape drawn by your code. If I reduce the parameters to 3 and 16, the shape does not turn out to be torus like. I was trying this solution and have posted a related question on this forum as well, but got no reply. Can you answer the query posted on this forum?

[Alex Pedcenko]

Well,

you then need to plot torus itself with no lines with the same grid as before, e.g.

surf(x,y,z); % plot surface
shading interp % dont show lines of the "mesh"

But afterwards generate coodinates of your three (circular) lines in x-y plane, plot them on top of the torus surface and then plot spheres (if you still need them) at the nodes of these three lines.

[Johannes Stoerkle]

Thanks for the great examples. I used the code of Alex Pedcenko and extended it with the implicit description of Toroidal Surfaces with a conic constant, as commonly used in optical simulation programs and raytracing tools. Furthermore I computed the normal unity vectors and I considered positive and negative radii. The used equations are derived and described in https://johannes.stoerkle.info/wp-content/uploads/2020StoerkleJohannes_TechnicalNotes_TorodialSurface_MathFormulation_Optics.pdf and the results are validated with ZEMAX OpticStudio.

R = 4;   % outer radius of torus (set to 0, for cylinder)
r = 1;   % inner tube radius
k = 0;   % conic constant
sd = 10; % semidia in x
sdY = 10; % semidia in y

isRadiSwitched = r<0 && R>0 || r>0 && R<0;
isSmallR = abs(R)<abs(1/(1+k)*r);
isNeg_r = r<0;
isCyl = R==0;
if isNeg_r && ~isSmallR || ~isNeg_r && isSmallR
    signSqrt = -1;
else
    signSqrt = 1;
end

% x-y field spanned by polar coordinates
% -------------------------------------
th = linspace(-pi/2*0.99,pi/2*0.99,18); % partitions along perimeter of the tube 
phi = linspace(pi*0.01,pi*0.99,18); % partitions along azimuth of torus
% phi=linspace(pi*0.35,pi*0.65,18); % partitions along azimuth of torus
% we convert our vectors phi and th to [n x n] matrices with meshgrid command:
[Phi,Th] = meshgrid(phi,th); 
% now we generate n x n matrices for x,y,z according to eqn of torus
xO = (R-r+r.*cos(Th)).*cos(Phi);
yO = 1/sqrt(1+k)*r.*sin(Th);
% implicit surface description of torus
zO = R-sign(r)*sqrt((R+1/(1+k)*(signSqrt*sqrt(r^2-(1+k)*yO.^2)-r)).^2-...
    double(~isCyl)*xO.^2) ...
    +double(isSmallR)*(1/(1+k)*r-R)*2 ...
    -double(isRadiSwitched)*R*2;

% x-y field spanned by cartesian coordinates
% ----------------------------------------
testSemidias = (R+1/(1+k)*(signSqrt*sqrt(r^2-...
            (1+k)*sdY^2)-r))^2-double(~isCyl)*sd^2;
if 0>testSemidias || ~isreal(testSemidias) 
    % auto determine
    sdYAllowd = abs(r/sqrt(1+k)*0.99);
    % define according to the bound
    sdAllowd=abs(abs(R)+1/(1+k)*(signSqrt*sqrt(r^2-...
            (1+k)*sdYAllowd^2)-abs(r)))+double(isCyl)*1e9;
else
    sdAllowd = sd;
    sdYAllowd = sdY;
end
xRange = [-1 1]*min([sd sdAllowd]);
yRange = [-1 1]*min([sdY sdYAllowd]);
% create cartesian mesh
xLineSpace = linspace(xRange(1),xRange(2),18);
yLineSpace = linspace(yRange(1),yRange(2),18);
[x_,y_] = meshgrid(xLineSpace,yLineSpace); 
z_ = R-sign(r)*sqrt((R+1/(1+k)*(signSqrt*sqrt(r^2-(1+k)*y_.^2)-r)).^2-...
    double(~isCyl)*x_.^2)...
    +double(isSmallR)*(1/(1+k)*r-R)*2 ...
    -double(isRadiSwitched)*R*2;
% plot
% ----
figure; hold on
surf(yO,zO,xO,'faceColor','r'); % plot surface
surf(y_,z_,x_,'faceColor','g'); % plot surface
title(sprintf('Toridal: R=%0.0f, r=%0.0f, k=%0.1f, sd=%0.0f, sdY=%0.0f',...
    R,r,k,min([sd sdAllowd]),min([sdY sdYAllowd])))
xlabel('Y');ylabel('Z');zlabel('X');
grid on
axis equal
%% normal vector
% x_ = xO; y_ = yO; z_ = zO; % comment this out for normals of cartesian coord
% partial derivatives
dV_dx = sign(r)*double(~isCyl)*x_./...
    sqrt((R+1/(1+k)*(signSqrt*sqrt(r^2-(1+k)*y_.^2)-r)).^2-x_.^2);
dV_dy = sign(r)*y_.*(R+1/(1+k)*(signSqrt*sqrt(r^2-(1+k)*y_.^2)-r))...
    ./(sqrt((R+1/(1+k)*(signSqrt*sqrt(r^2-(1+k)*y_.^2)-r)).^2-...
    double(~isCyl)*x_.^2).*sqrt(r^2-(1+k)*y_.^2)*signSqrt);
dV_dz = repmat(-1,size(x_));
% create normal unity vectors
nx = dV_dx./sqrt(dV_dx.^2 + dV_dy.^2 + dV_dz.^2); 
ny = dV_dy./sqrt(dV_dx.^2 + dV_dy.^2 + dV_dz.^2); 
nz = dV_dz./sqrt(dV_dx.^2 + dV_dy.^2 + dV_dz.^2); 
% plot unity vectors
quiver3(y_,z_,x_,ny,nz,nx)
view(54,32)
axis equal

[Alex Pedcenko]

How about 3D spiral?

R=5; % outer radius of torus
a=1; % inner tube smaller radius
b=1; % inner tube larger radius
p=0.5; % pitch in z-direction

N=10; %turns along z

th=linspace(0,2*pi,36); % e.g. 36 partitions along perimeter of the tube 
phi=linspace(0,N*2*pi,36*N); % e.g. 18 partitions along azimuth of torus
% we convert our vectors phi and th to [n x n] matrices with meshgrid command:
[Phi,Th]=meshgrid(phi,th); 
% now we generate n x n matrices for x,y,z according to eqn of torus
x=(R+a.*cos(Th)).*cos(Phi);
y=(R+b.*cos(Th)).*sin(Phi);
z=a.*sin(Th)+p*Phi;
surf(x,y,z); % plot surface
daspect([1 1 1]) % preserves the shape of torus
colormap('jet') % change color appearance 
%shading interp
title('Not a Torus')
xlabel('X');ylabel('Y');zlabel('Z');

[Stephen23]

Looks good... please animate and upload a GIF!

For inspiration: https://www.mathworks.com/matlabcentral/fileexchange/25086-extrude-a-ribbon-tube-and-fly-through-it

Answer 3

[DGM]

MATLAB may not have a built-in function, but that doesn't mean there aren't any functions out there that can conveniently do the work.

https://www.mathworks.com/matlabcentral/fileexchange/58413-plot-superquadratic-surfaces

I'm sure this isn't the only thing on the File Exchange, but it's the one I use. Syntax is similar to sphere() or ellipsoid(), returning three matrices which can be fed to surf() or mesh(). The input arguments are the center location, radii, order, and number of points.

center = [0 0 0];
radius = [1 1 1 3];
order = 2;
npoints = 100;

[x y z] = supertoroid(center,radius,order,npoints);

surf(x,y,z)
shading flat
axis equal
colormap(parula)
view(-16,27)
camlight

As axis orders are independent and user-defined, the profile and sections do not have to be circular, but can be any superellipse:

radius = [1 1 2 3];
order = [5 3];

radius = [1 1 1 3];
order = [0.8 4];

Also included is a generalized superellipsoid tool.