%%%%%%%%%%% BEGIN enc_code.m %%%%%%%%%%%%%%%%%
%%% this is the code called by the user edited gui set up in 'encounter.m'

Ns = str2num(get(h_Ns,'string')); 
Nt = str2num(get(h_Nt,'string')); 
speedNs = str2num(get(h_sNs,'string')); 
speedNt = str2num(get(h_sNt,'string')); 
search_radius = str2num(get(h_sr,'string')); 

a=figure
whitebg('k')
set(a,'Name', 'Encounter Rate Simulation')
set(a,'MenuBar','none')
set(a,'Color', 'k')
set(a,'NumberTitle','off')

%% space variables
domainX = 1000;      % size of the twoD arena
domainY= domainX;

%% time variables
dT=.1; % time increments
T_max=100; % simulation duration

%% position and movement variables
Theta_s = 2 * pi .*(0.5-rand(Ns,1)); % initial angle of movement searchers
Theta_t = 2 * pi .*(0.5-rand(Nt,1)); % initial angle of movement targets

Xpos_Nt=rand(Nt,1)*domainX;
Ypos_Nt=rand(Nt,1)*domainY;

Xpos_Ns=rand(Ns,1)*domainX;
Ypos_Ns=rand(Ns,1)*domainY;

%% place holder for encounter stats
encounters=[];

% time loop advancing organisms and scoring their encounter
for t=1:dT:T_max;
    
    Theta_s = 2 * pi .*(0.5-rand(Ns,1));
    Theta_t = 2 * pi .*(0.5-rand(Nt,1));
    
    us = cos(Theta_s).*speedNs;	%determine the movement in the X plane for the searcher
    vs = sin(Theta_s).*speedNs;	%determine the movement in the Y plane
    
    Xpos_Ns= Xpos_Ns + us *dT;
    Ypos_Ns= Ypos_Ns + vs *dT;
    
    Xpos_Ns = mod(Xpos_Ns,domainX);
    Ypos_Ns = mod(Ypos_Ns,domainY);
    
    ut = cos(Theta_t).*speedNt;	%determine the movement in the X plane for the target
    vt = sin(Theta_t).*speedNt;	%determine the movement in the Y plane
    
    Xpos_Nt= Xpos_Nt + ut .*dT;
    Ypos_Nt= Ypos_Nt + vt .*dT;
    
    Xpos_Nt = mod(Xpos_Nt,domainX);
    Ypos_Nt = mod(Ypos_Nt,domainY);
    
    
    %% quantifying encounters by tallying the number of targets within a
    %% searchers sensory radius 
    hits_sum=0;
    hits_list =[];
    
    for p=1:Ns
        hits = find((Xpos_Ns(p)-Xpos_Nt).^2+(Ypos_Ns(p)-Ypos_Nt).^2 < search_radius^2);
        hits_sum=hits_sum+length(hits);
        hits_list =[hits_list;hits];
    end
    
    if mod(t,1)==0    % results are plotted
        plot(Xpos_Ns, Ypos_Ns, 'ro')
        hold on
        plot(Xpos_Nt, Ypos_Nt, 'g*')
        plot(Xpos_Nt(unique(hits_list)), Ypos_Nt(unique(hits_list)), 'w*')
        xlabel('X domain')
        ylabel('Y domain')
        title(['Time =    ',num2str(t)]);
        axis([0,domainX,0,domainY])
        hold off
        drawnow
        shg
        encounters =[encounters;hits_sum];
    end
end

average_encounter_rate=round(mean(encounters)/Ns)

text(domainX/3,domainY/2,['Average encounter rate =   ',num2str(round(average_encounter_rate))])

pause(3)

%%% encounter rate is plotted over time
if sum(encounters) >0
    figure
    plot(1:T_max,encounters/Ns,'b*')
    hold on
    plot(1:T_max, mean(encounters)/Ns,'b-')
    %axis([0,T_max,0,max(encounters)+max(encounters)*.01])
    axis([0,T_max,0,max(encounters)/Ns])
    xlabel('Time')
    ylabel('Encounters')
end

%%%%%%%%%%% END enc_code.m %%%%%%%%%%%%%%%%%