%% plot data:
x = load('../data/WT_01.mat');
wtdata = x.data;
plotcurrents(wtdata.t, wtdata.I);

x = load('../data/A1622D_01.mat');
addata = x.data;
plotcurrents(addata.t, addata.I);

%% I-V curve:
[wtsteps, wtpeaks] = ivcurve(wtdata.steps, wtdata.t, wtdata.I, 100.0);

[adsteps, adpeaks] = ivcurve(addata.steps, addata.t, addata.I, 100.0);

figure();
plot(wtsteps, wtpeaks, '-b');
hold on;
plot(adsteps, adpeaks, '-r');
hold off;

%% reversal potential:
wtE = reversalpotential(wtsteps, wtpeaks);
adE = reversalpotential(adsteps, adpeaks);

%% activation curve:
wtg = wtpeaks./(wtsteps - wtE);
adg = adpeaks./(adsteps - adE);

wtinfty = wtg(wtsteps<40.0)/mean(wtg((wtsteps>=20.0)&(wtsteps<=40.0)));
adinfty = adg(adsteps<40.0)/mean(adg((adsteps>=20.0)&(adsteps<=40.0)));
wtsteps = wtsteps(wtsteps<40.0);
adsteps = adsteps(adsteps<40.0);

figure();
plot(wtsteps, wtinfty, '-b');
hold on;
plot(adsteps, adinfty, '-r');

%% boltzmann fit:
bf = @(p, v) 1.0./(1.0+exp(-p(1)*(v - p(2))));
p = lsqcurvefit(bf, [1.0, -40.0], wtsteps, wtinfty);
wtfit = bf(p, wtsteps);
p = lsqcurvefit(bf, [1.0, -40.0], adsteps, adinfty);
adfit = bf(p, adsteps);

plot(wtsteps, wtfit, '-b');
plot(wtsteps, adfit, '-r');
hold off;