//
// Crystal Structures
// - Sodium Chloride (NaCl) Structure
// - Caesium chloride (CsCl) Structure
// - Hexagonal Close-packed (hcp) Structure
// - Diamond Structure
// - Zincblende (ZnS) Structure
// - Perovskite Structure (SrTiO3)
//
// See also my website. There are all renderings.
// http://www.eva.hi-ho.ne.jp/nakajima2001/tadayuki/Scilab/Scilab.html
//
funcprot(0); //nothing special is done when a function is redefined
function [] = PlotAtom(mode, orig, ra)
    t = linspace(0, %pi, 20);
    p = linspace(0, 2*%pi, 40);
    x = ra .* (sin(t)' * cos(p));
    y = ra .* (sin(t)' * sin(p));
    z = ra .* (cos(t)' * ones(p));
    x1 = x + orig(1)*ones(x);
    y1 = y + orig(2)*ones(y);
    z1 = z + orig(3)*ones(z);
    if mode == 1 then
        [xf, yf, zf] = nf3d(x1', y1', z1');
        plot3d1(xf, yf, zf, theta = 300, alpha = 75, leg = "@@", flag = [-1 6 0]);
    elseif mode >= 2 then
        plot3d2(x1, y1, z1, theta = 300, alpha = 75, leg = "@@", flag = [1 6 0]);
        ce=gce(); // store the Compound containing the plot
        ce.foreground = (mode - 2)/8*ColormapSize;
        ce.hiddencolor = mode/8*ColormapSize;
    end
endfunction
function [] = PlotAtoms(mode, A, ra)
    T = [a1 a2 a3]';
    A1 = (A*T)';
    for i = 1:size(A1, "c")
        PlotAtom(mode, A1(:,i), ra);
    end
endfunction
function [A] = GetUnitCell(A0, n1, n2, n3)
    A1 = [];
    for i = n1
        for j = n2
            for k = n3
                A1 = [A1;LatticeTranslate(A0, [i j k])];
            end
        end
    end
    A = [];
    for i = 1:size(A1, "r")
        if and([(A1(i,1) <= max(n1)) (A1(i,2) <= max(n2)) (A1(i,3) <= max(n3))]) then
            A = [A;A1(i,:)];
        end
    end
endfunction
function [A] = RemoveAtoms(A0, RA)
    A = [];
    for i = 1:size(A0, "r")
        plota = %t;
        for j = 1:size(RA, "r")
            if and([(A0(i,1)==RA(j,1)) (A0(i,2)==RA(j,2)) (A0(i,3)==RA(j,3))]) then
                plota =%f;
            end
        end
        if plota then
            A = [A;A0(i,:)];
        end
    end
endfunction
function [A] = LatticeTranslate(A0, v)
    A = ones(A0)*0;
    for i = 1:size(A0, "r")
        A(i,:) = A0(i,:) + v;
    end
endfunction
function [] = ConnectAtoms(mode, A, r0, col, th)
    T = [a1 a2 a3]';
    A1 = (A*T)';
    for i = 2:size(A1, "c")
        for j = 1:i - 1
            xv = [A1(1,i);A1(1,i-j)];
            yv = [A1(2,i);A1(2,i-j)];
            zv = [A1(3,i);A1(3,i-j)];
            select mode,
            case 0 then
                B = [A(i,:) - A(i-j,:)];
            case 1 then
                B = [A1(:,i) - A1(:,i-j)];
            end
            r = sqrt(B(1)^2 + B(2)^2 + B(3)^2);
            if round(100*r)/100 == round(100*r0)/100 then
                xsegs(xv, yv, zv);
                ce=gce(); // store the Compound containing the plot
                ce.segs_color = col;
                ce.thickness = th;
            end
        end
    end
endfunction


//Sodium Chloride (NaCl) Structure
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
Cl = [[0 0 0];[1 1 0]*1/2;[1 0 1]*1/2;[0 1 1]*1/2];
Na = [[1 1 1]*1/2;[0 0 1]*1/2;[0 1 0]*1/2;[1 0 0]*1/2];

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
PlotAtoms(1, GetUnitCell(Cl, 0:1, 0:1, 0:1), 0.07);
PlotAtoms(3, GetUnitCell(Na, 0:1, 0:1, 0:1), 0.05);
NaCl = [GetUnitCell(Cl, 0:1, 0:1, 0:1); GetUnitCell(Na, 0:1, 0:1, 0:1)];
ConnectAtoms(0, NaCl, 1/2, 5/8*ColormapSize, 1);
xtitle('Sodium Chloride (NaCl) Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75, 300 + n];
    show_pixmap();
    n = n + 0.2; 
end
cf.pixmap = 'off';


//Caesium chloride (CsCl) Structure
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
Cs = [0 0 0];
Cl = [1 1 1]*1/2;

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
r = PlotAtoms(3, GetUnitCell(Cs, 0:1, 0:1, 0:1), 0.05);
PlotAtoms(1, GetUnitCell(Cl, 0:1, 0:1, 0:1), 0.07);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:1, 0:1, 0:1), 1, 5/8*ColormapSize, 1);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:1, 0:1, 0:1), sqrt(3), 5/8*ColormapSize, 1);
xtitle('Caesium chloride (CsCl) Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75,300 + n];
    show_pixmap();
    n = n + 0.2; 
end
cf.pixmap = 'off';


//Hexagonal Close-packed (hcp) Structure
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [cos(2/3*%pi);sin(2/3*%pi);0];
a3 = [0;0;sqrt(8/3)];
Ti1 = [0 0 0];
Ti2 = [2/3 1/3 1/2];

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
RA = [[0 2 0];[2 0 0];[0 2 1];[2 0 1]];
HCP1 = RemoveAtoms(GetUnitCell(Ti1, 0:2, 0:2, 0:1), RA);
PlotAtoms(1, HCP1, 0.07);
ConnectAtoms(0, HCP1, 1, 5/8*ColormapSize, 1);
ConnectAtoms(1, HCP1, 1, 5/8*ColormapSize, 1);
HCP2 = RemoveAtoms(GetUnitCell(Ti2, 0:2, 0:2, 0:1), Ti2 + [1 0 0]);
PlotAtoms(1, HCP2, 0.07);
xtitle('Hexagonal Close-packed (hcp) Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [70,290 + n];
    show_pixmap();
    n = n + 0.2; 
end
cf.pixmap = 'off';


//Diamond Structure
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
C1 = [[0 0 0];[0 1 1]*1/2;[1 0 1]*1/2;[1 1 0]*1/2];
C2 = C1 + [[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4];
C = [C1;C2];

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
PlotAtoms(1, GetUnitCell(C, 0:1, 0:1, 0:1), 0.07);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:1, 0:1, 0:1), 1, 5/8*ColormapSize, 1);
ConnectAtoms(0, GetUnitCell(C, 0:1, 0:1, 0:1), sqrt(3)/4, 3/8*ColormapSize, 4);
xtitle('Diamond Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75,300 + n];
    show_pixmap();
    n = n + 0.2; 
end
cf.pixmap = 'off';


//Zincblende (also called Sphalerite) (ZnS) Structure
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
Zn = [[0 0 0];[0 1 1]*1/2;[1 0 1]*1/2;[1 1 0]*1/2];
S = Zn + [[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4];

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
PlotAtoms(1, GetUnitCell(Zn, 0:1, 0:1, 0:1), 0.07);
PlotAtoms(3, GetUnitCell(S, 0:1, 0:1, 0:1), 0.05);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:1, 0:1, 0:1), 1, 5/8*ColormapSize, 1);
ZnS = [GetUnitCell(Zn, 0:1, 0:1, 0:1); GetUnitCell(S, 0:1, 0:1, 0:1)];
ConnectAtoms(0, ZnS, sqrt(3)/4, 3/8*ColormapSize, 4);
xtitle('Zincblende (ZnS) Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75,300 + n];
    show_pixmap();
    n = n + 0.2;
end
cf.pixmap = 'off';


//Zincblende (also called Sphalerite) (ZnS) Structure 2x2x2
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
Zn = [[0 0 0];[0 1 1]*1/2;[1 0 1]*1/2;[1 1 0]*1/2];
S = Zn + [[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4;[1 1 1]*1/4];

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
PlotAtoms(1, GetUnitCell(Zn, 0:2, 0:2, 0:2), 0.07);
PlotAtoms(3, GetUnitCell(S, 0:2, 0:2, 0:2), 0.05);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:2, 0:2, 0:2), 1, 5/8*ColormapSize, 1);
ZnS = [GetUnitCell(Zn, 0:2, 0:2, 0:2); GetUnitCell(S, 0:2, 0:2, 0:2)];
ConnectAtoms(0, ZnS, sqrt(3)/4, 3/8*ColormapSize, 4);
xtitle('Zincblende (ZnS) Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75,300 + n];
    show_pixmap();
    n = n + 0.2;
end
cf.pixmap = 'off';


//Perovskite Structure (SrTiO3)
a1 = [1;0;0]; //a1,a2,a3 :primitive axes
a2 = [0;1;0];
a3 = [0;0;1];
Ti = [0 0 0];
O =  [[0 0 1]*1/2;[0 1 0]*1/2;[1 0 0]*1/2];
Sr = [1 1 1]*1/2;

clf(); //To clear the previous plot
cf = gcf();
cf.pixmap = 'on';
cf.color_map = bonecolormap(32);
ColormapSize = size(cf.color_map, "r"); //get the colormap size
PlotAtoms(1, GetUnitCell(Ti, 0:1, 0:1, 0:1), 0.07);
PlotAtoms(3, GetUnitCell(O, 0:1, 0:1, 0:1), 0.05);
PlotAtoms(5, GetUnitCell(Sr, 0:1, 0:1, 0:1), 0.07);
ConnectAtoms(0, GetUnitCell([0 0 0], 0:1, 0:1, 0:1), 1, 5/8*ColormapSize, 1);
xtitle('Perovskite Structure');
ca = gca(); //get the current axes
ca.tight_limits = "on";
ca.title.font_foreground = 6/8*ColormapSize;
ca.title.font_style = 1; //1 to 10
ca.title.font_size = 4;
n = 0;
while n <= 360
    ca.rotation_angles = [75, 300 + n];
    show_pixmap();
    n = n + 0.2; 
end
cf.pixmap = 'off';