gtBoundaryPoles.m

function [output, save_angles]=gtBoundaryPoles(list, boundaries_structure, crack_vol)

%read(+plot), on a standard stereographic triangle(001, 101, 111), the poles of each boundary in
%the coordinates of the relevant grains.


save_angles=[];
output=[];

load('parameters.mat');
symm = gtCrystGetSymmetryOperators(parameters.cryst.crystal_system);

if isempty(list)
  list=1:length(boundaries_structure);
end


for k=1:length(list) %just look at the possible twin boundaries
  i=list(k)
%for i=2:1987  
%k=i;

%skip "open" boundaries (grain-air)
  if boundaries_structure(i).grain1==0 || boundaries_structure(i).grain2==0 || boundaries_structure(i).count<4
    continue
    
% else
 elseif boundaries_structure(i).sigma==3 
%disp('sigma 3')
%     %to use only those from a crack
%     elseif boundaries_structure(i).crack8_count>4 %& (boundaries_structure(i).crack8_count/boundaries_structure(i).count)>0.2%/boundaries_list(i,4)>0.5
%     
%     i
%     test=crack_vol==i;
%      [x,y,z]=ind2sub(size(test), find(test));
%     [origin, a, d, normd]=lsplane([x y z]);
%     crack_facet=a;
%     crack_facet(3)=-crack_facet(3);
%     
%     %angle between crack facet and grain boundary
%     test_angle=acos(dot(crack_facet, boundaries_structure(i).gb_norm))
%     save_angles=[save_angles; test_angle];
%     
%     if test_angle<30*pi/180 || test_angle>150*pi/180

      
   gb_normal=boundaries_structure(i).gb_norm'; 
   
   %playing with orientations here for id11 data
  % gb_normal(1)=-gb_normal(1);
    %  gb_normal(2)=-gb_normal(2);
   
   
   
    % Get the grain orientations
    R1 = boundaries_structure(i).grain1_R_vector;
    g1 = gtMathsRod2RotMat(R1);

    R2 = boundaries_structure(i).grain2_R_vector;
    g2 = gtMathsRod2RotMat(R2);
    
    gb_normal1 = inv(g1)*gb_normal;
    gb_normal2 = inv(g2)*gb_normal;
% 
% %should be able to read straight from the structure when not faffing
%   gb_normal1=boundaries_structure(i).gb_norm_grain1';
%   gb_normal2=boundaries_structure(i).gb_norm_grain2';
 
    % Need to search symmetry equivalents for the pole in the SST
    tmp1=[];tmp2=[];

    for ii=1:length(symm)
      tmp1(ii,:) = gb_normal1'*symm(ii).g3;
      tmp2(ii,:) = gb_normal2'*symm(ii).g3;
    end

%  output=[output; tmp1; tmp2];
  output = [ output                          ; ...
             repmat(i, length(symm), 1) tmp1 ; ...
             repmat(i, length(symm), 1) tmp2 ];

%    end %for crack condition
  
  end

  
end
  %end
  
%plot the poles

if 0 
figure(1)
a=sqrt(output(1:2:(end-1),1).^2 + output(1:2:(end-1),2).^2);
b=atan2(output(1:2:(end-1),2),output(1:2:(end-1),1));
plot(a.*cos(b), a.*sin(b), 'o');

a=sqrt(output(2:2:end,1).^2 + output(2:2:end,2).^2);
b=atan2(output(2:2:end,2),output(2:2:end,1));
plot(a.*cos(b), a.*sin(b), 'ro');

figure(2)
plot3(output(1:2:(end-1),1), output(1:2:(end-1),2), output(1:2:(end-1),3), 'bo')

plot3(output(2:2:end,1), output(2:2:end,2), output(2:2:end,3), 'ro')
end