% Script lattice constants LAMMPS
% Ambarish Kulkarni
clc
clear

mode = 2; %1: abc to lammps, 2: lammps to abc

if mode == 2
    % Used if mode = 2
    data_raw{1} = [0.00       0.00       0.00]'; % xlo_bound ylo_bound zlo_bound
    data_raw{2} = [2.2856     1.0658     0.7050]'; % xhi_bound yhi_bound zhi_bound
    data_raw{3} = [0.6428-1     0.6428     0.2997]'; % xy xz yz
    
%     data_raw{1} = [0.00       0.00       0.00]'; % xlo_bound ylo_bound zlo_bound
%     data_raw{2} = [1.0000     1.6428     0.7660]'; % xhi_bound yhi_bound zhi_bound
%     data_raw{3} = [0.0000     0.0000     0.6428-1]'; % xy xz yz
elseif mode == 1
    % Used if mode = 1
    a =1;
    b =1;
    c =1;
    alpha_deg = 50; % deg
    beta_deg = 50;
    gamma_deg = 50;
else
    error('Check mode\n');
end

if mode == 2
    xlo_bound = data_raw{1}(1);
    ylo_bound = data_raw{1}(2);
    zlo_bound = data_raw{1}(3);
    
    xhi_bound = data_raw{2}(1);
    yhi_bound = data_raw{2}(2);
    zhi_bound = data_raw{2}(3);
    
    xy = data_raw{3}(1);
    xz = data_raw{3}(2);
    yz = data_raw{3}(3);
    
    xlo = xlo_bound - min([0.0,xy,xz,xy+xz]);
    xhi= xhi_bound - max([0.0,xy,xz,xy+xz]);
    ylo = ylo_bound - min([0.0,yz]);
    yhi = yhi_bound - max([0.0,yz]);
    zlo = zlo_bound;
    zhi = zhi_bound;
    
    lx = xhi - xlo;
    ly = yhi - ylo;
    lz = zhi - zlo;
    
    a = lx;
    b = sqrt(ly^2 + xy^2);
    c = sqrt(lz^2 + xz^2 +yz^2);
    
    alpha = acos((xy*xz + ly*yz)/b/c);
    beta = acos(xz/c);
    gamma = acos(xy/b);
    
    alpha_deg = alpha * 180/pi;
    beta_deg = beta * 180/pi;
    gamma_deg = gamma * 180/pi;
    
    fprintf('Lattice constants = \n');
    fprintf('%10s %10s %10s %10s %10s %10s\n','a','b','c', 'alpha_deg','beta_deg','gamma_deg');
    fprintf('%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',a,b,c, alpha_deg, beta_deg, gamma_deg);
    
elseif mode == 1
    alpha = alpha_deg * pi/180;
    beta = beta_deg * pi/180;
    gamma = gamma_deg * pi/180;
    xy = b*cos(gamma);
    lx = a;
    xz = c*cos(beta);
    ly = sqrt(b^2 - xy^2);
    yz = (b*c*cos(alpha)- xy*xz)/ly;
    lz = sqrt(c^2 - xz^2 - yz^2);
    xlo =0;
    ylo =0;
    zlo = 0;
    xhi = lx;
    yhi = ly;
    zhi = lz;
    
    xlo_bound = xlo + min([0.0,xy,xz,xy+xz]);
    xhi_bound = xhi + max([0.0,xy,xz,xy+xz]);
    ylo_bound = ylo + min([0.0,yz]);
    yhi_bound = yhi + max([0.0,yz]);
    zlo_bound = zlo;
    zhi_bound = zhi; 
    
    fprintf('Lammps lo_hi = \n');
    fprintf('%10s %10s %10s %10s %10s %10s\n','xlo','ylo','zlo', 'xhi','yhi','zhi');
    fprintf('%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',xlo,ylo,zlo, xhi,yhi,zhi);
    fprintf('\n');
    
    fprintf('Lammps tits = \n');
    fprintf('%10s %10s %10s \n','xy','xz','yz');
    fprintf('%10.4f %10.4f %10.4f\n',xy,xz,yz);
    fprintf('\n');
    
    fprintf('Lammps bounds = \n');
    fprintf('%10s %10s %10s %10s %10s %10s\n','xlo_bound','ylo_bound','zlo_bound', 'xhi_bound','yhi_bound','zhi_bound');
    fprintf('%10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n',xlo_bound,ylo_bound,zlo_bound, xhi_bound,yhi_bound,zhi_bound);
    fprintf('\n');
else
    error('Check mode\n');
end

