function [S,V] = iafsim(I,R,C,Vt,Tref,tstep)
% [S,V] = iafsim(I,R,C,Vt,Tref,tstep)
% This is a routine to implement the integrate and fire model of a
% neuron. The inputs to the routine is the current injected into
% the neuron. The output of the routine is the spike train vector
% output of the neuron in response to the current.
%
% ***********  Description of Parameters and Variables ************
%
% V            = membrane voltage waveform (V)
% S            = spike waveform (0-1)
% Vt           = voltage threshold     (V)
% I            = input current vector  (A)
% R            = membrane resistance   (Ohm)
% C            = membrane capacitance  (F)
% Tref         = refractory period     (sec)
% tstep        = (optional) time step  (sec)
%
% The default timestep (eta) for the routine is 1 msec
% which can be overridden by the user by supplying tstep (in sec).

%Vthr = 1.0;

if (nargin < 6)
  eta = 0.001;
else
  eta = tstep;
end

L          = size(I,1);
trials     = size(I,2);
S          = zeros(L,trials);
V          = zeros(L,trials);

for n = 1:trials
   Tlast = -Tref;
   for j=2:L
     if(j*eta > Tlast + Tref)      % If the time past last spike > Tref

        %***************************************************
        %*                  EULER UPDATE RULE
        %***************************************************
	dV = -V(j-1,n)/R/C + I(j-1,n)/C;

        % FILL IN THE NEXT STEP OF THE EULER INTEGRATION:
        % V(j,n) = ?


        if(V(j,n) >= Vt)              % If the threshold voltage is exceeded
           V(j,n) = 0;                  % Set the voltage to Vrest for t <Tref
           Tlast  = j*eta;              % Record the spike time
           S(j,n) = 1;
        end
     end
     j=j +1;                           % Increment time counter
   end
end