%
% rcbe.m
%
% solve the rc single cell via backward euler
%
% GE' = -GE/sigE + sum SE_j sum delta(t-T_{j,k}) + FE sum delta(t-TF_k})
%                  Lj=E      k                         k
%
% GI' = -GI/sigI + sum SI_j sum delta(t-T_{j,k}) + FI sum delta(t-TF_k})
%                  Lj=I      k                         k
%
% V' = -GL*(V-epsL) - GE*(V-epsE) - GI*(V-epsI)
%
%  usage   rcbe(NE,NI,dt)
%
%    NE = number of excitatory cells
%    NI = number of inhibitory cells
%    dt = time step
%
%  example   rcbe(4,2,.1)

function rcbe(NE,NI,dt)

GL = 0.05;
epsL = 0;
epsT = 1;
epsR = 0;
tauref = 4;
sigE = 2;
sigI = 7;
epsE = 14/3;
epsI = -2/3;
SE = .15;
SI = .25;

rand('state',0);

TE = 1000*rand(NE);
TEsort = sort(reshape(TE,NE^2,1));

TI = 1000*rand(NI);
TIsort = sort(reshape(TI,NI^2,1));

Tfin = max(TEsort(end),TIsort(end));
tspike = -Tfin;

t = 0;
tcnt = 1;
tim(tcnt) = t;
spcnt = 0;

GE = zeros(ceil(Tfin/dt),1);
GI = GE;
V = GE;

while t < Tfin,

      t = t + dt;
      tcnt = tcnt + 1;
      tim(tcnt) = t;

      GE(tcnt) = SE*sum( exp( (TEsort(1:max(find(TEsort<t)))-t)/sigE ) );
   
      GI(tcnt) = SI*sum( exp( (TIsort(1:max(find(TIsort<t)))-t)/sigI ) );

      V(tcnt) = epsR;

      if t - tspike > tauref,

         num  = V(tcnt-1) + dt*(GL*epsL + GE(tcnt)*epsE + GI(tcnt)*epsI);
         den = 1 + dt*(GL + GE(tcnt) + GI(tcnt));
         V(tcnt) = num/den;

      end

      if V(tcnt) > epsT,

         tspike = t;
         spcnt = spcnt + 1;
         sptim(spcnt) = t;
         Vsp(spcnt) = V(tcnt);
         V(tcnt) = epsR;

      end

end

figure(1)
plot(tim(1:tcnt),GE(1:tcnt),'r')
hold on
plot(tim(1:tcnt),GI(1:tcnt))
hold off
xlabel('t  (ms)','fontsize',16)
ylabel('GE (red) and GI (blue)','fontsize',16)
figure(2)
plot(tim(1:tcnt),V(1:tcnt))
hold on
plot(sptim,Vsp,'r*')
hold off
xlabel('t  (ms)','fontsize',16)
ylabel('V  (1=thresh, 0=reset)','fontsize',16)

return