%function out=pulse_loop(Conc,tmax,kf,kr,n_pulse)
%
% solve ode's for a looped chemical reaction
% The loop equations are singular because the reactants sum to a constant (1)
% We reduce the n x n system to an n-1 x n-1 system with the transformation
% MR = M * P  where P = [eye(n-1) zeros(n-1,1); -1 -1 -1 ... 1];
% MR = MR(1:n-1, 1:n-1)    beta = MR(1:n-1,n)
%
Conc = 1.e-5;
tmax = 250;
kf = 1;
kr = .0001;
n_pulse = 10;

dt = 1/16;
n = 4;				% rank of the ode system
Y = zeros(n);
y = zeros(n-1);
M = zeros(n,n);
rate = zeros(4);

Vrest = -80;
Vtest = -20;

%Conc = 1e-5 ;	% forward velocity:  Conc * k = .004 /msec

%					Rates are in units of /msec
%
Y = [1.0 0 0 0]';
y = Y(1:n-1);
P = [eye(n-1) zeros(n-1,1); -1 -1 -1 1];

rest(1) = Y(1);
inactive(1) = Y(2);
in_block(1) = Y(3);
rest_block(1) = Y(4);

V(1) = Vtest;
V(2) = Vrest;
T(1) = 50;
T(2) = tmax -T(1);

for p=1:n_pulse
   for i=1:2
       M = gate(V(i),Conc,kf,kr);
       MR = M*P;
       beta = MR(1:n-1,n);
       MR = MR(1:n-1,1:n-1);
%      Minv = inv(eye(n-1) - dt*MR);
       for k=1:T(i)
           for j=1:1/dt
%             y = Minv*(y + dt*beta);
              y = (eye(n-1) - dt*MR)\(y + dt*beta);
           end
          rest(k+T(1)*(i-1)+1 + (p-1)*(T(1)+T(2))) = y(1);
    	   inactive(k+T(1)*(i-1)+1 + (p-1)*(T(1)+T(2))) = y(2);
    	   in_block(k+T(1)*(i-1)+1 + (p-1)*(T(1)+T(2))) = y(3);
    	   rest_block(k+T(1)*(i-1)+1 + (p-1)*(T(1)+T(2))) = 1 - sum(y(1:3));
 		end
    end
end
total_block = in_block + rest_block;

t = linspace(0,tmax*n_pulse,n_pulse*tmax+1);
plot(t,rest,'r-',t,inactive,'b-',t,in_block,'g-',t,rest_block,'m-',t,total_block,'c'); 
%plot(t,rest,'r-',t,inactive,'b-',t,in_block,'g-'); 
axis([0 tmax*n_pulse 0 1]);
xlabel('Time (msec)');
ylabel('Fraction of inactivated channels');

%out = [t' rest' inactive' in_block' rest_block' total_block'];