MATLAB Answers

How do I plot a toroid in MATLAB?

315 views (last 30 days)
I would like to plot a toroid in MATLAB but MATLAB does not have a built in function to do this.

Accepted Answer

MathWorks Support Team
MathWorks Support Team on 27 Jun 2009
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

  3 Comments

Alex Pedcenko
Alex Pedcenko on 15 Oct 2017
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
Stephen Cobeldick on 27 Sep 2019
"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:
or from Wolfram Mathematics http://mathworld.wolfram.com/Toroid.html :
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.

Sign in to comment.

More Answers (1)

Alex Pedcenko
Alex Pedcenko on 5 Nov 2017
Edited: Alex Pedcenko on 27 Sep 2019
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');
torus1.png

  10 Comments

Show 7 older comments
Johannes Stoerkle
Johannes Stoerkle on 26 Feb 2020
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
Alex Pedcenko on 5 Jul 2020
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');

Sign in to comment.

Tags

No tags entered yet.

Products

Community Treasure Hunt

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

Start Hunting!

Translated by