function [hfig, save_density, valid_poles_weight] = gtMakePoleFigure(varargin)
% GTMAKEPOLEFIGURE  Make a pole figure given a list of poles
%     [hfig, save_density, valid_poles_weight] = gtMakePoleFigure({varargin})
%     -----------------------------------------------------------------------
%
%     the poles should cover a larger angular area than required, to avoid
%     depleted pole density at the edges.  Any un-needed poles will be crops to
%     improve speed.
%
%     can pass in either poles or a density matrix.  If the density matrix, YOU
%     have to make sure the other variables are consistant with it's size.
%     if passing in density also need to pass valid_poles_weight for
%     normalising
%
%     Loads parameters.mat:cryst(phaseid):spacegroup,latticepar.
%
%     INPUT: parse parameter pairs format.
%
%       'phaseid'     = phase number <int> {1}
%
%       'poles'       = poles from gtReadBoundariesStructure
%                       [poles, weigths] = gtReadBoundariesStructure(boundaries_structure, vThr, DoSym)
%
%       'range'       = {'sst'}/'phi90'/'phi360'
%                       standard stereographic triangle, or quarter or whole hemisphere.
%
%       'label_poles' = 'all'/'bcc'/'main'/{'none'}
%                       all hkl poles or just the corners of the sst.
%
%       'mark_poles'  = {'all'}/'bcc'/'main'/'none'
%
%       'weights'     = can supply a vector the same length as poles to weight the
%                       pole data {[]}
%
%       'searchR'     = smoothing radius in pole figure {0.1 rad}
%
%       'increment'   = data point density in pole figure {0.01 rad}
%
%       'plot_figure' = can choose not to plot figure, and use the returned
%                       density values (non normalised) for something
%                       {true}
%
%     Version 002 26-06-2012 by LNervo
%       Add phaseid as input argument for multiple phases materials
%
%     Version 001 Sep-2009 by AKing
%       could be improved to pass in a list of which hkl poles to plot and to
%       label.


% read in the specs
appdefaults.phaseid            = 1;
appdefaults.poles              = [];
appdefaults.density            = [];
appdefaults.valid_poles_weight = [];
appdefaults.weights            = [];
appdefaults.range              = 'sst'; %or phi90, phi360
appdefaults.label_poles        = 'none';
appdefaults.mark_poles         = 'all';
appdefaults.increment          = 0.01;
appdefaults.searchR            = 0.1;
appdefaults.new_figure         = true;
appdefaults.plot_figure        = true;

app = parse_pv_pairs(appdefaults,varargin);

phaseid            = app.phaseid;
poles              = app.poles;
density            = app.density;
searchR            = app.searchR;
inc                = app.increment;
weights            = app.weights;
valid_poles_weight = app.valid_poles_weight;

% get the parameters file

parameters = [];
load('parameters.mat');
spacegroup = parameters.cryst(phaseid).spacegroup;

if ~isempty(poles)
    
    % if no weighting data passed in
    if isempty(weights)
        weights=ones(size(poles, 1),1);
    end
    
    % before building figure, take only those poles that fall in the range 0-phimax,
    % 0-psimax (allowing for seach radius)
    if strcmp(app.range, 'sst')
        % for standard stereographic triangle
        phimax=pi/4;
        psimax=atan(sqrt(2));
        dum=find(poles(:,2)>-sin(searchR) & poles(:,2)<poles(:,1)+(sqrt(2)*sin(searchR)) & poles(:,3)>poles(:,1)-(sqrt(2)*sin(searchR)));
        poles=poles(dum, :);
        weights=weights(dum);
        valid_poles_weight=sum(weights(find(poles(:,2)>=0 & poles(:,2)<=poles(:,1) & poles(:,3)>=poles(:,1))));
        total_solid_angle=pi/12;
    elseif strcmp(app.range, 'phi90')
        % for a 90 degree section of the pole figure
        phimax=pi/2;
        psimax=pi/2;
        dum=find(poles(:,1)>-sin(searchR) & poles(:,2)>-sin(searchR) & poles(:,3)>cos(psimax)-sin(searchR));
        poles=poles(dum, :);
        weights=weights(dum);
        valid_poles_weight=sum(weights(find(poles(:,1)>=0 & poles(:,2)>=0 & poles(:,3)>=0)));
        total_solid_angle=pi/2;
    elseif strcmp(app.range, 'phi360')
        % for 360 pole figure - reduce only using psi
        phimax=2*pi;
        psimax=pi/2;
        dum=find(poles(:,3)>cos(psimax)-sin(searchR));
        poles=poles(dum,:);
        weights=weights(dum);
        valid_poles_weight=sum(weights(find(poles(:,3)>=0)));
        total_solid_angle=2*pi;
    else
        disp('unrecognised plot range :-(')
        return
    end
    
    
    % analyse the pole data
    i=1;j=1;
    for phi=[0:inc:phimax phimax] %around z-axis
        j=1;
        for psi=[0:inc:psimax psimax] % out from z-axis
            
            z=cos(psi);
            x=sin(psi)*cos(phi);
            y=sin(psi)*sin(phi);
            test_normal=[x y z];
            
            angle_dev = acos((dot(poles, repmat(test_normal, size(poles,1), 1), 2)));
            density(i,j)=sum(weights(find(angle_dev<searchR)));
            j=j+1;
        end
        if mod(i, round(length([0:inc:phimax phimax])/10))==0
            fprintf('%d percent done\n', round(100*i/length([0:inc:phimax phimax])))
        end
        i=i+1;
    end
    
    % before normalising, save raw density values
    save_density=density;
    % normalise density into something meaningful
    % total solid angle covered 4pi for a whole sphere
    % we have reduced the search area
    search_solid_angle=pi*(searchR^2);
    % random density is
    random_density=valid_poles_weight*(search_solid_angle/total_solid_angle);
    density=density./random_density;
    
else
    % we need the range of phi and phi
    if strcmp(app.range, 'sst')
        % for standard stereographic triangle
        phimax=pi/4;
        psimax=atan(sqrt(2));
        total_solid_angle=pi/12;
    elseif strcmp(app.range, 'phi90')
        % for a 90 degree section of the pole figure
        phimax=pi/2;
        psimax=pi/2;
        total_solid_angle=pi/2;
    elseif strcmp(app.range, 'phi360')
        phimax=2*pi;
        psimax=pi/2;
        total_solid_angle=2*pi;
    else
        disp('unrecognised plot range :-(')
        return
    end
    
    % normalise density passed in using valid poles_weight
    search_solid_angle=pi*(searchR^2);
    random_density=valid_poles_weight*(search_solid_angle/total_solid_angle);
    density=density/random_density;
end

if app.plot_figure
    % build the figure
    % have: density as a function of phi and psi.
    % Can convert to x,y,z, and use surf.
    % modify for steorographic projection
    phi=[0:inc:phimax phimax]; %around z-axis
    psi=[0:inc:psimax psimax]; %away from z axis
    %radius distorted for stereographic proj
    r=tan(psi/2);
    [phi,r]=meshgrid(phi,r);
    [x,y]=pol2cart(phi,r);
    
    if app.new_figure
        hfig=figure();
    end
    
    surf(x,y,density')
    hold on
    shading('interp')
    axis equal
    
    if ~strcmp(app.mark_poles, 'none') %draw poles
        
        % Get symmetry operators
        symm = gtCrystGetSymmetryOperators([], spacegroup); 
        
        % cubic materials
        if spacegroup==225 || spacegroup==229
            
            if strcmp(app.mark_poles, 'main')
                uvw_poles=[0 0 1; 1 0 1; 1 1 1];
            elseif strcmp(app.mark_poles, 'bcc')
                uvw_poles=[0 0 1; 1 0 1; 1 1 1; 1 1 2; 2 1 3];
            elseif strcmp(app.mark_poles, 'all')
                uvw_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
                disp('unrecognised "mark_poles" option - should be all/main/none')
            end
            uvw_poles2=[];
            for i=1:length(symm)
                uvw_poles2=[uvw_poles2; uvw_poles*symm(i).g];
            end
            uvw_poles=[uvw_poles2; -uvw_poles2];
            
        % hexagonal materials
        elseif spacegroup==194 || spacegroup==167 || spacegroup==663
            uvw_poles_hex=[...
                1 0 -1 0;...
                0 0 0 1;...
                0 1 -1 2;...
                0 1 -1 1;...
                1 1 -2 3;...
                1 1 -2 0];
            
            % permute
            uvw_poles_hex2=[];
            for i = 1:size(uvw_poles_hex,1)
                uvw_poles_hex2 = [uvw_poles_hex2 ; gtGetReflections(uvw_poles_hex(i,:), spacegroup)];
            end
            uvw_poles_hex=uvw_poles_hex2;
            
            % convert to cartesian
            uvw_poles(:,1)= uvw_poles_hex(:,1) + 0.5 * uvw_poles_hex(:,2);
            uvw_poles(:,2)= 3^0.5/2 *uvw_poles_hex(:,2);
            uvw_poles(:,3)= uvw_poles_hex(:,4);
            % account for cell parameters
            uvw_poles(:,1)=uvw_poles(:,1)*2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1));
            uvw_poles(:,2)=uvw_poles(:,2)*2/(sqrt(3)*parameters.cryst(phaseid).latticepar(1));
            uvw_poles(:,3)=uvw_poles(:,3)/parameters.cryst(phaseid).latticepar(3);
            % add invert through origin
            uvw_poles=[uvw_poles; -uvw_poles];
        else
            disp('spacegroup not recognised...')
        end
        
        % get the relevent uvw_poles
        if strcmp(app.range, 'sst')
            %for standard stereographic triangle
            dum=find(uvw_poles(:,2)>=0 & uvw_poles(:,2)<=uvw_poles(:,1) & uvw_poles(:,3)>=uvw_poles(:,1));
            uvw_poles=uvw_poles(dum, :);
        elseif strcmp(app.range, 'phi90')
            %for a 90 degree section of the pole figure
            dum=find(uvw_poles(:,1)>=0 & uvw_poles(:,2)>=0 & uvw_poles(:,3)>=0);
            uvw_poles=uvw_poles(dum, :);
        elseif strcmp(app.range, 'phi360')
            %for 360 pole figure - reduce only using psi
            uvw_poles=uvw_poles(find(uvw_poles(:,3)>=0),:);
        end
        
        % plot these uvw/hkl poles
        phi_uvw=atan2(uvw_poles(:,2), uvw_poles(:,1));
        psi_uvw=acos(uvw_poles(:,3)./(sqrt(sum(uvw_poles.*uvw_poles, 2))));
        r_uvw=tan(psi_uvw/2);
        dummy=find(r_uvw>1);
        r_uvw(dummy)=[];
        phi_uvw(dummy)=[];
        [x_uvw,y_uvw]=pol2cart(phi_uvw, r_uvw);
        dummy=find(x_uvw<-1 | x_uvw>1 | y_uvw<-1 | y_uvw>1);
        x_uvw(dummy)=[];
        y_uvw(dummy)=[];
        
        plot3(x_uvw, y_uvw, ones(size(x_uvw))*(max(density(:))+1), 'k*')
        
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% plot rings on phi360 figures
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if strcmp(app.range, 'phi360')
        for ring=15:15:75
            ring_phi=(pi/180)*[0:3:360 0];
            ring_psi=(pi/180)*ring*ones(size(ring_phi));
            % convert to cartesian
            ring_r=tan(ring_psi/2);
            [ring_x,ring_y]=pol2cart(ring_phi,ring_r);
            ring_z=ones(size(ring_phi))*(max(density(:))+1);
            plot3(ring_x, ring_y, ring_z, 'w-');
            if 0 % label rings
                h=text(ring_x(1)+0.01, 0, (max(density(:))+1), num2str(ring));
                set(h, 'color', 'w');
            end
        end
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % tidy the figure
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    set(gca, 'XTickLabel','');
    set(gca, 'YTickLabel','');
    set(gca, 'GridLineStyle','none');
    set(gca, 'Ycolor', 'w');
    set(gca, 'Xcolor', 'w');
    colorbar();
    
    %change to x-y view
    set(gca, 'CameraPosition', [(min(x(:))+max(x(:)))/2 (min(y(:))+max(y(:)))/2 5*(max(density(:))+1)]);
    set(gca, 'CameraUpVector', [0 1 0]);
    
       
    if strcmp(app.range, 'sst')
        %for sst - apply a mask to the figure
        
        %line from 101 to 111 is z=x
        line_phi=[0:0.01:pi/4 pi/4];
        linex=cos(line_phi);
        linez=linex;
        liney=sin(line_phi);
        line_points=[linex' liney' linez'];
        %normalise
        line_points=line_points./repmat(sqrt(sum(line_points.*line_points, 2)), 1, 3);
        %convert to stereographic proj.
        line_psi=acos(line_points(:,3));
        line_r=tan(line_psi/2);
        [line_x,line_y]=pol2cart(line_phi', line_r);
        
        %white patch over unwanted bit of figure
        mask_x=[line_x ; sqrt(2)-1+0.045 ; sqrt(2)-1+0.045 ; line_x(1)];
        mask_y=[line_y ; line_y(end) ; 0 ; 0];
        mask_z=ones(size(mask_y))*(max(density(:))+0.5); %between the pole points and the surface
        p=patch(mask_x, mask_y, mask_z, 'w');
        set(p, 'edgecolor', 'w')
        
        %black outline on polefigure
        mask_x=[line_x ; 0 ; line_x(1)];
        mask_y=[line_y ; 0 ; 0];
        mask_z=ones(size(mask_y))*(max(density(:))+1);
        plot3(mask_x, mask_y, mask_z, 'k')
        
        %tweak axis limits and therefore adjust camera position
        set(gca, 'xlim', [-0.02 sqrt(2)-1+0.045])
        set(gca, 'ylim', [-0.025 0.3769])
        set(gca, 'CameraPosition', [(sqrt(2)-1+0.025)/2 (-0.025+0.3769)/2 (max(density(:))+2)]);
        set(gca, 'CameraUpVector', [0 1 0]);
        
        
        %redo the colorbar, and adjust position
        %warning('nice sst colorbar switched off')
        c=colorbar;
        set(c, 'Position', [0.18 0.3 0.04 0.6])
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % label poles
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if ~strcmp(app.label_poles, 'none')
        
        label_z=max(density(:))+1;
        
        if strcmp(app.label_poles, 'all') || strcmp(app.label_poles, 'main') || strcmp(app.label_poles, 'bcc')
            %label the poles son the corners of the sst
            text(0, 0-0.015, label_z, '(001)')
            text(0.366+0.01, 0.366-0.005, label_z, '(111)')
            text(0.4142-0.015, -0.015, label_z, '(101)')
        end
        
        if strcmp(app.label_poles, 'all')
            %label the long list of poles
            text(0.1213-0.04, 0.1213, label_z, '(114)')
            text(0.1583-0.04, 0.1583, label_z, '(113)')
            text(0.2247-0.04, 0.2247, label_z, '(112)')
            text(0.2808-0.04, 0.2808, label_z, '(223)')
            text(0.3052-0.045, 0.3052, label_z, '(334)')
            
            text(0.3845+0.01, 0.2884, label_z, '(434)')
            text(0.3901+0.01, 0.2601, label_z, '(323)')
            text(0.4000+0.01, 0.2000, label_z, '(212)')
            text(0.4077+0.01, 0.1359, label_z, '(313)')
            text(0.4105+0.01, 0.1026, label_z, '(414)')
            
            text(0.1231-0.02, -0.015, label_z, '(104)')
            text(0.1623-0.01, -0.015, label_z, '(103)')
            text(0.2361-0.015, -0.015, label_z, '(102)')
            text(0.3028-0.02, -0.015, label_z, '(203)')
            text(0.3333-0.01, -0.015, label_z, '(304)')
            
            text(0.2967+0.01, 0.1483, label_z, '(213)')
            text(0.2330+0.01, 0.1165, label_z, '(214)')
            text(0.3297+0.01, 0.1099, label_z, '(314)')
            text(0.3197+0.01, 0.2131, label_z, '(324)')
        elseif strcmp(app.label_poles, 'bcc')
            %label bcc slip planes
            text(0.2247-0.04, 0.2247, label_z, '(112)')
            text(0.3, 0.153, label_z, '(213)')
        end
    end
    
else
    hfig=[];
end

end % end of function
