gtReadBoundaryProperties.m

function gtReadBoundaryProperties(phaseid, boundaries_structure, vol_out, r_vectors)
% gtReadBoundaryProperties(phaseid, boundaries_structure, vol_out, r_vectors)

disp('have changed this to axis convention matlab1 // instX , matlab2 // instY, matlab3 // instZ')
disp('AK, 11/8/2010')
pause(3)

%
%continue to add to boundaries_structure.mat - orientations,
%disorientation, r-vectors etc.
parameters=[];
load('parameters.mat');
spacegroup=parameters.cryst(phaseid).spacegroup;

%find the max grainID - ugly but ca va
maxgrain=0;
for i=1:length(boundaries_structure)
    maxgrain=max(maxgrain, boundaries_structure(i).grain1);
    maxgrain=max(maxgrain, boundaries_structure(i).grain2);
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...')
    return
end


% loop through the boundaries
for i=1:size(boundaries_structure, 2)

    if boundaries_structure(i).count<4 %can still be useful to calc index planes where we have only one of the grains
        continue
    else

        %Read in the r-vector data for the two grains
        i
        grain1=boundaries_structure(i).grain1;
        grain2=boundaries_structure(i).grain2;

        if grain1==0
            g1=[];
            g1equiv=[];
        elseif grain1==-1
            continue
        else
            R1= r_vectors(grain1,2:4);
            boundaries_structure(i).grain1_R_vector=R1;
            g1=Rod2g(R1);
        end

        if grain2==0
            g2=[];
        elseif grain2==-1
            continue
        else
            R2= r_vectors(grain2,2:4);
            boundaries_structure(i).grain2_R_vector=R2;
            g2=Rod2g(R2);
        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
            
            boundaries_structure(i).misorientation_axis=misorientation_axis';
            boundaries_structure(i).misorientation_angle=misorientation_angle;
        else
            boundaries_structure(i).misorientation_axis=[];
            boundaries_structure(i).misorientation_angle=[];
        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);
                    boundaries_structure(i).eigVal=eigVal(p,p);
                end
            end
            
            boundaries_structure(i).misorientation_axis=misorientation_axis';
            boundaries_structure(i).misorientation_angle=misorientation_angle;
            
        else
            boundaries_structure(i).misorientation_axis=[];
            boundaries_structure(i).misorientation_angle=[];
            boundaries_structure(i).eigVal=[];
        end
    
    
end


        %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);

        boundaries_structure(i).gb_norm = gb_normal';
        boundaries_structure(i).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
            boundaries_structure(i).gb_norm_grain1=plane1';
            boundaries_structure(i).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

            boundaries_structure(i).gb_norm_grain1=plane1';
            boundaries_structure(i).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.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));

                boundaries_structure(i).gb_norm_hex_grain1=hkil;
            else
                boundaries_structure(i).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.latticepar(1)));
                hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));
                boundaries_structure(i).gb_norm_hex_grain2=hkil;
            else
                boundaries_structure(i).gb_norm_hex_grain2=[];
            end
% 
%             if ~isempty(boundaries_structure(i).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.latticepar(1)));
%                 hkil(4)   = hkil(4) / (1/parameters.cryst.latticepar(3));
%                 boundaries_structure(i).misorientation_axis_hex=hkil;
%             end

        else
            disp('crystallography not supported!')
        end

    end % skip small boundaries

end

save('boundaries_structure.mat', 'boundaries_structure')