% LIF.m
clear

Ntrials = 51;


E=-75e-3;
Rm=10e6;
tau_m=10e-3;
Vth=-45e-3;
Vreset=-80e-3;
Iappmax = 5e-9;

sigma = 1e-6;

dt=0.0002;
tmax=20;
tpulse=0.0;
lengthpulse=tmax-tpulse;
Nt=tmax/dt;
T=0:dt:tmax;

for trial = 1:Ntrials

    Iapp(trial) = Iappmax*(trial-1)/Ntrials

    I=zeros(size(T));
    V=zeros(size(T));

    noise = zeros(size(T));
    for i = 1:length(T)
        noise(i) = randn(1)*sigma*sqrt(dt);
    end

    V(1)=E;
    tref=0.002;
    lastspike=-2*tref;
    spikes=zeros(size(T));
    for i = tpulse/dt+1 : (tpulse+lengthpulse)/dt
        I(i) = Iapp(trial);
    end
    for i = 2:Nt;
        V_inf = E+Rm*(I(i-1)+noise(i-1));
        V(i) = V_inf + (V(i-1) - V_inf)*exp(-dt/tau_m);

        if ( T(i) < lastspike + tref )
            V(i) = Vreset;
        end
        if V(i) > Vth
            V(i) = Vreset;
            spikes(i) = 1;
            lastspike=T(i);
        end
    end

    rate(trial) = sum(spikes)/lengthpulse

end

plot(Iapp,rate)