function [xx, yy, zz] = gtCalculateDensityPhiPsi(poles, weights, phirange, psirange, searchR, N)
% Calculates the x,y,z coordinates for a density plot in the given phi/psi range
% all the angles are in radians

phimin = phirange(1);
phimax = phirange(2);
psimin = psirange(1);
psimax = psirange(2);

% calculate solid angle
total_solid_angle = (phimax - phimin)*(cos(psimin) - cos(psimax));
% increment in phi and psi
inc = pi/N;

ii=1;
for phi=[phimin:inc:phimax phimax] %around z-axis
    jj=1;
    for psi=[psimin:inc:psimax psimax] % out from z-axis

        tmp_z=cos(psi);
        tmp_x=sin(psi)*cos(phi);
        tmp_y=sin(psi)*sin(phi);
        test_normal=[tmp_x tmp_y tmp_z];
        test_normal = repmat(test_normal, size(poles,1), 1);

        angle_dev = acos(dot(poles, test_normal, 2));
        zz(ii,jj)=sum(weights(angle_dev<searchR));

        jj=jj+1;
    end
    ii=ii+1;
end
zz = zz';

phi=[phimin:inc:phimax phimax];
psi=[psimin:inc:psimax psimax];
% stereographic projection on plane
rr=tan(psi/2);
[phi,rr]=meshgrid(phi,rr);
[xx,yy]=pol2cart(phi,rr);

valid_poles_weight=sum(weights);
% 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
normalization=valid_poles_weight*(search_solid_angle/total_solid_angle);
zz = zz./normalization;


end