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!

[Stephen Cobeldick]

"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?

[Stephen Cobeldick]

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.

[Stephen Cobeldick]

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?

[Stephen Cobeldick]

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