function SRplot_potash_regions(SR)
% call with:
% SRplot_potash_regions(data.potash_alum_3.SR_part)
clc
clf
hold on
%% Settings
colors = bone(100);
color1 = colors(25,:);
color2 = colors(75,:);
color3 = colors(50,:);

%% Load data
Av = SR.domain.facetnormals;
% Av = Av([4 3 2 1],:);
cdef.A = SR.facetnormals;
cdef.G = SR.grouping;


%% Complete Cone
P_potash_cone = convertHR2PT([Av; 0 1 0], [0 0 0 0 1]');
ex_potash = extreme(P_potash_cone)';
ex_potash = ex_potash([3 1],1:4)';

%% Partition 1
P_potash_part1 = polytope([SR.se(SR.unified.se_partitions(1,:), :); 0 0 0]);
P_potash_part1 = and(P_potash_cone, 10*P_potash_part1);
ex_potash_part1 = extreme(P_potash_part1)'
ex_potash_part1 = ex_potash_part1([3 1],1:3)'

%% Partition 2
P_potash_part2 = polytope([SR.se(SR.unified.se_partitions(2,:), :); 0 0 0]);
P_potash_part2 = and(P_potash_cone, 10*P_potash_part2);
ex_potash_part2 = extreme(P_potash_part2)'
ex_potash_part2 = ex_potash_part2([3 1],1:3)'

%% Plot Partitions
P_reg = polytope(ex_potash);
P_part1 = polytope(ex_potash_part1);
P_part2 = polytope(ex_potash_part2);
viewPolyTope(P_part1, 'color', color2, 'shade', 1, 'linewidth', 1);
viewPolyTope(P_part2, 'color', color2, 'shade', 1, 'linewidth', 1);
viewPolyTope(P_reg, 'wire', 1, 'linewidth', 2, 'wirestyle', '-');

%% Plot Lines
[A, h] = double(P_part1);
%i = 2;plot([0 2], h(i)/A(i,2) - A(i,1)/A(i,2)*[0 2], 'k-');
i = 3;plot([0 2], h(i)/A(i,2) - A(i,1)/A(i,2)*[0 2], 'k-', 'linewidth', 1);
i = 1;plot([1 1]*h(i)/A(i,1), [0 2], 'k-', 'linewidth', 1);

[A, h] = double(P_part2);
i = 1;plot([0 2], h(i)/A(i,2) - A(i,1)/A(i,2)*[0 2], 'k-', 'linewidth', 1);
%i = 3;plot([0 2], h(i)/A(i,2) - A(i,1)/A(i,2)*[0 2], 'r-');
i = 2;plot([1 1]*h(i)/A(i,1), [0 2], 'k-', 'linewidth', 1);

%% Regions
X = [0.6 0.9;0.6 0.35;1.5 1.05;1.45 1.9;...
     1.0 1.5;0.57 0.6;1 0.55;1.5 1.5;...
     1.2 1.2;0.75 0.82];
labels = {'I', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VIII', 'IX', 'X'};
for i=1:size(X,1)
    text(X(i,1), X(i,2), labels{i}, 'FontSize', 8);
end
labels = {'', '', '', '', '', '', '', '', '', ''};

lim = 2;
xlim([0 lim]);ylim([0 lim]);
set(gca, 'FontSize', 8);
set(1,'renderer', 'painters', ...
      'PaperUnits', 'centimeters', ...
      'PaperSize', [9, 9], ...
      'PaperPosition', [0 0 9 9]);
xlabel('Relative growth rate R_{110}/R_{100}');
ylabel('Relative growth rate R_{111}/R_{100}');
box on
print -depsc potash_regions.eps

%% Print Structuring Elements
% close all
%se_dash = SR.se([1 2 4 3],:);
se_dash = SR.se([2 4 3 1],:);
scale = [0.5 0.5 0.5 0.5];
for i = 1:size(se_dash, 1)
    hC = se_dash(i, :)';
    internal_print(['potash_se_' num2str(i)], cdef, hC, color3, labels{i}, scale(i));
end

%% Print Domain Examples
scale = [0.5 0.5];
for i = 1:size(SR.unified.se_partitions, 1)
    outfile = ['potash_part_' num2str(i)];
    hC = SR.se(SR.unified.se_partitions(i,:), :)';
    hC = hC * ones([size(hC, 2) 1]);
    internal_print(outfile, cdef, hC, color3, labels{i+8}, scale(i));
end
% 
%% Print combinations of two
ij = [1 4; 1 2; 2 3; 3 4];
scale = [0.5 0.5 0.5 0.5];
for i = 1:size(SR.se, 1)
    outfile = ['potash_sesum_' num2str(ij(i,1)) num2str(ij(i,2))];
    hC = se_dash([ij(i,1),ij(i,2)], :)';
    hC = hC * ones([size(hC, 2) 1]);
    internal_print(outfile, cdef, hC, color3, labels{i+4}, scale(i));
end

%% 3D-cone
close all
ex1 = SR.se(SR.simplex.se_partitions(1,:),:)
ex2 = SR.se(SR.simplex.se_partitions(2,:),:)
ex3 = [1 0 0;0 1 0;0 0 1];
ex3 = [ex3-2e-3;ex3-1e-3];
n0 = [1 1 1]/3*sum(ex1(1,:));
scale = 3;
% morph cone
for i = 1:3
    ex1(i,:) = (1-scale)*n0 + scale*ex1(i,:);
    ex2(i,:) = (1-scale)*n0 + scale*ex2(i,:);
    ex1(i,:) = ex1(i,:)/sum(ex2(end,:));
    ex2(i,:) = ex2(i,:)/sum(ex2(end,:));
end

% View Polytopes
P1 = polytope([ex1; [0 0 0]]);
P2 = polytope([ex2; [0 0 0]]);
P3 = polytope(ex3);
viewPolyTope(P1, 'color', color1)
hold on
viewPolyTope(P2, 'color', color3)
viewPolyTope(P3, 'color', color2, 'shade', 0.2);

% View SE's
shift = [0.11 -0.04 -0.11;0 -0.1 0.07;0 0.02 0.02;0.04 0 -0.04];
for i = 1:3
    p = [ex1(i,1), ex1(i,2), ex1(i,3)]*(1+2e-2);
    plot3(p(1), p(2), p(3), 'k.', 'MarkerSize', 20);
    text(p(1)+shift(i,1), p(2)+shift(i,2), p(3)+shift(i,3), ['S' num2str(i)]);
end
plot3(ex2(2,1), ex2(2,2), ex2(2,3), 'k.', 'MarkerSize', 20);
text(ex2(2,1)+shift(4,1), ex2(2,2)+shift(4,2), ex2(2,3)+shift(4,3), 'S4');


view([50 15]);box on
set(gca, 'xtick', [0 1], 'ytick', [0 1], 'ztick', [0 1], ...
         'xlim', [0 1], 'ylim', [0 1], 'zlim', [0 1]);
xlabel('h1');ylabel('h2');zlabel('h3');
set(1,'renderer', 'painters');
set(gca, 'FontSize', 20, 'linewidth', 1);
saveas(1, 'pseudo_potash_cone.eps');

end

function internal_print(outfile, cdef, hC, color, label, scale)
% takes filename and polytope
% (box fitting) > enlarge axes > print jpg
    hf = figure;
    cmap = bone(4+size(cdef.G, 2));
    cmap = cmap(2+(1:size(cdef.G, 2)), :);
    P = convertHR2PT(cdef.A, cdef.G*hC);
    viewCrystal(cdef, ...
            cdef.G*hC,...
            'colormap', cmap, ...
            'colorbar', 0, ...
            'annotation', 'none', ...
            'flatness', 1e-3);
%     ratio_a = [1 1];
%     ratio_f = [1 1];
    if exist('scale', 'var')
        lim = get(gca, {'xlim', 'ylim', 'zlim'});
        viewCrystal(cdef, ...
            cdef.G*hC*scale,...
            'colormap', cmap, ...
            'colorbar', 0, ...
            'annotation', 'none', ...
            'flatness', 1e-3);
%         for i = 1:3
%             tmp = (lim{i}(2) - lim{i}(1));
%             tmp = tmp*scale - tmp;
%             lim{i}(1) = lim{i}(1) + tmp/2;
%             lim{i}(2) = lim{i}(2) - tmp/2;
%         end
        set(gca, 'xlim', lim{1}, 'ylim', lim{2}, 'zlim', lim{3});
    end
    if exist('label', 'var')
        [lower upper] = getBoxPT(P);
        X = ([5 5 9]'.*lower + [5 5 1]'.*upper)/10;
        text(X(1), X(2), X(3), ['      ' label], 'FontSize', 12);
    end
%     set(hf, 'Position', [50 50,ratio_f*300]);
%     set(hf, 'PaperSize', ratio_f*8,...
%         'PaperPosition', [0 0, ratio_f*8], 'Color', 'none');
%     set(hf, 'Renderer', 'opengl')
%     set(gca, 'PlotBoxAspectRatio', [ratio_a, 1], 'Visible', 'off', 'Box', 'off',...
%         'Position', [-0 -0 1 1],...
%         'Clipping', 'on');
    view([-35 20]);
    box off,axis off
    set(gcf, 'PaperUnits', 'centimeters', ...
        'PaperSize', [4, 4], ...
        'PaperPosition', [0 0 5 4],...
        'renderer', 'painter');
    set(gca, 'Position', [0 0 1 1]);
    eval(['print -depsc ' outfile '.eps']);
    close(hf);
end
