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

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

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 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.
Alex Pedcenko on 27 Sep 2019
you're right

Alex Pedcenko on 5 Nov 2017
Edited: Alex Pedcenko on 27 Sep 2019
R=5; % outer radius of torus
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'); 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 ...
% 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 ...
% 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 on 5 Jul 2020
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
title('Not a Torus')
xlabel('X');ylabel('Y');zlabel('Z'); Stephen Cobeldick on 5 Jul 2020

### Tags

No tags entered yet.

### Community Treasure Hunt

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

Start Hunting!