function f = firefly_synch_neighbor(N,NN,NR,AnimFlag)

%syntax: 
%f = firefly_synch(N,NN,AnimFlag);
%N - the number of fireflies
%NN - the number of (nearest) neighbors to look at for synchronization
%NR - the number of random neighbors to look at for synchronization
%AnimFlag - animation flag, if it's 1 create animation, otherwise don't

%Dr. Zoran Nenadic, UC Irvine
%02/01/2007

close all

%parameters of the simulation time vector
T = 20;     %duration
dt = 0.01;  %integration step
Period = 0.9; %the period between flashes


%set the flashing frequencies of individual fireflies 
Omega = 2*pi./(Period*[1 + 0.1*rand(1,N)]);

%set the duration of flashing
Frac = 0.25;
Duration = Frac * Period;

%initialize the phases of flashing (theta = 2*j*pi triggers flashing)
Theta0 = [0 rand(1,N-1)]'*2*pi;

%create the position of the fireflies
x = [rand(1,N)-0.5];
y = [rand(1,N)-0.5];

disp('calculating distances ...')
Dist = squareform(pdist([x' y']));
%exclude the distance to itself by making it very large
Dist(Dist == 0) = 10^6;
[temp J] = sort(Dist,2);
clear temp

if NR ~= 0 
    ind_nr = ceil(rand(N,NR) * [N-1-NN]) + NN;
else 
    ind_nr = [];
end

if NN ~= 0 
    ind_nn = 1:NN;
else
    ind_nn = [];
end

if isempty(ind_nr)
    Jn = J(:,ind_nn);
else
    for i = 1:N
        Jn(i,:) = J(i,[ind_nn ind_nr(i,:)]);
    end
end

disp('solving ...')
OPT = odeset('RelTol',1e-6,'AbsTol',[1e-6]);
[t,Theta] = ode45(@phase_dyn,[0:dt:T],Theta0,OPT,Omega,Jn); 

%flashing indicator sequence (1-flash, 0-do not flash)
Seq = zeros(N,length(t)); 

for i = 1:N
    
    %find samples of Theta which are negative
    ind_neg = find(sin(Theta(:,i)) <= 0);
    
    %form the difference vector
    dind = diff(ind_neg);
    
    %number that's not 1 indicates that zero crossing with positive
    %slope--> must be a flashing point
    flash_ind = ind_neg(dind~=1);
    
    for k = 1:length(flash_ind)
        
        %extend the duration beyond the flashing point
        End_pt = min([flash_ind(k)+round(Duration/dt), length(t)]);
        Seq(i,flash_ind(k):End_pt) = 1;
    end
end


SumSeq = sum(Seq,1);

%heuristic measure of synchrony
SumSeqSort = fliplr(sort(SumSeq));
Sync = mean(SumSeqSort(1:round(Frac*T/dt)))/[N];




if AnimFlag == 1
    
    plot(t,SumSeq/N)
    xlabel('Time')
    ylabel('Population Average Behavior')

    fh = figure;
    set(fh,'Position',[100 50 420 420])
    th1 = text(-0.8,-0.8,['Time: ' num2str(0) ' [s]'],'Color','y', ...
        'FontWeight','b');
    th2 = text(0.4,-0.8,['Sync.: ' num2str(Sync*100,'%2.2f') ' [%]'],'Color','y', ...
        'FontWeight','b');
    set(gcf,'DoubleBuffer','on','ToolBar','none','MenuBar','none')
    hold on
pause
    for i = 1:N
        ph(i) = plot(x(i),y(i),'ro','MarkerFaceColor','r','MarkerSize',4, ...
            'Visible','off');
    end
    set(gca,'Xlim',[-1 1],'YLim',[-1 1],'Color','k', ...
        'XTickLabel',[],'YTickLabel',[],'Position',[0 0 1 1], ...
        'DataAspectRatio',[1 1 1])


    for k = 1:length(t)

        set(th1,'String',['Time: ' num2str(dt*(k-1),'%2.2f') ' [s]'])
        ind = find(Seq(:,k)==1);
        ind_c = setdiff(1:N,ind);
        set(ph(ind),'Visible','on')
        set(ph(ind_c),'Visible','off')
        %drawnow
        pause(0.01)

    end
else
    %do nothing
end

f = Sync;

function dTheta = phase_dyn(t,Theta,Omega,J)

A = 1.0;
N = length(Theta);

dTheta = zeros(N,1);

for i = 1:N
    dTheta(i) = Omega(i) + A*ones(1,size(J,2))*sin(Theta(J(i,:))-Theta(i));
end