function [err, FId, FCd] = SAkin3origM(Gmc,Gse,rs,N,k)
% [err, FId, FCd] = SAkin3origM(Gmc,Gse,rs,N,k)
% 3 state kinetic trap model for original BCA tiles
%   forward rate f   = k*exp(-Gmc)       k = 20*10^6
%   reverse rate r_b = k*exp(-b*Gse)
%   freezing rate rs for each tile addition
%   assumes a BCA tileset with N+1 symbols
%
%   uses matrix inversion to do the calculation, rather than algebra.

f  = k*exp(-Gmc);
r2 = k*exp(-2*Gse);
r1 = k*exp(-Gse);

%                 | 1
%       rs   2Nf  V   f    rs
%   FI <-- M <--> E <--> C --> FC
%             r1     r2  


%       E           M        C
M = [ -(2*N+1)*f    r1       r2    ; ...
      2*N*f      -(r1+rs)    0     ; ...
        f           0    -(r2+rs)   ];

xss = inv(M) * [-1 0 0]';   % steady-state, dumping 1 unit/sec into E

FId = rs*xss(2);
FCd = rs*xss(3);
err = FId/(FId+FCd);     % this should be exactly equal to FId