function out = gtTwinTest(phaseid, parameters, grain1, grain2, grain)
% out = gtTwinTest(phaseid, parameters, grain1, grain2, grain)
pxsize = (parameters.labgeo.pixelsizeu + parameters.labgeo.pixelsizev)/2;
spacegroup = parameters.cryst(phaseid).spacegroup;
% find the max grainID
maxgrain = max(0, grain1);
maxgrain = max(maxgrain, grain2);
r_vectors = [];
for ii=1:length(grain)
r_vectors(ii, 1) = ii;
r_vectors(ii, 2:4) = grain{ii}.R_vector;
end
% get sym operators
if spacegroup==225 || spacegroup==229
% cubic material
sym = gtGetCubicSymOp();
elseif spacegroup==167 || spacegroup==663 || spacegroup==194
% hexagonal material - 4x4 element sym operators
sym = gtGetHexagonalSymOp_sab();
else
disp('unrecognised spacegroup! quitting...')
out = [];
return
end
if grain1==0
g1=[];
g1equiv=[];
elseif grain1~=-1 && grain1~=grain2
R1 = r_vectors(grain1,2:4);
g1 = Rod2g(R1);
else
out = [];
return;
end
if grain2==0
g2 = [];
elseif grain2~=-1 && grain2~=grain1
R2 = r_vectors(grain2,2:4);
g2 = Rod2g(R2);
else
out = [];
return;
end
% this only works for cubic systems until someone works out how to pdo te
% permutations of the g1 for the hexagonal system
if spacegroup==225 || spacegroup==229
% calculate the misorientation angle and axis between these grains
if ~isempty(g1) && ~isempty(g2)
%if we have both grains
%need to search symmetry equivelents for the minimum misorientation
rot_offset=[];
% warning('g*sym.g maybe should be sym.g*g')
for j=1:24
g1equiv = g1*sym(j).g;
netg = g1equiv\g2;
rot_offset(j) = acos((trace(netg)-1)/2);
end
dummy = find(rot_offset ==min(rot_offset));
g1equiv = g1*sym(dummy).g;
netg = g1equiv\g2;
misorientation_angle = acosd((trace(netg)-1)/2);
[eigVec, eigVal]=eig(netg);
for p=1:3
if isreal(eigVec(:,p)) %find the real eigenvector - should be just one
misorientation_axis=eigVec(:,p);
end
end
end
elseif spacegroup==167 || spacegroup==663 || spacegroup==194
% same calculation but using hexagonal sym operations
% calculate the misorientation angle and axis between these grains
if ~isempty(g1) && ~isempty(g2)
%if we have both grains
%need to search symmetry equivelents for the minimum misorientation
rot_offset=[];
% warning('g*sym.g maybe should be sym.g*g')
for j=1:length(sym)
g1equiv = g1*sym(j).g3;
netg = g1equiv\g2;
rot_offset(j) = acos((trace(netg)-1)/2);
end
dummy = find(rot_offset ==min(rot_offset));
g1equiv = g1*sym(dummy).g3;
netg = g1equiv\g2;
misorientation_angle = acosd((trace(netg)-1)/2);
[eigVec, eigVal]=eig(netg);
for p=1:3
if isreal(eigVec(:,p)) %find the real eigenvector - should be just one
misorientation_axis=eigVec(:,p);
eigVal=eigVal(p,p);
end
end
end
end
out.grain1 = grain1;
out.grain2 = grain2;
out.grain1_R_vector = r_vectors(grain1,2:4);
out.grain2_R_vector = r_vectors(grain2,2:4);
out.mis_angle = misorientation_angle;
out.mis_axis = misorientation_axis';
hkl=misorientation_axis';
hkil=[];
hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
hkil(2) = (2/sqrt(3))*hkl(2);
hkil(3) = -(hkil(1) + hkil(2));
hkil(4) = hkl(3);
% account for lattice vector lengths
hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
out.mis_axis_hex = hkil;
out.sym = sym;
out.angle = rad2deg(atan2(norm(cross(R1,R2)), dot(R1,R2)));
out.center1 = grain{grain1}.center/pxsize;
out.center2 = grain{grain2}.center/pxsize;
out.dist_centers = norm(out.center1-out.center2);
out.radius1 = min(grain{grain1}.stat.bbxsmean,grain{grain1}.stat.bbysmean);
out.radius2 = min(grain{grain2}.stat.bbxsmean,grain{grain2}.stat.bbysmean);
out.distf = parameters.index.strategy.m.end.distf;
out.is_inside1 = (out.dist_centers < out.distf*out.radius1);
out.is_inside2 = (out.dist_centers < out.distf*out.radius2);
gtDBConnect();
[int, x, y] = mym(['select avintint, avbbXsize, avbbYsize from ' parameters.acq.pair_tablename ' where grainID = ',num2str(grain1)]);
out.avgint1 = mean(int);
[int, x, y] = mym(['select avintint, avbbXsize, avbbYsize from ' parameters.acq.pair_tablename ' where grainID = ',num2str(grain2)]);
out.avgint2 = mean(int);
out.ratio1 = out.dist_centers/(out.radius1*out.distf);
out.ratio2 = out.dist_centers/(out.radius2*out.distf);
out.table = [out.mis_angle out.dist_centers out.distf out.distf*out.radius1 out.distf*out.radius2];
out.table_header = 'mis_angle dist_centers distf distf*radius1 distf*radius2';
%
% %even if we have only one grain, can deal with the boundary plane
% test = vol_out==i;
% [x,y,z]=ind2sub(size(test), find(test));
% [origin, a]=lsplane([x y z]);
%
%
% gb_normal=a;
% % gb_normal(3)=-gb_normal(3);
%
% gb_norm = gb_normal';
% gb_centre = origin';
%
% if ~isempty(g1equiv)
% plane1=g1equiv\gb_normal;
% elseif ~isempty(g1)
% plane1=g1\gb_normal;
% else
% plane1=[];
% end
%
% if ~isempty(g2)
% plane2=g2\gb_normal;
% else
% plane2=[];
% end
%
% if spacegroup==225 || spacegroup==229
% % if cubic then this is the index plane
% gb_norm_grain1=plane1';
% gb_norm_grain2=plane2';
%
% elseif spacegroup==167 || spacegroup==663 || spacegroup==194
% % if hexagonal, keep cartesian vector, but also
% % convert the cartesian plane normal to hexagonal four index notation
% % ditto for misorientation axis
%
% gb_norm_grain1=plane1';
% gb_norm_grain2=plane2';
%
% if ~isempty(plane1)
% hkl=plane1;
% hkil=[];
% hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
% hkil(2) = (2/sqrt(3))*hkl(2);
% hkil(3) = -(hkil(1) + hkil(2));
% hkil(4) = hkl(3);
% % account for lattice vector lengths
% hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
% hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
%
% gb_norm_hex_grain1=hkil;
% else
% gb_norm_hex_grain1=[];
% end
%
% if ~isempty(plane2)
% hkl=plane2;
%
% hkil=[];
% hkil(1) = hkl(1) - (1/sqrt(3))*hkl(2);
% hkil(2) = (2/sqrt(3))*hkl(2);
% hkil(3) = -(hkil(1) + hkil(2));
% hkil(4) = hkl(3);
% % account for lattice vector lengths
% hkil(1:3) = hkil(1:3) / (2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1)));
% hkil(4) = hkil(4) / (1/parameters.cryst(phaseid).latticepar(3));
% gb_norm_hex_grain2=hkil;
% else
% gb_norm_hex_grain2=[];
% end
%
%
%
% else
% disp('crystallography not supported!')
% end
end