
Ncycles=50;
Nmembers=200;
Ndiseases=4;
killed_by_deasese=zeros(Ncycles,Ndiseases);

member_stat=zeros(1,Nmembers);
% let's make them all live
member_stat=member_stat+4;

member_stat_map=zeros(Ncycles,Nmembers);

% let's disease a few unlucky individuals
Ndiseased2=20;
for i=1:Ndiseased2
	m=ceil(Nmembers*rand());
	member_stat(m)=2; % hard to heal but weakly contagious
end
Ndiseased3=20;
for i=1:Ndiseased3
	m=ceil(Nmembers*rand()); % easy to heal but contagious
	member_stat(m)=3;
end

member_stat_map(1,:) = member_stat; % first day situation recorded

%plot the initial colony status
figure(1);
xlim([1, Nmembers]);
bar(member_stat);
figure(2);
bar(killed_by_deasese);
drawnow('expose');
%pause; %wait for user

for c=2:Ncycles % on cycle one we just initialize the colony
	if c~=1
		killed_by_deasese(c,:)=killed_by_deasese(c-1,:);
	end
	% spread diseases
	for i=1:Nmembers
		disease = member_stat(i);

		[self_action, ...
		disease_left_neighbor, ...
		disease_right_neighbor] = disease_outcome(disease);
		
		if (i-1 >= 1)
			% we have left neighbor
			if ( disease_left_neighbor == true)
				if(member_stat(i-1) ~= 1)
					% only alive guys can catch a disease
					member_stat(i-1)=disease;
				end
			end
		end

		if (i+1 <= Nmembers)
			% we have right neighbor
			if ( disease_right_neighbor == true)
				if(member_stat(i+1) ~= 1)
					% only alive guys can catch a disease
					member_stat(i+1)=disease;
				end
			end
		end

		if ( (self_action == 1) && (disease ~=1) ) % we should not count already dead
			%calculate death toll
			killed_by_deasese(c,disease)=killed_by_deasese(c,disease)+1;
		end
		member_stat(i)=self_action;
	end

	% update member stat vs day map
	member_stat_map(c,:) = member_stat;

	figure(1); % only real illness count 
	bar(killed_by_deasese(:,[2,3]),'stacked'); % notice array slicing
	legend('death by hard to heal but not contagious', 'death by easy to heal but highly contagious');	
	xlabel('cycle number');
	title('death toll');
end

figure(1); % call figure 2 to attention
print('-depsc', 'death_toll.eps');

% plot the evolution of the colony
figure(2);
colony_colorscheme=[ ...  
	% color coded as RGB triplet
	0,0,0; % black is for dead
	0,0,1; % blue is for weakly contagious
	1,0,0; % red is for highly contagious
	0,1,0; % green is for healthy
	];
image( member_stat_map ); 
colormap(colony_colorscheme);
xlabel('member position');
ylabel('cycle');
title('Colony evolution');





		
