Skip to content
Snippets Groups Projects
Commit 5624d9ad authored by Laura Nervo's avatar Laura Nervo
Browse files

remove duplicate function gtTwinTest


Signed-off-by: default avatarLaura Nervo <laura.nervo@esrf.fr>
parent d58ae665
No related merge requests found
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 = gtMathsRod2OriMat(R1);
else
out = [];
return;
end
if grain2==0
g2 = [];
elseif grain2~=-1 && grain2~=grain1
R2 = r_vectors(grain2,2:4);
g2 = gtMathsRod2OriMat(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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment