function out_status=simau1

%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

global run_status id_step id_protoc 
global h_main
global cmap

%protocole 1=nb_step 2=p_stim  4=nb_stim 5= nb_disc
global dim_protoc protocole liste_stim liste_disc



global dim_x dim_y
global etat  %actual state in etat(dim_x,dim_y)
 

%state_space
global nb_state_ex nb_state_arp

%neighbourhoo
global threshold mat_vois wei_vois

%stepping
global stepping
if(isempty(stepping))stepping=0;end;

global etat  %actual state in etat(dim_x,dim_y)

%trace 
global liste_tr pile_tr id_tr

if(strcmp(run_status,'stop'))break;end;
% Initialize
% protocole and step


out_status=0;
encore=1;
init_prot=1;

disp(run_status);

	if(strcmp(run_status,'start'))
	
	run_status='continue';
	id_protoc=1;
	id_step=0
		if(~isempty(liste_tr))
		id_tr=0;
		pile_tr=[];
                ii=find(liste_tr(1,:) <= dim_x & liste_tr(2,:) <= dim_y);
                    if(length(ii) < length(liste_tr(1,:)))
                    a=liste_tr(:,ii);
                    liste_tr=[];liste_tr=a;clear a; 
                    end
                          
		end;
	end;
disp('begin');
% neighbour matrix

	nb_neigh=length(mat_vois(:,1));
	for i=1:nb_neigh

	mat_neigh(i)=-mat_vois(i,1)-dim_x*mat_vois(i,2);
	end;
	mat_neighx=-mat_vois(:,1);
	mat_neighy=-mat_vois(:,2);

%non excitable state and excitable state;
	ii=find(threshold == 0);
	non_excs=ii-1;

	ii=find(threshold >  0);
	excs=ii-1;

%matrix to trace the patch
clear cx cy
for j=1:dim_y
for i=1:dim_x
%cx(j+(i-1)*dim_y,:)=[i-1 i i i-1];
%cy(j+(i-1)*dim_y,:)=[j-1 j-1 j j];
cx(i+(j-1)*dim_x,:)=[i-1 i i i-1];
cy(i+(j-1)*dim_x,:)=[j-1 j-1 j j];
end;
end

%color and first draw and identification

figure=h_main
ha=findobj(h_main,'type','axes');

   if(~isempty(ha))
      cla;
      
   else
     ha=axes('units', 'norm', 'pos', [0, 0, 1 1],'vis','off');
   end

cmap=hot(nb_state_ex+nb_state_arp+1);
colormap(cmap);
vect_colmap=colormap;

set(ha, 'xlim', [0 dim_x], 'ylim', [ 0 dim_y]);
caxis([1 nb_state_ex+nb_state_arp+1]);


for i=1:dim_x*dim_y
pile_patch(i)=patch(cx(i,:),cy(i,:),etat(i)+1, 'eras', 'none','UserData',i);
end;
drawnow







while(strcmp(run_status,'continue') & encore )


% interrupt

	if(strcmp(run_status,'stop'))
	out_status=2;  % interrupt
	return;
	end;

%if change protocole
	if(id_step==protocole(1,id_protoc))

		% normal termination
		if(id_protoc == length(protocole(1,:)))
		out_status=1; 
		run_status='stop';
		return;

		else
		% prepare for next protocole
		id_protoc=id_protoc+1;
		id_step=0;
		init_prot=1;
		end
	end;

%if begin protocole

if(init_prot)
	nb_step=protocole(1,id_protoc);		
	p_stim=protocole(2,id_protoc);
	nb_stim=protocole(3,id_protoc);
	amp_stim=protocole(5,id_protoc);

		if(nb_stim ==0)amp_stim=0;end;
		if(amp_stim ==0)nb_stim=0;end;
		if(nb_stim ==0 & amp_stim==0)
		p_stim=nb_step+1;
		else

		clear id_stimule
		ii=find(liste_stim(3,:) == id_protoc);	

			for i=1:length(ii)		
%			id_stimule(i)=liste_stim(2,ii(i))+...
%                	dim_y*(liste_stim(1,ii(i))-1);
			id_stimule(i)=liste_stim(1,ii(i))+...
                	dim_x*(liste_stim(2,ii(i))-1);

			end

		end





	nb_disc=protocole(4,id_protoc);

	if(nb_disc > 0)
		clear id_disc
		ii=find(liste_disc(3,:) == id_protoc);

			for i=1:length(ii)		
%			id_disc(i)=liste_disc(2,ii(i))+...
%                	dim_y*(liste_disc(1,ii(i))-1);
			id_disc(i)=liste_disc(1,ii(i))+...
                	dim_x*(liste_disc(2,ii(i))-1);



			end

	end

end; % init_prot
% make one step

passe=etat;




% treat disconnected cell first such that they cannot
% be stimulated nor stimulate


	if(nb_disc > 0)
	passe(id_disc)=-1*ones(length(id_disc),1);
	end

% treat all non-excitable states first;

for id_state=1:length(non_excs)

	ii=find(passe == non_excs(id_state));
	if(length(ii) > 0)etat(ii)=passe(ii)+1;end
end


% treat all excitable states

%step without stimulation

if((rem(id_step+1,p_stim) ~= 1 & p_stim ~= 1) | amp_stim==0)

for id_state=1:length(excs)


	id_r=find(passe == excs(id_state));


	for i=1:length(id_r)
%find neighbourhood



	ii=find(ceil(id_r(i)/dim_x)+mat_neighy > 0 & ...
		   ceil(id_r(i)/dim_x)+mat_neighy <= dim_y &...
                   rem(id_r(i)-1,dim_x)+1+mat_neighx > 0 & ...
                   rem(id_r(i)-1,dim_x)+1+mat_neighx <= dim_x);


	act_neigh=id_r(i)+mat_neigh(ii);
        poids=wei_vois(ii);



%	if( length(...
%	find(passe(act_neigh) > 0 & passe(act_neigh) <= nb_state_ex))...
%	>=threshold(excs(id_state)+1) )

	ii=find(passe(act_neigh) > 0 & passe(act_neigh)<= nb_state_ex);
	
        if(sum(poids(ii)) >=threshold(excs(id_state)+1))

	etat(id_r(i))=1;
	elseif(passe(id_r(i)) > 0)
 	etat(id_r(i))=passe(id_r(i))+1;
	end

	end;
		
end;% exc state



%step with stimulation
else
for id_state=1:length(excs)


	id_r=find(passe(:) == excs(id_state));


	for i=1:length(id_r)
%find neighbourhood




	ii=find(ceil(id_r(i)/dim_x)+mat_neighy > 0 & ...
		   ceil(id_r(i)/dim_x)+mat_neighy <= dim_y &...
                   rem(id_r(i)-1,dim_x)+1+mat_neighx > 0 & ...
                   rem(id_r(i)-1,dim_x)+1+mat_neighx <= dim_x);
	act_neigh=id_r(i)+mat_neigh(ii);

        poids=wei_vois(ii);

%	if( length(...
%        find(passe(act_neigh) > 0 & passe(act_neigh) <= nb_state_ex))+...
%	(~isempty(find(id_stimule == id_r(i))))*amp_stim... 
%	>=threshold(excs(id_state)+1))

	ii=find(passe(act_neigh) > 0 & passe(act_neigh)<= nb_state_ex);
	
        if(sum(poids(ii))+(~isempty(find(id_stimule == id_r(i))))*amp_stim >=threshold(excs(id_state)+1))

	etat(id_r(i))=1;
	elseif(passe(id_r(i)) > 0)
 	etat(id_r(i))=passe(id_r(i))+1;
	end

	end;
		
end;% exc state

end; % with and without stimulation

% retreat disconnected cell
% they can evolve but cannot be stimulated

	if(nb_disc > 0)
	passe(id_disc)=etat(id_disc);
		if(~isempty(find(passe(id_disc) > 0)))

	 	ii=find(passe(id_disc) > 0);
	 	etat(id_disc(ii))=passe(id_disc(ii))+1;
		end
	end	

% return all state over total state to 0

	ii=find(etat > (nb_state_ex+ nb_state_arp));
	if(length(ii) > 0)etat(ii)=zeros(length(ii),1);end;


%affiche nouvel etat
ii=find(passe(:) ~= etat(:));
for i=1:length(ii)
   set(pile_patch(ii(i)),'cdata', etat(ii(i))+1,'facecolor',...
        vect_colmap(etat(ii(i))+1,:),'edgecolor',[0 0 0]);
end




id_step=id_step+1

%if trace on
	if(~isempty(liste_tr))
	id_tr=id_tr+1;
	for i=1:length(liste_tr(1,:))
	pile_tr(id_tr,i)=etat(liste_tr(1,i),liste_tr(2,i));
	end;
	end;

if(stepping ==1)
k = waitforbuttonpress;
end;
drawnow
end;


	
