%
%  alphafit.m
%
%  Gabbiani & Cox, Mathematics for Neuroscientists
%
%  Fit the transmitter model to an alpha function
%
%  usage:  ap = alphafit(tran,rec)
%
%  tran = struct('T0',1,'t1',0,'t2',4,'tfin',20,'dt',.1)  
%  rec = struct('kp',1.1,'km',.18,'G',100);
%
%   ap = [gmax tau_alpha]
%

function ap = alphafit(tran,rec)

[t,g] = gsyn(tran,rec);     % get the data to be fit

ap0 = [100 2];   		% initialize the guess

options = optimset('display','iter');

ap = lsqcurvefit(@afun,ap0,t,g,[],[],options);  % find best fit

hold on

plot(t,afun(ap,t),'r')		% plot best fit

hold off
box off
xlabel('t  (ms)','fontsize',14)
ylabel('g_{syn}  (mS/cm^2)','fontsize',14)
legend('transmitter model','alpha function')

return

function val = afun(ap,t)		% generic alpha function
tim = t/ap(2);
val = ap(1)*tim.*exp(1-tim);

%
% gsyn.m
%
% Gabbiani & Cox, Mathematics for Neuroscientists
%
%  gsyn(tran,rec)
%
%  tran = struct('T0',1,'t1',1,'t2',2,'tfin',10,'dt',.1) 
%  rec = struct('kp',1.1,'km',.19,'G',100);
%

function [t,g] = gsyn(tran,rec)

t = 0:tran.dt:tran.tfin;
R = tran.T0/(tran.T0+rec.km/rec.kp);
R = R*exp(rec.km*min(0,tran.t2-t));
R = R.*(1-exp((rec.kp*tran.T0+rec.km)*(tran.t1-min(t,tran.t2)))).*(t>tran.t1);
g = R*rec.G;
plot(t,g,'k')
return