%
% Gabbiani & Cox,  Mathematics for Neuroscientists
% 
% 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)
%  or
%    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 trapsyndrive

close all

gsyn = struct('gmax', 0.2, 'taua', 2, 't1', 5, 'Vsyn', 0);
[t,V1,g1] = trapsyn(0.01,35,gsyn);

gsyn = struct('gmax', 0.2*[1 1], 'taua', 2*[1 1], 't1', [4 5], 'Vsyn', [-68 0]);
[t,V2,g2] = trapsyn(0.01,35,gsyn);

figure(1)
plot(t,g1,'k')
hold on
plot(t,g2(:,1),'r',t,g2(:,2),'r--')
hold off
box off
legend('excitatory, solo','inhibitory, paired','excitatory, paired','location','best')
xlabel('t  (ms)','fontsize',14)
ylabel('g_{syn}  (\muS/cm^2)','fontsize',14)

figure(2)
plot(t,V1,'k')
hold on
plot(t,V2,'r')
hold off
box off
legend('excitatory only','inhibitory and excitatory','location','best')
xlabel('t  (ms)','fontsize',14)
ylabel('V  (mV)','fontsize',14)

function [t,V,g] = 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
