gtFindPolesInRange.m

function [uvw_poles, uvw_weights, max_phi, max_psi] = gtFindPolesInRange(poles, weights, range, mark_poles)
% #
% GTFINDPOLESINRANGE  Finds the poles in the right phi/psi range basing on the crystal symmetry.
%     [uvw_poles, uvw_weights, max_phi, max_psi] = gtFindPolesInRange(poles, weights, range, mark_poles)
%     --------------------------------------------------------------------------------------------------
%     Given a phi/psi range, it can be possible to mark/draw the poles in the pole figure. 
%     Options for "mark_poles" are main/bcc/all for cubic, main/all for hexagonal. 
%     Only cubic and hexagonal crystals are taken into account.
%
%     INPUT:
%       poles       = <double>   If given, only the poles in the specified region 
%                                are chosen; if empty, the theoretical poles are 
%                                calculated basing on the symmetry.
%       weigths     = <double>   Weights of the poles (can be empty)
%       range       = <string>   Phi range for the pole figure
%                                'phi360','phi90','sst','hex'
%       mark_poles  = <string>   Option to draw poles in the pole figure. 
%                                Possibilities are: 'main','all' (for cubic also 'bcc')
%
%     OUTPUT:
%       uvw_poles   = <double>   Poles <uvw> in the phi/psi given range
%       uvw_weigths = <double>   Weights of the poles in phi/psi given range
%
%     Version 001 17-01-2013 by LNervo
% #



if isempty(poles) % finds poles [hkl] and convert them into the cartesian space
    parameters = [];
    load('parameters.mat');
    cryst = parameters.cryst;
    clear parameters
    spacegroup     = cryst.spacegroup;
    latticepar     = cryst.latticepar;
    crystal_system = cryst.crystal_system;
    
    if strcmp(crystal_system, 'cubic')
        if strcmp(mark_poles, 'main')
            poles = [0 0 1; 1 0 1; 1 1 1];
        elseif strcmp(mark_poles, 'bcc')
            poles = [0 0 1; 1 0 1; 1 1 1; 1 1 2; 2 1 3];
        elseif strcmp(mark_poles, 'all')
            poles = [0 0 1; 1 0 1; 1 1 1;...
                     2 0 1; 2 1 1; 2 2 1;...
                     3 0 1; 3 1 1; 3 2 0; 3 2 1; 3 2 2; 3 3 1; 3 3 2;...
                     4 0 1; 4 1 1; 4 2 1; 4 3 0; 4 3 1; 4 3 2; 4 3 3; 4 4 1; 4 4 3];
        else
            gtError('gtFindPolesInRange:wrongrgument','unrecognized "mark_poles" option - should be all/main/bcc')
        end
    elseif strcmp(crystal_system, 'hexagonal')
        if strcmp(mark_poles,'main')
            poles = [0 0 0 1; 1 0 -1 0; 2 -1 -1 0];
        elseif strcmp(mark_poles, 'all')
            poles = [0 0 0 1; 1 0 -1 0; 2 -1 -1 0; 1 0 -1 1; 2 -1 -1 1];
        else
            gtError('gtFindPolesInRange:wrongrgument','unsupported "mark_poles" option - should be all/main')
        end
    else
        gtError('gtFindPolesInRange:wrongArgument','The crystal system is not recognised...Quitting')
    end

    % permute
    poles_2 = [];
    for ii = 1:size(poles,1)
        poles_2 = [poles_2; gtCrystSignedHKLs(poles(ii,:), gtCrystGetSymmetryOperators(crystal_system, spacegroup))];
    end

    % convert to cartesian
    poles = gtCrystHKL2Cartesian(poles_2', gtCrystHKL2CartesianMatrix(latticepar))';
    poles = [poles; -poles];
end

if size(poles,2) == 4 % first column should be grainID in this case; removed
    poles(:,1) = [];
end

if isempty(weights) || size(weights,1) ~= size(poles,1)
    weights = ones(size(poles,1),1);
end

% converts poles from cartesian to spherical coordinates
[phi_uvw,psi_uvw,r_uvw]=cart2sph(poles(:,1),poles(:,2),poles(:,3));
% bring phi in [0,2pi] range from [-pi,+pi] range
phi_uvw=phi_uvw+pi;
% bring psi in [0,pi] range from [-pi/2,pi/2] range
psi_uvw=pi/2-psi_uvw;

% find the right phi-psi range: allowed are: 0 <= phi <= 2pi; 0 <= psi <= pi
if strcmp(range, 'sst')
    % for cubic standard stereographic triangle
    phimax=pi/4;
    psimax=atan(sqrt(2));
elseif strcmp(range, 'phi90')
    % for a 90 degree section of the pole figure
    phimax=pi/2;
    psimax=pi/2;
elseif strcmp(range, 'phi360')
    % for 360 pole figure
    phimax=2*pi;
    psimax=pi/2;
elseif strcmp(range, 'hex')
    % for hexagonal sst
    phimax=pi/6;
    psimax=pi/2;
else
    gtError('gtFindPolesInRange:wrongArgument','Unrecognised plot range...Quitting')
end

dum=find(phi_uvw <= phimax & phi_uvw >= 0 & psi_uvw <= psimax & psi_uvw >= 0);
uvw_poles = poles(dum, :);

if nargout > 1
    uvw_weights = weights(dum);
    if nargout > 2
        max_phi = phimax;
        if nargout > 3
            max_psi = psimax;
        end
    end
end

end % end of function