% Program to generate cartoons of different crystal sample types 
% Single crystal, twinned, mosaic spread, textured, powder, and
% nanocrystalline. 

% Uses an adaptation of the function platonic_solid, see 
% https://ch.mathworks.com/matlabcentral/fileexchange/28213-platonic_solid 
% Courtesy Kevin Moerman 

clear; close all;

disp('Choose one of the following crystal types: Single (S), Twinned (T), Mosaic spread (M),');
prompt = 'teXtured (X), Powder (P), nanocrystalline (N) ';
str = input(prompt,'s');

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
str1 = 'CrystalType_';
str3 = '.mp4';
fileName = [str1,str,str3]; % Name of file depends on type of crystal sample 
vid = VideoWriter(fileName,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
set(gca, 'Projection','perspective');
lp = [-1 -1 0.1];
myBlue = [0.5 0.55 0.88];

% Single crystal 
if (str == 'S') || (str == 's') 
    for phi = 0:2:358
        newplot
        hold off 
        LL = 1.6;
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis square
        axis off
        
        [V,F] = platonic_solid(3,1); % Octahedron, size = 1
        ps = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction = [0 0 1];
        rotate(ps,direction,phi,[0 0 0]) 
        material shiny 
        hold on
        plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant

        set(gca, 'Projection','perspective');
        set(gca,'View',[10,20]);
        light('Position',lp,'Style','infinite');
        frame = getframe(gcf);
        writeVideo(vid,frame);
        hold off 
    end
end

% Twinned crystal 
if (str == 'T') || (str == 't') 
    for phi = 0:2:358
        newplot
        hold off 
        LL = 1.6;
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis square
        axis off
        
        [V,F] = platonic_solid(3,1); % Octahedron, size = 1 
        V(:,3) = V(:,3)+0.586;
        ps1 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction1 = [0 1 0];
        rotate(ps1,direction1,35,[0 0 -0.414])
        direction2 = [0 0 1];
        rotate(ps1,direction2,phi,[0 0 0]) 
        [V,F] = platonic_solid(3,1); % Octahedron, size = 1 
        V(:,3) = V(:,3)+0.586;
        ps2 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction1 = [0 1 0];
        rotate(ps2,direction1,-35,[0 0 -0.414])
        direction2 = [0 0 1];
        rotate(ps2,direction2,phi,[0 0 0])
        [V,F] = platonic_solid(3,1); % Octahedron, size = 1 
        V(:,3) = V(:,3)-0.586;
        ps3 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction1 = [0 1 0];
        rotate(ps3,direction1,35,[0 0 0.414])
        direction2 = [0 0 1];
        rotate(ps3,direction2,phi,[0 0 0])
        [V,F] = platonic_solid(3,1); % Octahedron, size = 1 
        V(:,3) = V(:,3)-0.586;
        ps4 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
            'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
        direction1 = [0 1 0];
        rotate(ps4,direction1,-35,[0 0 0.414])
        direction2 = [0 0 1];
        rotate(ps4,direction2,phi,[0 0 0])
        
        material shiny
        hold on
        plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant

        set(gca, 'Projection','perspective');
        set(gca,'View',[10,20]);
        light('Position',lp,'Style','infinite');
        frame = getframe(gcf);
        writeVideo(vid,frame);
        hold off 
    end
end

% Crystal with Mosaic spread 
if (str == 'M') || (str == 'm')
    randOrient = rand(5,5,5,3); % 5x5x5 mosaic blocks, each with its own three tilt angles
    for phi = 0:2:358
        newplot
        hold off
        LL = 4.5;
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis square
        axis off
                
        for xx = -2:2
            for yy = -2:2
                for zz = -2:2
                    [V,F] = platonic_solid(2,0.88); % Cube, size = 0.88
                    V(:,1) = V(:,1)+xx;
                    V(:,2) = V(:,2)+yy;
                    V(:,3) = V(:,3)+zz;
                    ps1 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
                        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
                    direction1 = [1 0 0];
                    rotate(ps1,direction1,7*(randOrient(xx+3,yy+3,zz+3,1)-0.5),[xx yy zz])
                    direction2 = [0 1 0];
                    rotate(ps1,direction2,7*(randOrient(xx+3,yy+3,zz+3,2)-0.5),[xx yy zz])
                    direction3 = [0 0 1];
                    rotate(ps1,direction3,7*(randOrient(xx+3,yy+3,zz+3,3)-0.5),[xx yy zz])
                    
                    direction4 = [0 0 1];
                    rotate(ps1,direction4,phi,[0 0 0])
                                        
                    material shiny
                    hold on
                end
            end
        end
        plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant
        
        set(gca, 'Projection','perspective');
        set(gca,'View',[10,20]);
        light('Position',lp,'Style','infinite');
        frame = getframe(gcf);
        writeVideo(vid,frame);
        hold off
        
    end % end for rotation loop
end

% teXtured sample 
if (str == 'X') || (str == 'x') 
    randOrient = rand(7,7,11); % 7x7x11 texture blocks, randomly oriented in x-y plane 
    for phi = 0:2:358
        newplot
        hold off
        LL = 5;
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis square
        axis off
                
        for xx = -3:3
            for yy = -3:3
                for zz = -5:5
                    [V,F] = platonic_solid(2,0.8); % Cube, size = 0.8
                    V(:,1) = V(:,1)+xx;
                    V(:,2) = V(:,2)+yy;
                    V(:,3) = (7/11)*(V(:,3)+zz);
                    ps1 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
                        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
                    direction3 = [0 0 1];
                    rotate(ps1,direction3,90*(randOrient(xx+4,yy+4,zz+6)-0.5),[xx yy 7*zz/11])
                    
                    direction4 = [0 0 1];
                    rotate(ps1,direction4,phi,[0 0 0])
                                        
                    material shiny
                    hold on
                end
            end
        end
        plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant
        
        set(gca, 'Projection','perspective');
        set(gca,'View',[10,20]);
        light('Position',lp,'Style','infinite');
        frame = getframe(gcf);
        writeVideo(vid,frame);
        hold off
    end % end for rotation loop
end

% Powder or nanocrystalline sample
if (str == 'P') || (str == 'p') || (str == 'N') || (str == 'x') 
    if (str == 'P') || (str == 'p')
        LL = 6;
    end
    if (str == 'N') || (str == 'n')
        LL = 10;
    end

    randOrient = rand(2*LL-1,2*LL-1,2*LL-1,7); % (2LL+1)^3 mosaic blocks, each with its own 
    % three tilt angles, crystal size, and coordinates x, y, and z
    for phi = 0:2:358
        newplot
        hold off
        xlim([-LL LL]);
        ylim([-LL LL]);
        zlim([-LL LL]);
        
        axis equal
        axis square
        axis off
                
        for xx = -LL+1:LL-1
            for yy = -LL+1:LL-1
                for zz = -LL+1:LL-1 
                    xxx = randOrient(xx+LL,yy+LL,zz+LL,5) - 0.5; 
                    yyy = randOrient(xx+LL,yy+LL,zz+LL,6) - 0.5; 
                    zzz = randOrient(xx+LL,yy+LL,zz+LL,7) - 0.5; 
                    radius = (xxx^2 + yyy^2 + zzz^2)^0.5;
                    if (radius <= 0.525)
                    xstalSize = 0.14+0.37*randOrient(xx+LL,yy+LL,zz+LL,4);
                    [V,F] = platonic_solid(3,xstalSize); % Octahedron, random size  
                    V(:,1) = V(:,1)+2*(LL-1)*xxx;
                    V(:,2) = V(:,2)+2*(LL-1)*yyy;
                    V(:,3) = V(:,3)+2*(LL-1)*zzz;
                    ps1 = patch('Faces',F,'Vertices',V,'FaceColor',myBlue,'FaceAlpha',1, ...
                        'EdgeColor','none','FaceLighting','flat','DiffuseStrength',1,'AmbientStrength',1,'SpecularStrength',0);
                    direction1 = [1 0 0];
                    rotate(ps1,direction1,90*(randOrient(xx+LL,yy+LL,zz+LL,1)-0.5),[xxx yyy zzz])
                    direction2 = [0 1 0];
                    rotate(ps1,direction2,90*(randOrient(xx+LL,yy+LL,zz+LL,2)-0.5),[xxx yyy zzz])
                    direction3 = [0 0 1];
                    rotate(ps1,direction3,90*(randOrient(xx+LL,yy+LL,zz+LL,3)-0.5),[xxx yyy zzz])
                    
                    direction4 = [0.1 0.1 1];
                    rotate(ps1,direction4,phi,[0 0 0])
                                        
                    material shiny
                    hold on 
                    end 
                end
            end
        end
        plot3([-LL LL],[-LL LL],[-LL -LL],'color','white', 'LineWidth', 0.1); % dummy line to keep FOV constant
        phi % Monitor on screen progression, as this can take some time for LL = 10 
        set(gca, 'Projection','perspective');
        set(gca,'View',[10,20]);
        light('Position',lp,'Style','infinite');
        frame = getframe(gcf);
        writeVideo(vid,frame);
        hold off
        
    end % end for rotation loop
end

% Output the movie as an mpg file
close(vid);

