function [Xs,Ys,Zs,F1] = projectEIGS(vecs,sres,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PROJECTEIGS  projects unit vectors onto a unit sphere
% calculates the spherical KSDE of a population of unit vectors (place field eigenvectors or positional unit vectors)
% for every point on the unit sphere we calculate the spherical distance between every vector end point and the sphere
% point. These are then Gaussian weighted and summed. The final KDE is a map of these summed values.
%
% USAGE:
%         [Xs,Ys,Zs,F1] = projectEIGS(vecs,sres,sigma)
%
% INPUT:
%         vecs - unit vectors to be projected (Mx3)
%         sres - resolution (number of points) in spherical grid
%         sigma - the sigma/smoothing/SD of the Gaussian weighting, higher means a smoother output
%
% OUTPUT:
%    Xs,Ys,Zs - coordinates of the spherical grid, these are in a format that can be used directly by surf
%    F1 - the KDE map, values correspond to the coordinates in Sp
%
% EXAMPLES:
%
%     npoints = 1000; % number of random points
%     t = 2*pi*rand(npoints,1); % generate random azimuth
%     p = asin(-1+2*rand(npoints,1)); % generate random pitch
%     [X,Y,Z] = sph2cart(t,p,ones(npoints,1)); % convert to XYZ
%     [Xs,Ys,Zs,F1] = projectEIGS([X,Y,Z],64,10); % project these onto sphere
%     figure % plot the result
%     subplot(1,2,1)
%     plot3(X,Y,Z,'ko')
%     daspect([1 1 1])
%     axis xy tight
%     subplot(1,2,2)
%     surf(Xs,Ys,Zs,F1,'EdgeColor','none');
%     daspect([1 1 1])
%     axis xy tight
%
% See also: ACOSD DOT SPHERE

% HISTORY:
% version 1.0.0, Release 03/05/18 Initial release
% version 2.0.0, Release 14/11/18 Modified to use spherical gaussian
%
% Author: Roddy Grieves
% UCL, 26 Bedford Way
% eMail: r.grieves@ucl.ac.uk
% Copyright 2018 Roddy Grieves

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT ARGUMENTS CHECK
% deal with input variables
inps = {'sres','sigma'};
vals = {'64','10'};
reqd = [0 0];
for ff = 1:length(inps)
    if ~exist(inps{ff},'var')
        if reqd(ff)
            error('ERROR: vargin %s missing... exiting',inps{ff});            
        end        
        eval([inps{ff} '=' vals{ff} ';']);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FUNCTION BODY
% get sphere coordinates
[Xs,Ys,Zs] = sphere(sres);

F1 = NaN(size(Xs));
if isempty(vecs)
    return
end

sig1 = (sqrt(2*pi) .* sigma);
sig2 = 1/sig1;

vpoints = vecs;
vpoints(any(isnan(vpoints),2),:) = [];
F1 = NaN(size(Xs));
for i = 1:numel(Xs)
    vsphere = [Xs(i) Ys(i) Zs(i)];
    vsphere = repmat(vsphere,length(vpoints(:,1)),1);
    dp = acosd(dot(vsphere,vpoints,2));
    dp_norm = (exp(-0.5 * (dp./sigma).^2) ./ sig1) ./ sig2;
    F1(i) = sum(dp_norm);
end

% figure
% subplot(1,2,1)
% plot3(eigs(:,1),eigs(:,2),eigs(:,3),'ko')
% daspect([1 1 1])
% axis xy tight
% 
% subplot(1,2,2)
% surf(Xs,Ys,Zs,F1,'EdgeColor','none');
% daspect([1 1 1])
% axis xy tight