function RGBcube(val,num)
%RGBcube.m displays sections of the colorcube

if nargin~=2    %Initialize
	ga=0.3*sqrt(3); %ga must be in the open interval (0,1) 
    al=0;
	th=0;
    
	fig=figure('color','k',...
      'Position',[232 150 560 560],...%'resize','off',...
      'Menubar','none',...
      'Doublebuffer','on',...
      'userdata',[ga,al,th]);

	%Button for gamma
	gammaslide=uicontrol('style','slider',...  
              'pos',[20 20 40 20],...
              'min',0,...
              'max',sqrt(3),...
              'fontsize',12,...
              'userdata','gamma',...
              'value',ga,...
              'callback','RGBcube(get(findobj(''Parent'',gcf,''userdata'',''gamma''),''value''),1)');
	%Button for alpha
	alphaslide=uicontrol('style','slider',...
              'pos',[20 40 40 20],...
              'min',-0.5*pi,...
              'max',0.5*pi,...
              'fontsize',12,...
              'userdata','alpha',...
              'value',al,...
              'callback','RGBcube(get(findobj(''Parent'',gcf,''userdata'',''alpha''),''value''),2)');
	%Button for theta
	thetaslide=uicontrol('style','slider',...
              'pos',[20 60 40 20],...
              'min',0,...
              'max',2*pi,...
              'fontsize',12,...
              'userdata','theta',...
              'value',th,...
              'callback','RGBcube(get(findobj(''Parent'',gcf,''userdata'',''theta''),''value''),3)');
     
	%Perform computations
	[vert,color,cb,ct,c1,c2,c3,c4,c5,c6]=compute(ga,al,th);

	%Plot
	slice=patch(vert(1,:),vert(2,:),vert(3,:),color);
    set(slice,'userdata','patch');
	line('color','w', ...
     'Linestyle','-',...
     'userdata',1,...
     'xdata',[cb(1) c1(1)],...
     'ydata',[cb(2) c1(2)],...
     'zdata',[cb(3) c1(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',2,...
     'xdata',[cb(1) c2(1)],...
     'ydata',[cb(2) c2(2)],...
     'zdata',[cb(3) c2(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',3,...
     'xdata',[cb(1) c3(1)],...
     'ydata',[cb(2) c3(2)],...
     'zdata',[cb(3) c3(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',4,...
     'xdata',[ct(1) c4(1)],...
     'ydata',[ct(2) c4(2)],...
     'zdata',[ct(3) c4(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',5,...
     'xdata',[ct(1) c5(1)],...
     'ydata',[ct(2) c5(2)],...
     'zdata',[ct(3) c5(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',6,...
     'xdata',[ct(1) c6(1)],...
     'ydata',[ct(2) c6(2)],...
     'zdata',[ct(3) c6(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',7,...
     'xdata',[c1(1) c4(1)],...
     'ydata',[c1(2) c4(2)],...
     'zdata',[c1(3) c4(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',8,...
     'xdata',[c4(1) c2(1)],...
     'ydata',[c4(2) c2(2)],...
     'zdata',[c4(3) c2(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',9,...
     'xdata',[c2(1) c5(1)],...
     'ydata',[c2(2) c5(2)],...
     'zdata',[c2(3) c5(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',10,...
     'xdata',[c5(1) c3(1)],...
     'ydata',[c5(2) c3(2)],...
     'zdata',[c5(3) c3(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',11,...
     'xdata',[c3(1) c6(1)],...
     'ydata',[c3(2) c6(2)],...
     'zdata',[c3(3) c6(3)]);
	line('color','w', ...
     'Linestyle','-',...
     'userdata',12,...
     'xdata',[c6(1) c1(1)],...
     'ydata',[c6(2) c1(2)],...
     'zdata',[c6(3) c1(3)]);	
	view(3)
	axis equal
	axis off
else
    a=get(gcf,'userdata');
    a(num)=val;
    set(gcf,'userdata',a);
    [vert,color,cb,ct,c1,c2,c3,c4,c5,c6]=compute(a(1),a(2),a(3));
    %Plot
    set(findobj('Parent',gca,'userdata','patch'),...
     'xdata',vert(1,:),...
     'ydata',vert(2,:),...
     'zdata',vert(3,:),...
     'cdata',color);
	set(findobj('Parent',gca,'userdata',1),...
     'xdata',[cb(1) c1(1)],...
     'ydata',[cb(2) c1(2)],...
     'zdata',[cb(3) c1(3)]);
	set(findobj('Parent',gca,'userdata',2),...
     'xdata',[cb(1) c2(1)],...
     'ydata',[cb(2) c2(2)],...
     'zdata',[cb(3) c2(3)]);
	set(findobj('Parent',gca,'userdata',3),...
     'xdata',[cb(1) c3(1)],...
     'ydata',[cb(2) c3(2)],...
     'zdata',[cb(3) c3(3)]);
	set(findobj('Parent',gca,'userdata',4),...
     'xdata',[ct(1) c4(1)],...
     'ydata',[ct(2) c4(2)],...
     'zdata',[ct(3) c4(3)]);
	set(findobj('Parent',gca,'userdata',5),...
     'xdata',[ct(1) c5(1)],...
     'ydata',[ct(2) c5(2)],...
     'zdata',[ct(3) c5(3)]);
	set(findobj('Parent',gca,'userdata',6),...
     'xdata',[ct(1) c6(1)],...
     'ydata',[ct(2) c6(2)],...
     'zdata',[ct(3) c6(3)]);
	set(findobj('Parent',gca,'userdata',7),...
     'xdata',[c1(1) c4(1)],...
     'ydata',[c1(2) c4(2)],...
     'zdata',[c1(3) c4(3)]);
	set(findobj('Parent',gca,'userdata',8),...
     'xdata',[c4(1) c2(1)],...
     'ydata',[c4(2) c2(2)],...
     'zdata',[c4(3) c2(3)]);
	set(findobj('Parent',gca,'userdata',9),...
     'xdata',[c2(1) c5(1)],...
     'ydata',[c2(2) c5(2)],...
     'zdata',[c2(3) c5(3)]);
	set(findobj('Parent',gca,'userdata',10),...
     'xdata',[c5(1) c3(1)],...
     'ydata',[c5(2) c3(2)],...
     'zdata',[c5(3) c3(3)]);
	set(findobj('Parent',gca,'userdata',11),...
     'xdata',[c3(1) c6(1)],...
     'ydata',[c3(2) c6(2)],...
     'zdata',[c3(3) c6(3)]);
	set(findobj('Parent',gca,'userdata',12),...
     'xdata',[c6(1) c1(1)],...
     'ydata',[c6(2) c1(2)],...
     'zdata',[c6(3) c1(3)]);
end
      
function [vert,color,cb,ct,c1,c2,c3,c4,c5,c6]=compute(ga,al,th);
%Reference configuration
e1=[sqrt(2/3);0;1/sqrt(3)];
e2=[-1/sqrt(6);1/sqrt(2);1/sqrt(3)];
e3=[-1/sqrt(6);-1/sqrt(2);1/sqrt(3)];

%Rotate according to al and th
R1=[cos(th) sin(th) 0;-sin(th) cos(th) 0; 0 0 1];
R2=[cos(al) 0 -sin(al);0 1 0;sin(al) 0 cos(al)];
e1=R2*R1*e1;
e2=R2*R1*e2;
e3=R2*R1*e3;

%Get upward pointing vector in terms of ei
E3=[e1(3);e2(3);e3(3)];

%Construct cube corners
u=(e1+e2+e3)/sqrt(3);   %unit vector along diagonal
cb=-ga*u;
c1=e1-ga*u;
c2=e2-ga*u;
c3=e3-ga*u;
c4=e1+e2-ga*u;
c5=e2+e3-ga*u;
c6=e1+e3-ga*u;
ct=(sqrt(3)-ga)*u;

%Find polygon vertices
vert=[];
%first from the bottom
if sign(cb(3))~=sign(c1(3))
    vert=[vert cb+(c1-cb)*(abs(cb(3))/(abs(cb(3))+abs(c1(3))))];
end
if sign(cb(3))~=sign(c2(3))
    vert=[vert cb+(c2-cb)*(abs(cb(3))/(abs(cb(3))+abs(c2(3))))];
end
if sign(cb(3))~=sign(c3(3))
    vert=[vert cb+(c3-cb)*(abs(cb(3))/(abs(cb(3))+abs(c3(3))))];
end
%then from the top
if sign(ct(3))~=sign(c4(3))
    vert=[vert ct+(c4-ct)*(abs(ct(3))/(abs(ct(3))+abs(c4(3))))];
end
if sign(ct(3))~=sign(c5(3))
    vert=[vert ct+(c5-ct)*(abs(ct(3))/(abs(ct(3))+abs(c5(3))))];
end
if sign(ct(3))~=sign(c6(3))
    vert=[vert ct+(c6-ct)*(abs(ct(3))/(abs(ct(3))+abs(c6(3))))];
end
%and lastley around the middle
if sign(c1(3))~=sign(c4(3)) & c1(3)*c4(3)~=0
    vert=[vert c1+(c4-c1)*(abs(c1(3))/(abs(c1(3))+abs(c4(3))))];
end
if sign(c4(3))~=sign(c2(3)) & c4(3)*c2(3)~=0
    vert=[vert c4+(c2-c4)*(abs(c4(3))/(abs(c4(3))+abs(c2(3))))];
end
if sign(c2(3))~=sign(c5(3)) & c2(3)*c5(3)~=0
    vert=[vert c2+(c5-c2)*(abs(c2(3))/(abs(c2(3))+abs(c5(3))))];
end
if sign(c5(3))~=sign(c3(3)) & c5(3)*c3(3)~=0
    vert=[vert c5+(c3-c5)*(abs(c5(3))/(abs(c5(3))+abs(c3(3))))];
end
if sign(c3(3))~=sign(c6(3)) & c3(3)*c6(3)~=0
    vert=[vert c3+(c6-c3)*(abs(c3(3))/(abs(c3(3))+abs(c6(3))))];
end
if sign(c6(3))~=sign(c1(3)) & c6(3)*c1(3)~=0
    vert=[vert c6+(c1-c6)*(abs(c6(3))/(abs(c6(3))+abs(c1(3))))];
end

%Order the vertices
for n=4:size(vert,2)
    ind=1;
    ang=0;
    %between current and next up to where next is last.
    for k=1:n-2
        dottemp=(vert(:,k)-vert(:,n))'*(vert(:,k+1)-vert(:,n));
        if dottemp  %only divide if nonzero
            dottemp=dottemp/norm(vert(:,k)-vert(:,n));
            dottemp=dottemp/norm(vert(:,k+1)-vert(:,n));
        end
        newang=acos(dottemp);
        if newang>ang
            ang=newang;
            ind=k;
        end
    end
    %between last and first.
    dottemp=(vert(:,n-1)-vert(:,n))'*(vert(:,1)-vert(:,n));
    if dottemp  %only divide if nonzero
        dottemp=dottemp/norm(vert(:,n-1)-vert(:,n));
        dottemp=dottemp/norm(vert(:,1)-vert(:,n));
    end
    newang=acos(dottemp);
    if newang>ang
        ind=n-1;
    end
    vert=[vert(:,1:ind) vert(:,n) vert(:,ind+1:n-1) vert(:,n+1:end)];
end    

%Compute Colors
for n=1:size(vert,2)
    color(n,1,1:3)=[e1'*(vert(:,n)+ga*u);
                    e2'*(vert(:,n)+ga*u);
                    e3'*(vert(:,n)+ga*u)];
end