function BondOrder=bond_order_parameter(center_positions, N, l)
%Bond_order_prm=bond_order_parameter(center_positions,N)
%
%Calculates the bond order parameter given the center positions of
%particles. This quantitity should be unity for a perfect hexagonal
%lattice. Bond_order_prm uses spherical harmonics Ylm(theta,phi) in the
%computation and the bond order of each particle is given for positive
%values of m - m=0:l.
%Ref: http://www.seas.harvard.edu/weitzlab/bondorder.html
%Q_{lm}(i)=1/N * \sum_j=1^N Ylm(\hat{Rij}) where N=number of neighbors of
%the ith particle.
%For bugs, comments, etc. please contact Kapil Krishan :
%kkrishan@uci.edu
%www.physics.uci.edu/~foams
%version: I

%% Full 3D version, but optimized for 2D
BondOrder=zeros(size(center_positions,1),length(0:l));

for ic=1:length(center_positions)
    cp0=center_positions(ic,:);
    cpN=center_positions(N(ic,:)==1,:);
    
    xyzN=cpN-ones(size(cpN,1),1)*cp0;
    
    %2D optimization below:
    if size(xyzN,2)==2
        xyzN(:,3)=zeros(size(xyzN,1),1);
    end
    
    [th,phi]=cart2sph(xyzN(:,1),xyzN(:,2),xyzN(:,3));
    
    YlmSUM=0;
    for icc=1:size(th,1)
         Ylm=spherical_harmonics(l,phi(icc),th(icc));
         YlmSUM=YlmSUM+Ylm;
    end
    BondOrder(ic,:)=YlmSUM/size(th,1);
end

%% Earlier code based on input from EW that I dont get.
% phi_i=zeros(size(N,1),1);
% for ic=1:size(N,1)
%     nbrs=find(N(ic,:)==1);
%     v_ij=center_positions(nbrs,[1 2])-ones(size(nbrs,2),1)*center_positions(ic,[1 2]);
%     alpha_ij=angle(v_ij(:,1)+sqrt(-1)*v_ij(:,2));
%     phi_i(ic)=sum(exp(6*sqrt(-1)*alpha_ij))/size(alpha_ij,1);
% end
% 
% Bond_order_prm=(abs(sum(phi_i))/size(phi_i,1))^2;