%% general settings:
tau = 10.0;
dt = 0.1;

%% test passive membrane:
tmax = 50.0;
time = [0:dt:tmax];
E = zeros(length(time), 1);
V = passivemembrane(1.0, E, tau, dt);
expfun = exp(-time/tau);
figure()
plot(time, V, 'b', 'linewidth', 2);
hold on;
plot(time, expfun, 'r');
hold off;

%% step input:
tmax = 120.0;
time = [0:dt:tmax];
E = zeros(length(time), 1);
E(20/dt:70/dt) = 10;
V = passivemembrane(0.0, E, tau, dt);
figure()
plot(time, E, 'r');
hold on;
plot(time, V, 'b', 'linewidth', 2);
hold off;

%% sine waves:
tmax = 1000.0;
time = [0:dt:tmax];
freqs = [1.0 10.0 30.0 100.0];
figure();
for k = 1:length(freqs)
    f = freqs(k);
    E = sin(2*pi*0.001*f*time);
    V = passivemembrane(0.0, E, tau, dt);
    subplot(4, 1, k);
    plot(time, E, 'r');
    hold on;
    plot(time, V, 'b', 'linewidth', 2);
    hold off;
end


%% transfer function:
tmax = 1000.0;
time = [0:dt:tmax];
taus = [3.0 10.0 30.0];
figure();
for tau = taus
    freqs = [1.0:1.0:100.0];
    rates = zeros(length(freqs), 1);
    for k = 1:length(freqs)
        f = freqs(k);
        E = sin(2*pi*0.001*f*time);
        V = passivemembrane(0.0, E, tau, dt);
        rates(k) = max(V(100/dt:end));
    end
    plot(freqs, rates);
    hold on;
end
hold off;


%% leaky IaF:
tau = 10.0;
tmax = 50.0;
time = [0:dt:tmax];
E = zeros(length(time), 1) + 1.5;
[spikes, V] = lifspikes(0.0, E, tau, dt);
figure()
plot(time, V, 'b', 'linewidth', 2);
hold on;
plot(spikes, ones(length(spikes), 1), 'or');
hold off;


%% leaky IaF and sine input:
tau = 10.0;
tmax = 500.0;
time = [0:dt:tmax];
f = 10.0;
E = 2.0*sin(2*pi*0.001*f*time);
[spikes, V] = lifspikes(0.0, E, tau, dt);
figure()
subplot(2, 1, 1);
plot(time, V, 'b', 'linewidth', 2);
hold on;
plot(spikes, ones(length(spikes), 1), 'or');
hold off;


%% transfer function of LIF spikes:
tmax = 1000.0;
time = [0:dt:tmax];
taus = [3.0 10.0];
figure();
for tau = taus
    freqs = [5.0:0.1:30.0];
    rates = zeros(length(freqs), 1);
    for k = 1:length(freqs)
        f = freqs(k);
        E = 2.0*sin(2*pi*0.001*f*time);
        [spikes, V] = lifspikes(0.0, E, tau, dt);
        spikes = spikes(spikes>100.0);
        if length(spikes) > 2
            rates(k) = 1000.0*(length(spikes)-1)/(spikes(end)-spikes(1));
        else
            rates(k) = 0.0;
        end
    end
    plot(freqs, rates);
    hold on;
end
hold off;