function [y,m] = IandF(in,simufreq,C,timeconst,initial_state)
%function [y,m] = IandF(in,simufreq,C,timeconst,initial_state)
%
% simulates an integrate and fire neuron
% in: the input current signal
% simufreq: the simulation frequency (should be a good deal higher than any
% input signal of interest)
% C: the neuron's capacitance
% timeconst: the neurons time constant (equivalent to R*C if you like)
% initial_state: the membrane voltage one timestep before the first
% m: the resulting membrane voltage
% y: the resulting spike train (spikes are 1 for one simulation timestep)

m=zeros(1, length(in));
y=zeros(1, length(in));
alpha=1-exp(-1/timeconst/simufreq);

m(1)=(1-alpha)*initial_state+in(1)/C/simufreq;
if m(1)>=1
    y(1)=1;
    m(1)=0;
end; %if
for i=2:length(in)
    m(i)=(1-alpha)*m(i-1)+in(i)/C/simufreq;
    if m(i)>=1
        y(i)=1;
       m(i)=0;
    end; %if
end; %for i
subplot(2,1,2);
plot((1:length(in))/simufreq,y);
subplot(2,1,1);
plot((1:length(in))/simufreq,m);