gtTwinTest.m

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