%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Simulation of (leaky) integrate-and-fire neuron
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear all;
 close all;
 clf;

 % parameters of the model
 tau_inv=0.1;    % inverse time constant
 uu=0;           % initial membrane voltage
 I_ext=500*rand(100);       % constant external input
 %I_ext=50*sin(linspace(0,6*3.14,100));       % constant external input
 tspan=[0 100];  % integration interval
 theta=10;       % firing threshold

 %Integration with Euler method
 t_step=0; for it=1:100;
     t_step=t_step+1;
     x=uu<theta;
     uu=x*(1-tau_inv)*uu+tau_inv*I_ext(it);%+randn;
     u(t_step)=uu;
     t(t_step)=it;
     s(t_step)=1-x;
 end

 subplot('position',[0.13 0.13 1-0.26 0.6])
   plot(t,u);
   hold on; plot([0 100],[10 10],'--');
   axis([0 100 0. 12])
   xlabel('time [\tau]');
   ylabel('u(t)')

 subplot('position',[0.13 0.8 1-0.26 0.1])
   plot(t,s,'.','markersize',20);
   axis([0 100 0.5 1.5])
   set(gca,'xtick',[],'ytick',[])
   ylabel('spikes')


% Compute gain function
tref = [0 3 10];
theta = 10;
ures = 5;
tau = 10;
for k = 1:3
   for RI = 11:30
       r(RI,k) = (tref(k) - tau*log(((theta-RI)/(ures-RI))))^-1;
   end
end
figure,plot(r)