%2014-01-23

% Passive membrane equation

clear all
close all

C = 1;
El = -65;
Gl = 0.1;
Iapp = 2;
D = 0;

Vth = -50;
Vrst = -65;


T = 1000;
dt = 0.1;
t = 0:dt:T;

ti = dt;
tf = 10000;
H = zeros(1,length(t));
H(floor(ti/dt):floor(tf/dt))=1;

V = zeros(1,length(t));
V(1) = El;

for j=1:length(t)-1
    kv1 = (-Gl*(V(j)-El)+Iapp*H(j))/C;
    av = V(j)+kv1*dt;
    kv2 = (-Gl*(av-El)+Iapp*H(j+1))/C;
    V(j+1) = V(j) + (kv1+kv2)*dt/2;
    if V(j+1)>=Vth
       V(j)=60;
       V(j+1)=Vrst;
   end
end


tspk = zeros(1);
ISI = zeros(1);
Vn = zeros(1,length(t));
Vn(1) = V(1);

spkcnt = 0;
for j=1:length(t)-1
    eta = randn;
    Vnaux = Vn(j)+sqrt(2*D*dt)*eta;
    kv1 = (-Gl*(Vnaux-El)+Iapp*H(j))/C;
    av = Vn(j)+kv1*dt;
    av = av + sqrt(2*D*dt)*eta;
    kv2 = (-Gl*(av-El)+Iapp*H(j+1))/C;
    Vn(j+1) = Vn(j) + (kv1+kv2)*dt/2;
    Vn(j+1) = Vn(j+1) + sqrt(2*D*dt)*eta;
    if Vn(j+1)>=Vth
       Vn(j)=60;
       Vn(j+1)=Vrst;
       spkcnt = spkcnt+1;
       tspk(spkcnt)=t(j);
   end
end

for k=1:length(tspk)-1;
    ISI(k) = tspk(k+1)-tspk(k);
end

figure
hold on
%plot(t,V,'r','linewidth',2);
plot(t,Vn,'b','linewidth',2);
axis([0 T -80 80]);
set(gca,'fontsize',20);
xlabel('t')
ylabel('V');

if D>0
    figure
    hist(ISI,20);
    set(gca,'fontsize',20);
    xlabel('ISI')
    title('Histogram');
end


savefile='IntegrateFireData.dat'
Data = [t; Vn]';
save(savefile,'Data','-ascii');

savefile='Isi.dat'
Data = ISI';
save(savefile,'Data','-ascii');

