function Soma(idx)
% Soma cube, June 10, 2022. 
% See 
% https://blogs.mathworks.com/cleve/2016/03/28/piet-hein-super-ellipses-and-soma-cubes
% https://blogs.mathworks.com/cleve/2022/04/12/digital-simulation-of-rubiks-cube-with-qube/
% \Program Files\MATLAB\Rxxxx\toolbox\matlab\demos\soma.m
% \Program Files\MATLAB\Rxxxx\toolbox\matlab\demos\somasols.m

% Copyright 2022 Cleve Moler

    if nargin == 0
        somafigure
        soln = findobj(gcf,'tag','soln');
        idx = str2double(soln.String);
    else
        soln = findobj(gcf,'tag','soln');
        soln.String = int2str(idx);
    end
    S = somasols;
    S = S(idx,:);
    somashow(S); 
    for k = 1:7
        somashow(S,k);
    end
end

function Q = Q0
    % From Qube, the Digital Rubik's Cube.
    v = 0.90*qzero;
    Q = cell(3,3,3);
    for z = [-2 0 2]
        for y = [-2 0 2]
            for x = [-2 0 2]
                Q{x/2+2,y/2+2,z/2+2} = v + [x y z];
            end
        end 
    end
end

function q0 = qzero
    % Unit cubelet.
    q0 = [-1 -1 -1
          -1 -1  1
          -1  1 -1
          -1  1  1
           1 -1 -1
           1 -1  1
           1  1 -1 
           1  1  1];
end

function somashow(S,k)
% somashow(S)  Display solution S. 

    Q = Q0;
    F = [1 5 7 3
         3 7 8 4
         1 3 4 2
         2 4 8 6
         1 2 6 5
         5 6 8 7];
    S = [0 S];
    C = get(0,'defaultaxescolororder');
    if nargin == 1
        for j = 1:27
            k = ceil(find(S==j)/4);
            patch( ...
                Vertices = Q{j}, ...
                Faces = F, ...
                FaceColor = C(k,:), ...
                LineWidth = 1);
        end
    else
        switch k
            case 1, jk = [13 14 16];    % V
            case 2, jk = [13 14 15 16]; % L
            case 3, jk = [13 14 15 17]; % T
            case 4, jk = [13 14 17 18]; % Z
            case 5, jk = [14 15 17 24]; % S
            case 6, jk = [14 17 18 27]; % R
            case 7, jk = [14 15 17 23]; % Y
        end
        t = 2*pi*k/7;
        shift = 15*[sin(t) cos(t) 0];
        for j = 1:length(jk)
            patch( ...
                Vertices = .75*Q{jk(j)} + shift, ...
                Faces = F, ...
                FaceColor = C(k,:), ...
                LineWidth = 1);
        end
    end
end

function somafigure
% Initialize figure window for Soma cube
        clf
        idx = randi(240);
        set(gcf,'name','Soma', ...
            'numbertitle','off', ...
            'color',0.75*[1 1 1])
        uicontrol('Style','pushbutton', ...
            'String','+', ...
            'Units','normalized', ...
            'Position',[.56 .93 .08 .04], ...
            'Fontsize',12, ...
            'Fontweight','bold', ...
            'Callback',@somacb)
        uicontrol('Style','pushbutton', ...
            'String','-', ...
            'Units','normalized', ...
            'Position',[.36 .93 .08 .04], ...
            'Fontweight','bold', ...
            'Fontsize',12, ...
            'Callback',@somacb)
        uicontrol( ...
            'units','normalized', ...
            'position',[.46 .93 .08 .04], ...
            'style','pushbutton', ...
            'string',int2str(idx), ...
            'Fontweight','bold', ...
            'Fontsize',9, ...
            'tag','soln', ...
            'callback',@somacb);
        set(gca,'position',[0 0 1 1], ...
            'clipping','off')
        axis(14*[-1 1 -1 1 -1 1])
        axis('square','off','vis3d')
        rotate3d on
        view(3)
        shg
end

function somacb(arg,~)
    soln = findobj(gcf,'tag','soln');
    idx = str2double(soln.String);
    switch arg.String
        case '+', idx = mod(idx,240)+1;
        case '-', idx = mod(idx-2,240)+1;
        otherwise, idx = randi(240);
    end
    soln.String = num2str(idx);
    Soma(idx)
end