% This script takes in input the average orientation of a grain expressed
% in Euler angles (or directly in miller indexes) and returns the
% coordinates for the 2D projection defined in Daviddi et al, Chem. Sci.
% 2020

%% Input
% Define low index grains
g1=[0 0 1];
g2=[0 1 1];
g3=[1 1 1];
% Input Euler angles, they should be in a 3 column array, with each column
% being a different Euler angle(conventionally ordered) and each row
% representing a grain; It will be called "euler"
%% Calculation of Coords of low index planes 
g=[g1;g2;g3];
P=zeros(3,3);
for i=1:3
    for j=1:3
        if i==j
        P(i,j)=0;
        else
        clear ang
ang=(g(i,1)*g(j,1)+g(i,2)*g(j,2)+g(i,3)*g(j,3))/(sqrt(g(i,1)^2+g(i,2)^2+g(i,3)^2)*sqrt(g(j,1)^2+g(j,2)^2+g(j,3)^2));
ang=acos(ang);
P(i,j)=real(rad2deg(ang));
        end
    end
end
tf = issymmetric(P);
if tf ==0
    disp ('The low index matrix i not symmetric, this calculation cannot work');
    return
else
end
a=P(1,2);
b=P(1,3);
c=P(2,3);
%% Calculation of parameters
u=a*c+b*c-c^2;
v=c*b+a*b-b^2;
w=a*c+a*b-a^2;
k=-2*a*b*c;

mu1=sqrt(c^2+b^2+2*a^2-2*b*c);
mu2=sqrt(u^2+v^2+w^2);

e11=a/mu1;
e12=-a/mu1;
e13=(c-b)/mu1;
e21=(w*(c-b)+v*a)/(mu1*mu2);
e22=(-u*(c-b)+v*a)/(mu1*mu2);
e23=-(a*(u+w))/(mu1*mu2);

%% Calculation of hkl indexes from Euler angles
E=euler(:,2:3);
[siz,~]=size(E);
hkl=zeros(siz,3);
hkl2=zeros(siz,3);
for l=1:siz
 hkl(l,1)=sin(deg2rad(E(l,1)))*sin(deg2rad(E(l,2)));  
 hkl(l,2)=sin(deg2rad(E(l,1)))*cos(deg2rad(E(l,2)));  
 hkl(l,3)=cos(deg2rad(E(l,1)));
 hkl2(l,:)=sort(hkl(l,:));
end

PE=zeros(siz,3);

for l=1:siz
 
PE(l,1)=(hkl2(l,1)*g(1,1)+hkl2(l,2)*g(1,2)+hkl2(l,3)*g(1,3))/(sqrt(hkl2(l,1)^2+hkl2(l,2)^2+hkl2(l,3)^2)*sqrt(g(1,1)^2+g(1,2)^2+g(1,3)^2));
PE(l,2)=(hkl2(l,1)*g(2,1)+hkl2(l,2)*g(2,2)+hkl2(l,3)*g(2,3))/(sqrt(hkl2(l,1)^2+hkl2(l,2)^2+hkl2(l,3)^2)*sqrt(g(2,1)^2+g(2,2)^2+g(2,3)^2));
PE(l,3)=(hkl2(l,1)*g(3,1)+hkl2(l,2)*g(3,2)+hkl2(l,3)*g(3,3))/(sqrt(hkl2(l,1)^2+hkl2(l,2)^2+hkl2(l,3)^2)*sqrt(g(3,1)^2+g(3,2)^2+g(3,3)^2));

PE(l,1)=real(rad2deg(acos(PE(l,1))));
PE(l,2)=real(rad2deg(acos(PE(l,2))));
PE(l,3)=real(rad2deg(acos(PE(l,3))));
end

C1=e11*PE(:,1)+e12*(PE(:,2)-a)+e13*(PE(:,3)-b);
C2=e21*PE(:,1)+e22*(PE(:,2)-a)+e23*(PE(:,3)-b);
lowC1=e11*P(:,1)+e12*(P(:,2)-a)+e13*(P(:,3)-b);
lowC2=e21*P(:,1)+e22*(P(:,2)-a)+e23*(P(:,3)-b);

C1_C2=[C1,C2];
C1_C2_low_index=[lowC1,lowC2];

dlmwrite('C1_C2.tsv',C1_C2)
dlmwrite('C1_C2low_index.tsv',C1_C2_low_index);