function  [t, v, spk]  = liandfr(t_stop, t_step, i_inj)
%
% [t, v, spk] = liandfr(t_stop, t_step, i_inj)
%
%liandfr: returns the membrane potential and spike train vector 
%for a leaky integrate-and-fire neuron receiving 
%a constant current i_inj as input (units: nA).

%variable initialization
t = 0:t_step:t_stop;
n_t = length(t);
v = zeros(1,n_t);
spk = zeros(1,n_t);

%LIF parameters
vr = -70; %in mV 
vth = -54; %in mV
tref = 2; %in msec
tau = 30; %in msec
rin = 20; %in MOhms

vth0 = vth - vr; %threshold with respect to rest
tref0 = round(tref/t_step); %refractory period in units of the time step

%save time by computing in advance 2 integration
%constants for forward Euler
f1 = 1 - t_step/tau;
f2 = (rin*t_step)/tau;

%Constant for backward Euler
b1 = 1/(1+t_step/tau);

v(1) = 0; %initial condition
ref_time = 0; %no refractory period at start

for i = 2:n_t
    if ( ref_time <= 0 )
        %refractory period over, integrate membrane
        %potential with forward or backward Euler
        
        %v(i) = f1*v(i-1) + f2*i_inj; %forward Euler
        v(i) = b1*(v(i-1) + f2*i_inj); %backward Euler
        
        if (v(i) > vth0)
            %time to spike
            spk(i) = 1;
            v(i) = 0; %reset membrane potential 
            ref_time = tref0; %arm refractory period counter
        end;
    else
      %refractory period not over yet
      %decrement refractory period counter
      %wait until refractory period is over
      ref_time = ref_time - 1;
    end;
end;