% illustration of a triple torus, obtained as an isosurface
function main()

   % big and small radii of the torus
   R = 3; r = 1; 

   % c controls the transition from one ring to the other
   c = 1.3*pi/2;
   
   Kb = R+r;
  
   h = 0.1; % h is the grid size. Smaller h means prettier picture.
   
   X = (-Kb-h):h:(5*Kb+h); m = length(X);
   Y = (-Kb-h):h:(Kb+h);   n = length(Y);
   Z = (-r-h):h:(r+h);     k = length(Z);
 
   W = zeros(m, n, k); % the zero level set of this function will be the desired shape
 
   for i=1:m
      for j=1:n
         x = X(i); x = my_map(x, Kb, c);   % map from two torii to one torus
         y = Y(j); 
         W(i, j, :) = (sqrt(x^2+y^2)-R)^2 + Z.^2-r^2; % torus eqn, vectorize in Z
      end
   end

   figure(4); clf; hold on; axis equal; axis off;

   H = patch(isosurface(W, 0));
   isonormals(W, H);
      
   light_green=[184, 224, 98]/256;

   % set some propeties
   set(H, 'FaceColor', light_green, '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))
      
% viewing angle
   view(-146, 42);

% add in a source of light
   camlight (-50, 54); lighting phong;

% save as png
  print('-dpng', '-r500', sprintf('Triple_torus_illustration.png'));
   
% This function constructs the second and third rings in the triple torus
% by mapping from the first one.
function y=my_map(x, K, c)

   if x > 2*K
      y = -my_map(4*K - x, K, c);
   else

      if x > K
         x = 2*K - x;
      end
      
      if x < K-c
         y = x;
      else
         y = (K-c) + sin((x - (K-c))*(pi/2/c));
      end
      
   end