function [err, FId, FCd] = SAkin9proofM(Gmc,Gse,rs,N,k)
% [err, FId, FCd] = SAkin9proofM(Gmc,Gse,rs,N,k)
% 9 state kinetic trap model for 2x2 proofreading tiles
%   forward rate f   = k*exp(-Gmc)       k = 20*10^6
%   reverse rate r_b = k*exp(-b*Gse)
%   freezing rate rs for complete 2x2 blocks
%   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     f       2f      f      4Nf  V    f     2f       f       f     rs
%   FI <-- M4 <--> M3 <--> M2 <--> M1 <--> E <--> C1 <--> C2 <--> C3 <--> C4 --> FC
%             2r2      r1     2r2      r1     r2      r2      2r2     r2 

% compute steady-state values for each intermediate

%      E           M1      M2      M3     M4       C1     C2     C3     C4
M = [ -(4*N+1)*f   r1      0       0      0        r2     0      0      0; ...
        4*N*f   -(f+r1)   2*r2     0      0        0      0      0      0; ...
        0          f   -2*(f+r2)   r1     0        0      0      0      0; ...
        0          0      2*f   -(f+r1)  2*r2      0      0      0      0; ...
        0          0       0       f  -(2*r2+rs)   0      0      0      0; ...
        f          0       0       0      0    -(2*f+r2)  r2     0      0; ...
        0          0       0       0      0       2*f  -(f+r2)  2*r2    0; ...
        0          0       0       0      0        0      f  -(f+2*r2)  r2; ...
        0          0       0       0      0        0      0      f   -(r2+rs) ];

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

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