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