function [X1,X2,kap]=steady_states(diff_orthants,circuitsBperp,circuitsM,Bperp,M)
xx1=diff_orthants;
X1=[];
X2=[];
Ma=zeros(12,6);
Ma(1,1)=1;
Ma(2,1)=1;
Ma(4,2)=1;
Ma(5,2)=1;
Ma(7,3)=1;
Ma(8,3)=1;
Ma(10,4)=1;
Ma(11,4)=1;
Ma(1,5)=1;
Ma(3,5)=1;
Ma(4,5)=1;
Ma(6,5)=1;
Ma(12,6)=1;
Ma(10,6)=1;
Ma(9,6)=1;
Ma(7,6)=1;
kap=[];
for k=1:size(xx1,1)
  V2=circuitsBperp(xx1{k,2},:); %circuits of Bperp contained in the orthant
  for i=1:size(V2,1)
    zz=sign(V2(i,:)).*xx1{k,1};
    ww=find(zz);
    if zz(ww(1))==-1
      V2(i,:)=-V2(i,:); %we make the nonzero signs of the circuit and the orthant coincide
    endif
  endfor %i
  W2=circuitsM(xx1{k,3},:);  %circuits of M contained in the orthant
  for i=1:size(W2,1)
    zz=sign(W2(i,:)).*xx1{k,1};
    ww=find(zz);
    if zz(ww(1))==-1
      W2(i,:)=-W2(i,:);%we make the nonzero signs of the circuit and the orthant coincide
    endif
  endfor %i
  a=zeros(1,10); %this will be a vector in the intersection of the orthogonal of the "binomial" subspace and the orthant
  s=zeros(1,10); %this will be a vector in the intersection of the stoichiometric subspace and the orthant
  for i=1:size(V2,1)
    a=a+1*V2(i,:); %we obtain "a" as the sum of the circuits of Bperp contained in the orthant
  end
  for i=1:size(W2,1)
   s=s+1*W2(i,:);%we obtain "s" as the sum of the circuits of M contained in the orthant
  end
  %X1, first steady state:
  for i=1:10
    if a(i)~=0
      X1(k,i)=s(i)/(exp(a(i))-1);
    else 
      X1(k,i)=1; %we set X1_i=1 if a_i=s_i=0;
    end
  end
  D=diag(exp(a));
 %X2, second steady state:
  X2(k,:)=(D*X1(k,:)')';
 % reaction constants:
  M1=diag([X1(k,1)*X1(k,9),X1(k,5),X1(k,5),X1(k,2)*X1(k,10),X1(k,6),X1(k,6),X1(k,2)*X1(k,3),X1(k,7),X1(k,7),X1(k,4)*X1(k,10),X1(k,8),X1(k,8)]);
  lam=[1,1,1,1,1,1];%we can set this vector as any positive vector
  kap=[kap;(inv(M1)*Ma*lam')'];
end
