% The following code enables the user to visualize a crystal structure
% with an underlying  cubic lattice in MATLAB. Once MATLAB reads the input
% file, and it creates the coordinates of the
% Primitive Unit Cell (PUC) of the crystal. Upon doing this, MATLAB then
% expands these coordinates into the desired dimensions, as specified in
% the input file to create the full crystal lattice.
% Author: Rajiv Uttamchandani
% Institution: California State University, Northridge


% CODE BEGINS HERE

% Taking input from file
clear all;
filename = input('Enter the file name of the input file: ', 's');
Title1 = input('Enter the title for graph of entire crystal: ', 's');
Title2 = input('Enter the title for graph of the Primitive Unit Cell: ', 
's');
fid = fopen(filename,'r');
H = fscanf(fid,'%g',[3 inf]);
F =H';

[rx,ry,rz] = sphere(30);
rx = rx/5; ry = ry/5; rz = rz/5;

% Creating PUC
count = 1;
for i=1:2
    for j=1:2
        for k=1:2
             T(count,:) = i*F(3,:)+j*F(4,:)+k*F(5,:);
             count = count +1;
        end
    end
end

% Additional units for the basis
for i=1:F(2,1)
    T(count,:) = F(i+5,:);
    count = count+1;
end

N = length(T);

% The prodcedure below enables the translation of the PUC.

% Translation in the a1-axis
count = 1;

for j = 1:N
        M(count,:) = T(j,:);
        count = count+1;
end

for i = 1:F(1,1)
    for j =1:N
      M(count,:) = T(j,:) + i*F(3,:);
      count = count+1;
    end
end

m = length(M);

% Translation in the a2-axis
for j = 1:F(1,2)
    for i = 1:m
        M(count,:) = M(i,:) + j*F(4,:);
        count = count + 1;
    end
end

m = length(M);

% Translation in the a3-axis
for j = 1:F(1,3)
    for i = 1:m
        M(count,:) = M(i,:) + j*F(5,:);
        count = count + 1;
    end
end

% Plotting the points (basis coordinates)
figure(1)
for i=1:length(M)
    surf(rx+M(i,1), ry+M(i,2), 
rz+M(i,3),'FaceColor','red','EdgeColor','none')
    hold on;
end
camlight left; lighting phong;
xlabel('x'), ylabel('y'), zlabel('z'), title(Title1)

% Plotting the points (Bravais Lattice Coordinates)
figure(2)
for i=1:length(T)
    surf(rx+T(i,1), ry+T(i,2), 
rz+T(i,3),'FaceColor','green','EdgeColor','none')
    hold on;
end
camlight left; lighting phong;
xlabel('x'), ylabel('y'), zlabel('z'), title(Title2)