% illustration of a double 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:(3*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(-165, 42);  % add in a source of light    camlight (-50, 54); lighting phong;  % save as png   print('-dpng', '-r500', sprintf('Double_torus_illustration.png'));     % This function constructs the second ring in the double torus % by mapping from the first one. function y=my_map(x, K, c)     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