
Ncycles=50;
Nmembers=200;
Ndiseases=4;
killed_by_deases=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 cvontageous
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);
bar(member_stat);
figure(2);
bar(killed_by_deases);
drawnow('expose');
%pause; %wait for user

for c=2:Ncycles % on cycle one we just initialize the colony
	if c~=1
		killed_by_deases(c,:)=killed_by_deases(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_deases(c,disease)=killed_by_deases(c,disease)+1;
		end
		member_stat(i)=self_action;
	end

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

	%% plot the current colony members status
	% we show the status of the members in the bar chart 
	% with the same hight for every member 
	% but different color corresponding to a disease
	figure(1);
	hold off;
	x=(1:Nmembers);
	%dummy plot to reserve enough space in bar chart
	bar(x, (x*0+1));
	x_dead=x(member_stat == 1);
	if (length(x_dead) >=1 )
		hold on;
		bar( x_dead, 1+0*member_stat(x_dead), 'k');
	end
	x_weakly_contagious=x(member_stat == 2);
	if (length(x_weakly_contagious) >=1 )
		hold on;
		bar( x_weakly_contagious, 1+0*member_stat(x_weakly_contagious), 'r');
	end
	x_highly_contagious=x(member_stat == 3);
	if (length(x_highly_contagious) >=1 )
		hold on;
		bar( x_highly_contagious, 1+0*member_stat(x_highly_contagious), 'b');
	end
	x_alive=x(member_stat == 4);
	if (length(x_alive) >=1 )
		hold on;
		bar( x_alive, 1+0*member_stat(x_alive), 'g');
	end
	title('health level among members at the given day');
	hold off;


	figure(2); % only real illness count 
	bar(killed_by_deases(:,[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

% plot the evolution of the colony
figure(3);
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');




		
