% 
% Solve the Passive Membrane Equation subject to a
% synaptic input, via the trapezoid rule
%
%  usage   trapsyn(dt,Tfin,gsyn)
%
%    gsyn = struct('gmax', 0.02, 'taua', 2, 't1', 1, 'Vsyn', 0)
%    gsyn = struct('gmax', 0.02*[1 1], 'taua', 2*[1 1], 't1', [1 3], 'Vsyn', [-68 0])
%
%  e.g.    trapsyn(0.01,40,gsyn)
%

function [t,V] = trapsyn(dt,Tfin,gsyn)

VCl = -68;	    % mV
Cm = 1;         % micro F/cm^2
gCl = 0.3;      % mS/cm^2
z = 2/dt;

Nt = ceil(1+Tfin/dt);                    % number of time steps
V = zeros(Nt,1);  t = V;  	 % preallocate space
g = zeros(Nt,length(gsyn.gmax));

j = 1;
V(1) = VCl;
a0 = gCl/Cm;
b0 = gCl*VCl/Cm;

for j=2:Nt,

   t(j) = (j-1)*dt;

   g(j,:) = gsyn.gmax.*((t(j)-gsyn.t1)./gsyn.taua).*exp(1-(t(j)-gsyn.t1)./gsyn.taua).*(t(j)>gsyn.t1);
   a1 = (gCl+sum(g(j,:)))/Cm;
   tmp = g(j,:)*(gsyn.Vsyn)';
   b1 = (gCl*VCl + g(j,:)*(gsyn.Vsyn)')/Cm;
   
   V(j) = ( (z-a0)*V(j-1) + b0+b1 ) / (z+a1);
   
   a0 = a1;
   b0 = b1;

end

if 1==0
[ax,h1,h2] = plotyy(t,1000*g(:,1),t,V)
hold on
plot(ax(1),t,1000*g(:,2),'r--')
hold off
box off
set(ax(2),'ylim',[-68 -65])
set(ax(2),'ytick',[-68 -67 -66 -65])
xlabel('t  (ms)','fontsize',14)
ylabel(ax(1),'g_{syn}  (\muS/cm^2)','fontsize',14)
ylabel(ax(2),'V  (mV)','fontsize',14)
end