function main()  % illustration of a triple torus      % torus radii    r = 1;   R = 3;      Kb = R+r;   Ks = R-r;   % Km controls the smoothness of the transition from one ring to the others   Km = 0.5125*Kb;   L = 2.5*(r+R);   h = 0.2; % grid size      X = (-L):h:L; m = length(X);   Y = (-L):h:L; n = length(Y);   Z = (-1.1*r):h:(1.1*r); k = length(Z);   W = zeros(m, n, k);   for i=1:length(X)      for j=1:length(Y)         x = X(i);         y = Y(j);                  [x, y] = triple_torus_function (x, y, r, R, Kb, Km);         val = (sqrt(x^2+y^2)-R)^2-r^2;         W(i, j, :) = val + Z.^2;      end   end   figure(1); clf; hold on;   axis equal; axis off;   light_green=[184, 224, 98]/256; % light green   H = patch(isosurface(X, Y, Z, W, 0));   isonormals(X, Y, Z, W, H);   mycolor = light_green;   set(H, 'FaceColor', mycolor, 'EdgeColor','none', 'FaceAlpha', 1);   set(H, 'SpecularColorReflectance', 0.1, 'DiffuseStrength', 0.8);   set(H, 'FaceLighting', 'phong', 'AmbientStrength', 0.3);   set(H, 'SpecularExponent', 108);   daspect([1 1 1]);   axis tight;   colormap(prism(28))   view(-12, 40);   % add in a source of light   camlight (-50, 54); lighting phong;   print('-dpng', '-r400',  ...         sprintf('Triple_torus_illustration.png'));   function [x, y] = triple_torus_function (x, y, r, R, Kb, Km)% a deformation in the plane, which, when comosed with a torus will give% a triple torus      % center of one of the torii   O = [-Kb, -Kb/sqrt(3)];       angle = 2*pi/3;   Mat = [ cos(angle)  -sin(angle); sin(angle)   cos(angle)  ];      p =[x, y]';   phi = atan2(y, x);      if phi >= pi/6 & phi <= 5*pi/6      p = Mat*p; % rotate 120 degree counterclockwise   elseif phi >= -pi/2 & phi < pi/6      p = Mat*p; p = Mat*p; % rotate 240 degrees counterclockwise   end      x=p(1); y = p(2);   % reflect against a line, to merge two cases in one   if y > x/sqrt(3)            p = [x, y];      v = [cos(2*pi/3), sin(2*pi/3)];            p = p - 2*v*dot(p, v)/dot(v, v);      x = p(1); y = p(2);         end      if x > O(1)      % project to the y axis, to a point B      if y < O(2)                  A = [O(1), y];         B = [0, y];      else         A = O;         p = [x, y];         rho = norm(p-O);         B = O+(Kb/rho)*(p-O);         %         t = -O(1)/(x-O(1));%         B = [0, O(2)+t*(y-O(2))];               end            p = [x, y];            d=norm(p-A);      q = norm(B-A);            d = my_map(d, q, Km);      p = (d/q)*B+(1-d/q)*A;      x=p(1); y=p(2);         end   % shift towards the origin   x = x-O(1);   y = y-O(2);function y = my_map(x, Kb, Km)      if x > Kb      y = Km + 1;   elseif x < Km      y = x;   else      y = Km+sin((pi/2)*(x-Km)/(Kb-Km));   end