function Bouncing_Bucky(f)
    % Bouncing_Bucky, no arguments, random self-initializing infinite loop.
    % Bouncing_Bucky(1), initialize and take one step.
    % Bouncing_Bucky(f), f > 1, take one step.
    % https://blogs.mathworks.com/cleve

    % Frames per second
    fps = 24;

    if nargin == 0 || f == 1
        init_graphics;

        % Adjacency matrix and vertices of truncated icosahedron
        [B,V] = bucky;
        V = V/8;

        % Color pentagons and hexagons
        blue_gray = 0.9*[0 0 .5; 1 1 1];
        C = allcycles(graph(B), maxCycleLength = 6);
        L = cellfun(@length, C);
        F = NaN(height(C), max(L));
        for k = 1:height(C)
            F(k,1:L(k)) = C{k};
        end
        FaceColor = blue_gray((L == max(L))+1, :);

        % Handle Graphics Transform object.
        H = hgtransform(Parent = gca);

        % Patches
        patch(Parent = H, ...
            Faces = F, ...
            Vertices = V, ...
            EdgeColor = blue_gray(1,:), ...
            FaceVertexCData = FaceColor, ...
            FaceColor = 'flat')

        % Initial position {x,y,dx,dy} 
        if nargin > 0
            initial = {0.5, 0, 1/fps, 1/fps};
        else
            initial = {rand, rand, randn/fps, randn/fps};
        end
        set(gca,'userdata',initial)
    end
    
    while bouncing

        % One step
        g = get(gca,'userdata');
        [x,y,dx,dy] = deal(g{:});
        x = x + dx;
        y = y + dy;

        % Bounce
        if x < 0 || x > 1
            dx = -dx;
        end
        if y < 0 || y > 1
            dy = -dy;
        end
    
        % Transform matrix with fake rotation
        beta = 90*(2*x+y);
        M = Ty(y) * Tx(x) * Ry(beta) * Rx(2*beta);
        H = get(gca,'children');
        H.Matrix = M;
    
        % Save position for next frame
        set(gca,'userdata',{x,y,dx,dy})

        if nargin > 0
            break
        else
            pause(1/fps)
        end
    end

    % -------------------------------------------------------

    function R = Rx(t)
        R = [1        0        0  0
             0  cosd(t)  sind(t)  0
             0 -sind(t)  cosd(t)  0
             0        0        0  1];
    end

    function R = Ry(t)
        R = [cosd(t)  0  sind(t)  0
             0        1        0  0
            -sind(t)  0  cosd(t)  0
             0        0        0  1];
    end

    function T = Tx(x)
        T = eye(4,4); 
        T(1,4) = x;
    end

    function T = Ty(y)
        T = eye(4,4); 
        T(2,4) = y;
    end

    function init_graphics
        clf
        axis([0 1 0 1])
        axis square 
        box on
        set(gca,'xtick',[],'ytick',[])
        uicontrol('style','toggle', ...
            'string','<=>', ...
            'fontsize',12, ...
            'fontweight','bold', ...
            'units','normalized', ...
            'position',[.90 .88 .07 .07], ...
            'callback',@restart);
    end
     
    function b = bouncing
       b = ~get(findobj('string','<=>'),'value');
    end
     
    function restart(arg,~)
       set(arg,'value',0)
       Bouncing_Bucky
    end
 end