gtTwinTest.m

function out = gtTwinTest(phaseid, parameters, grain1, grain2, grain, vol_out, id)
% out = gtTwinTest(phaseid, parameters, grain1, grain2, grain, vol_out, id)
% grain1 = grainID for grain 1
% grain2 = grainID for grain 2
% grain  = grain structure <cell>
% vol_out = output from gt_find_grain_boundaries

crystal_system = parameters.cryst(phaseid).crystal_system; 
spacegroup     = parameters.cryst(phaseid).spacegroup;
latticepar     = parameters.cryst(phaseid).latticepar;

% 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 symmetry operators
symm = gtCrystGetSymmetryOperators(crystal_system, spacegroup); 

if grain1==0
    g1=[];
    g1equiv=[];
elseif grain1~=-1 && grain1~=grain2
    R1 = r_vectors(grain1,2:4);
    g1 = gtMathsRod2RotMat(R1);
else
    out = [];
    return;
end

if grain2==0
    g2 = [];
elseif grain2~=-1 && grain2~=grain1
    R2 = r_vectors(grain2,2:4);
    g2 = gtMathsRod2RotMat(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

% same calculation but using hexagonal symm operations
% calculate the misorientation angle and axis between these grains
% TODO: replace by a MISORIENTATION function
if ~isempty(g1) && ~isempty(g2)
    % if we have both grains
    % Need to search symmetry equivalents for the minimum misorientation
    rot_offset=[];
    % warning('g*symm.g maybe should be symm.g*g')
    for jj=1:length(symm)
        g1equiv = g1*symm(jj).g3;
        netg = g1equiv\g2;
        rot_offset(jj) = acos((trace(netg)-1)/2);
    end

    dummy = find(rot_offset ==min(rot_offset));
    g1equiv = g1*symm(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    


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';

out.mis_axis_hex    = gtCart2Hex(misorientation_axis', latticepar);
out.symm            = symm;


out.center1         = grain{grain1}.center;
out.center2         = grain{grain2}.center;
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.mis_axis_hex out.mis_axis out.grain1 out.grain2];
out.table_header    = 'mis_angle  mis_axis_hex mis_axis grain1 grain2';

if ~isempty(vol_out)
    %even if we have only one grain, can deal with the boundary plane
    test = vol_out==id;
    [x,y,z]=ind2sub(size(test), find(test));
    [origin, a]=lsplane([x y z]);


    gb_normal=a;
   % gb_normal(3)=-gb_normal(3);

    out.gb_norm = gb_normal';
    out.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
        out.gb_norm_grain1=plane1';
        out.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

        out.gb_norm_grain1=plane1';
        out.gb_norm_grain2=plane2';

        if ~isempty(plane1)
            hkl=plane1';
            hkil = gtCart2Hex(hkl, parameters.cryst(phaseid).latticepar);
            out.gb_norm_hex_grain1=hkil;
        else
            out.gb_norm_hex_grain1=[];
        end
        if ~isempty(plane2)
            hkl=plane2';
            hkil = gtCart2Hex(hkl, parameters.cryst(phaseid).latticepar);
            out.gb_norm_hex_grain2=hkil;
        else
            out.gb_norm_hex_grain2=[];
        end

    else
        disp('crystallography not supported!')
    end
end

end