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).

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,n) = V(j-1,n) + eta * dV;


        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
   end
end



