% Movie schematic from a single unit cell (of either NaCl, diamond, GaAs, or
% benzyl penicillin) to a crystal

clear all; close all;

prompt = 'Which unit cell?: diamond (d/D), GaAs (g/G), NaCl (s/S), or penicillin (p/P)? ';
str = input(prompt,'s');

if (str == 'P') || (str == 'p')
    ucData = importdata('penicillinCoordinates.dat'); % (23 x 5) matrix
    C = [1 1 0; 1 0 0; 1 0 0; 1 0 0; 1 0 0; 0 0 1; 0 0 1; 0.4 0.4 0.4; ...
        0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ...
        0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; ...
        0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4; 0.4 0.4 0.4];
    noAtoms = 23;
    ucxStep = 1;
    ucyStep = 0.5;
    uczStep = 0.5;
    ucMax = max(ucData(:,1));
    ucMin = min(ucData(:,1));
    ucSize = max(ucData(:,5)) + (ucMax - ucMin);
    youChoose = 0.74; % How much smaller than u.c. the molecule is
    sphereRes = 16;
    azStep = 0.4; % Causes camera to swing azimuthally around the crystal as it
    % is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'S') || (str == 's')
    ucData = importdata('NaClCoordinates.dat'); % (8 x 5) matrix
    % Fuschia for Na, green for Cl
    C = [0.631 0.376 0.918; 0.631 0.376 0.918; 0.631 0.376 0.918; ...
        0.631 0.376 0.918; 0.443 0.925 0.302; 0.443 0.925 0.302; ...
        0.443 0.925 0.302; 0.443 0.925 0.302];
    noAtoms = 8;
    ucxStep = 1;
    ucyStep = 1;
    uczStep = 1;
    ucSize = 5.64; % u.c. lattice constant of NaCl
    youChoose = 1.00; % How much smaller than u.c. the molecule is
    sphereRes = 35;
    azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
    % is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'D') || (str == 'd')
    ucData = importdata('diamondCoordinates.dat'); % (8 x 5) matrix
    C = [0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; ...
        0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25; 0.25 0.25 0.25];
    noAtoms = 8;
    ucxStep = 1;
    ucyStep = 1;
    uczStep = 1;
    ucSize = 3.567; % u.c. lattice constant of diamond
    youChoose = 1.00; % How much smaller than u.c. the molecule is
    sphereRes = 16;
    azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
    % is being formed. Rotates by (10 noUC + 3)*azStep degrees.
elseif (str == 'G') || (str == 'g')
    ucData = importdata('GaAsCoordinates.dat'); % (8 x 5) matrix
    % Brown for Ga, purple for As
    C = [0.395 0.215 0.215; 0.395 0.215 0.215; 0.395 0.215 0.215; 0.395 0.215 0.215; ...
        0.773 0.502 0.949; 0.773 0.502 0.949; 0.773 0.502 0.949; 0.773 0.502 0.949];
    noAtoms = 8;
    ucxStep = 1;
    ucyStep = 1;
    uczStep = 1;
    ucSize = 4.066; % u.c. lattice constant of GaAs
    youChoose = 1.00; % How much smaller than u.c. the molecule is
    sphereRes = 16;
    azStep = 0.7; % Causes camera to swing azimuthally around the crystal as it
    % is being formed. Rotates by (10 noUC + 3)*azStep degrees.
end

vid = VideoWriter('uc2crystal.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);

noUC = 4;

ii = 0;
az = -125;

pol = 20;

% Ensure uc within 1 x 1 x 1 volume of plotted space
ucMax = max(max(ucData(:,1:3)));
ucMin = min(min(ucData(:,1:3)));
scaleFac = youChoose/ucSize; % Reduces length of molecule along x-axis to 1.0*youChoose

[x,y,z] = sphere(140);

% Read unit-cell data
for i = 1:noAtoms
    ii = ii+1;
    rad = ucData(i,5)*scaleFac;
    a = ucData(i,1)*scaleFac;
    b = ucData(i,2)*scaleFac;
    c = ucData(i,3)*scaleFac;
    atomCol = C(i,:);
    % Create concatenated surface image of the unit cell
    HkB(ii) = surf(-noUC+rad*x+a,-noUC+rad*y+b,noUC+rad*z+c,...
        'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
        'none','FaceLighting','flat',... %'gouraud',...
        'DiffuseStrength',1,'Visible','off');
    hold on
end
set(HkB,'Visible','on')

lp = [-1 -0 0.];
light('Position',lp,'Style','infinite');

% Pan out in Ken Burns style
zoomOutFrames = 50;

for kenB = 0:zoomOutFrames
    hold off
    % frame limits increasing from centering on molecule to it being at the
    % corner of a crystal
    xlo = -noUC + min(ucData(:,1))/ucSize - max(ucData(:,5))/ucSize;
    xhi = -noUC + max(ucData(:,1))/ucSize + max(ucData(:,5))/ucSize + 2*noUC*kenB/zoomOutFrames;
    ylo = -noUC + min(ucData(:,2))/ucSize - max(ucData(:,5))/ucSize;
    yhi = -noUC + max(ucData(:,2))/ucSize + max(ucData(:,5))/ucSize + 2*noUC*kenB/zoomOutFrames;
    zlo = noUC + min(ucData(:,3))/ucSize - max(ucData(:,5))/ucSize - 2*noUC*kenB/zoomOutFrames;
    zhi = noUC + max(ucData(:,3))/ucSize + max(ucData(:,5))/ucSize;
    
    axis square
    axis equal
    axis off
    set(gca,'View',[az,pol]);
    set(gca, 'Projection','perspective');
    xlim([xlo xhi]);
    ylim([ylo yhi]);
    zlim([zlo zhi]);
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

[x,y,z] = sphere(sphereRes); % Reduce resoution of spheres now they're tiny
hold off
% Generate first line of unit cells
for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
    ii = 0;
    for i = 1:noAtoms % Loop to draw atoms of one unit cell
        ii = ii+1;
        rad = ucData(i,5)*scaleFac;
        a = ucData(i,1)*scaleFac;
        b = ucData(i,2)*scaleFac;
        c = ucData(i,3)*scaleFac;
        atomCol = C(i,:);
        H1d(ii) = surf(no1d+rad*x+a,-noUC+rad*y+b,noUC+rad*z+c,...
            'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
            'none','FaceLighting','flat',... %'gouraud',...
            'DiffuseStrength',1,'Visible','off');
        hold on
    end
    set(H1d,'Visible','on')
    
    axis square
    axis equal
    axis off
    
    
    set(gca, 'Projection','perspective');
    xlim([xlo xhi]);
    ylim([ylo yhi]);
    zlim([zlo zhi]);
    if (no1d == -noUC)
        light('Position',lp,'Style','infinite');
    end
    
    for jj = 1:6
        %az = az+azStep; % commented out if you don't want the line to
        %rotate
        set(gca,'View',[az,pol]);
        % Store the frame
        frame = getframe(gcf);
        writeVideo(vid,frame);
    end
end % End of drawing one line

hold off
% Generate first plane of unit cells
for no2d = -noUC:ucyStep:noUC % Loop to replicate molecule over a single plane
    ii = 0;
    for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
        for i = 1:noAtoms % Loop to draw atoms of one unit cell
            ii = ii+1;
            rad = ucData(i,5)*scaleFac;
            a = ucData(i,1)*scaleFac;
            b = ucData(i,2)*scaleFac;
            c = ucData(i,3)*scaleFac;
            atomCol = C(i,:);
            H2d(ii) = surf(no1d+rad*x+a,no2d+rad*y+b,noUC+rad*z+c,...
                'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
                'none','FaceLighting','flat',... %'gouraud',...
                'DiffuseStrength',1,'Visible','off');
            hold on
        end % End of drawing one molecule
        if (no1d == -noUC) && (no2d == -noUC)
            light('Position',lp,'Style','infinite');
        end
    end % End of drawing one line
    set(H2d,'Visible','on')
    
    axis square
    axis equal
    axis off
    set(gca, 'Projection','perspective');
    xlim([xlo xhi]);
    ylim([ylo yhi]);
    zlim([zlo zhi]);
    
    for jj = 1:6
        az = az+azStep;
        set(gca,'View',[az,pol]);
        % Store the frame
        frame = getframe(gcf);
        writeVideo(vid,frame);
    end
    
end % End of drawing one plane


hold off
newplot
% Generate multiple planes of unit cells
for no3d = noUC:-uczStep:-noUC % Loop to replicate planes
    newplot
    ii = 0;
    for no2d = -noUC:ucyStep:noUC % Loop to replicate molecule over a single plane
        for no1d = -noUC:ucxStep:noUC % Loop to replicate molecule along a single line
            for i = 1:noAtoms % Loop to draw atoms of one unit cell
                ii = ii+1;
                rad = ucData(i,5)*scaleFac;
                a = ucData(i,1)*scaleFac;
                b = ucData(i,2)*scaleFac;
                c = ucData(i,3)*scaleFac;
                atomCol = C(i,:);
                H3d(ii) = surf(no1d+rad*x+a,no2d+rad*y+b,no3d+rad*z+c,...
                    'FaceAlpha',1.0,'FaceColor', atomCol,'LineStyle',...
                    'none','FaceLighting','flat',... %'gouraud',...
                    'DiffuseStrength',1,'Visible','off');
                hold on
            end % End of drawing one molecule
            if (no1d == -noUC) && (no2d == -noUC) && (no3d == noUC)
                light('Position',lp,'Style','infinite');
            end
        end % End of drawing one line
        
    end % End of drawing one plane
    set(H3d,'Visible','on')
    set(0, 'DefaultFigureRenderer', 'opengl');
    axis square
    axis equal
    axis off
    set(gca, 'Projection','perspective');
    xlim([xlo xhi]);
    ylim([ylo yhi]);
    zlim([zlo zhi]);
    
    for jj = 1:6
        az = az+azStep;
        set(gca,'View',[az,pol]);
        % Store the frame
        drawnow
        frame = getframe(gcf);
        writeVideo(vid,frame);
    end
    clear H3d
end % End of drawing multiple planes

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

