function [bonds,status] = FMR_cause(bonds);

% Causality for Flow-modulated R component

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


status = -1;

% There are 2 ports; extract the information
e_1 = bonds(1,1);
f_1 = bonds(1,2);
e_2 = bonds(2,1);
f_2 = bonds(2,2);

% e_2 must be effort causality (imposes zero effort)
 if e_2 == 1 				% Conflict
   status = 1;
 else 					% Do the rest of the causality
   if e_2==0 				% Set to the fixed causality
     e_2 = -1;
   end;
   
   number_set = sum(sum([e_1 f_1 f_2]~=zeros(1,3) ));
   if number_set<2 			% Under causal
     status = -1;
   elseif number_set==2 		% Set the causality
     if f_2 == -1 			% Unicausal
       if e_1 == 0
 	e_1 = f_1;
       else
	 f_1 = e_1;
       end;
     elseif f_2 == 1 			% Bicausal
       e_1 = 1;
       f_1 = -1;
     elseif f_2 == 0
       if e_1==f_1
	 f_2 = -1;                      % Unicausal
       else
	 f_2 = 1;                      % Bicausal
       end;
     end;
     status = 0;
   elseif number_set==3 		% Check the causality
     if ( (f_2==-1)&(e_1~=f_1) )|( (f_2==1)&(e_1==f_1) )
       status = 1;
     else
       status = 0;
     end;
   end;
 end;
 
 bonds = [e_1 f_1
           e_2 f_2];

function [bonds,status] = juncause(bonds,jun,cause)
% [bonds,status] = juncause(bonds,jun,cause)

% Causality for  either effort or flow on  either zero or one junctions

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


j =  (3-cause)/2; % j is 1 for effort, 2 for flow

[n_bonds,junk] = size(bonds);
[causing_bond, n,other_bonds,m] = getindex(bonds(:,j),jun);
if n>1 % over causal
  status = 1;
elseif n==1 %causal
  status = 0;
  bonds(other_bonds(:,1),j) = -jun*ones(m,1);
else  % undercausal - try other way
  [causing_bond, n,other_bonds,m] = getindex(bonds(:,j),-jun);
  if n==n_bonds % over causal
    status = 1;
  elseif n==n_bonds-1 %causal
    status = 0;
    bonds(other_bonds(:,1),j) = jun*ones(m,1);
  else  % undercausal
    status = -1;
  end;
end;