%% Copyright Radu Dogaru, 1995, 2006,
%  last update December 8, 2006
%
% e-mail: radu_d@ieee.org
%-------------------------------------------------------
% Appendix to Book:  
% R. Dogaru: "Universality and Emergent Computation in CNNs, World Scientific 2003", World Scientific, 2003
% http://www.worldscibooks.com/chaos/5053.html
%----------------------------------------------------------------------------------- 
% Returns a matrix ed1 of size N,N corresponding to 
% up to 13 different possible behaviors of the cell parameter points, 
% as obtained using the edge of chaos tests in equatiuons (19),(21), page 105 
% of the abovementioned book (one-diffusion coefficient, i.e. D1 non-zero). 
% The program can be easily extended to the case of "two difussion coefficients", i.e. D2 non zero.
% using (22) and (24) at page 108. 
% Here a FitzHugh-Nagumo model with parameters a,b,eps and c is considered 
% See eq. (3) pp. 99 in book 
% A parameter point (a,b) is associated with an (i,j) entry of the 
% ed1 matrix. The value of  ed1(i,j) can be as follows:
% 1 - cell is passive (P)
% 2 - cell is SA (stable and active) in at least one equilibrium point (EQ) 
%                                    typically there is only one EQ 
% 3- cell has 2 equilibria that are both SA, and one that is P - unlikely 
% 4- cell has 3 equilibria that are all SA - unlikely 
% 5- cell has only one equilibria and it is UA (Unstable ans Active)
% 6 - cell has an EQ which is P, an EQ which is SA and the last is UA 
% in the ablve, the cell can be considered stable and unstable, active and passive !
% 7 - cell has two EQ that are SA and one which is SU 
% 8 - unlikely
% 9 - cell has two EQ that are SU and one which is P
% 10 - cell has two EQ that are SU and one which is SA
% 11 
% 12
% 13 - cell has all three EQ unstable 
%---------------------------------------------------------------------------
% N is the number of sampples in the interval [amin,amax] and [bmin,bmax]

function [ed1]=edge_fhn(a,b,e,c,N)
% the first two parameters are given as ranges, 
% the others are fixed: 
% For instance  [ed1]=edge_fhn([-10 10],[-10 10],-0.1,1,100); vis_edge;
% interesting sub-regions within the (a,b) space may be further investigated 
% using the same program. The second function displays the section within the 
% (a,b) plane assigning certain colors (defined inside vis_edge.m) to each of the 
% above mentioned numbers (from 0 to 12) indicating various dynamic qualities 
% for the cell 


amin=a(1); amax=a(2);
bmin=b(1); bmax=b(2);

i=1:N; j=i;

a(j)=amin+(amax-amin)*(j-1)/(N-1);
b(i)=bmin+(bmax-bmin)*(i-1)/(N-1);


ed1=zeros(N,N);

for i=1:N
 for j=1:N
     
	  [x0]=roots([-b(i), 0, (3*c*b(i)-3), -3*a(j)]);  
     %find the equilibria of FitzHugh Nagumo cell
     ix=find(imag(x0)==0);
     x0=x0(ix);  
     % keep only the real valued solutions 
	  nsol=size(x0,1);
     % compute the elements of the Jacobian matrix that are independent 
     % of the number of root FitzHugh-Nagumo model 
     a12=-1;
     a21=-e;
     a22=e*b(i);
     %--------------------------------------------------------------
     % The scalars U, T and DEL
     U=[]; T=[]; DEL=[];
         
     % 1 diffusion case 
     typ=[];
     for k=1:nsol
         % compute the T,DEL and U coefficients required by the local activity tests
         % and the Jacobian parameter a11 which depends on the particular root 
      	a11(k)=c-x0(k)^2;
         T(k)=a11(k)+a22;
         U(k)=(a12-a21);
         DEL(k)=a11(k)*a22-a12*a21;
     		%-------------
           % test SA, P and SU conditions (19), (21) at p.105
           % and assign a number to the quality of the dynamic phenomena
           if ( (a22<0 & (T(k)<0 & T(k)>a22)) | (a22>0 & T(k)<0) ) & (DEL(k)>0) 
              typ(k)=1;  
           elseif (DEL(k)>0) & (T(k)<a22) & (T(k)<0) & (a22<0)
              typ(k)=0;  
           else
              typ(k)=4;
           end        	
      end
     % produce a unique number to identify the quality of the dynamic phenomena 
     % using a scalar  
     styp=sum(typ);
     ed1(i,j)=styp+1;
   end
end

%======================salveaza sectiunile==========================
disp('Save the result');
a=[amin amax];
b=[bmin bmax];
save sec0 ed1 a b e c
% display the result 
vis_edge;