%
%  Gabbiani & Cox, Mathematics for Neuroscientists
%
% trapezoid on the fork with synapses
%
%  usage:   trapforksyn(stim,pinc)
%
% stim = struct('t1',1,'tau',1,'Gsyn',5e-7,'vsyn',70,'loc',100,'Tfin',20)
% pinc = 40
%
% stim = struct('t1',[1 3],'tau',[1/2 1/2],'Gsyn',1e-6*[1 1],'vsyn',...
%                    70*[1 1],'loc',[100 350],'Tfin',10)
%

function [t,vrec] = trapforksyn(stim,pinc)

a = 1e-4*[1 1 1];
ell = [2.5 2.5 2.5]/100;
dx = .0001;
dt = 0.05;
N = ell/dx;
A3 = 2*pi*a(3)*dx;
As = 4*pi*1e-6;
rho = A3/As;
R2 = 0.3; %0.034;
gL = 1/15; %0.3;
Cm = 1;
tau = Cm/gL;
lam = a/(2*R2*gL)/dx^2;   % lambda^2
r = a/a(3);
Hd = [2*lam(1)*ones(1,N(1)) 2*lam(2)*ones(1,N(2)) 2*lam(3)*ones(1,N(3)+1)];
Hd(1) = lam(1);
Hd(N(1)+1) = lam(2);
Hd(N(1)+N(2)+1) =  lam*r';
Hd(end) = rho*lam(3);
Hlen = length(Hd);

Hu = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)) -lam(3)*ones(1,N(3))];
Hl = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)-1) -r(2)*lam(2) -lam(3)*ones(1,N(3))];
Hl(end) = rho*Hl(end);

H = spdiags( [[Hl 0]' Hd' [0 Hu]'], -1:1, Hlen, Hlen);

H(N(1)+N(2)+1,N(1)) = -r(1)*lam(1);
H(N(1),N(1)+N(2)+1) = -lam(1);

I = speye(Hlen);

B = I + (I+H)*(dt/tau/2); 

eloc = round(Hlen*1e-4*stim.loc/sum(ell));
bloc = (eloc-1)*(Hlen+1) + 1;
dBe = B(bloc);

x3 = 0:dx:ell(3);
x1 = ell(3):dx:ell(3)+ell(1);
x2 = ell(3):dx:ell(3)+ell(2);

v = zeros(Hlen,1);         % initial conditions
rhs = zeros(Hlen,1);

t = 0;
tcnt = 0;
x = linspace(0,ell(3)+max(ell(1:2)),Hlen);
if pinc
   figure(1)
   v1 = fliplr(v(1:N(1))');
   v2 = fliplr(v(N(1)+1:N(1)+N(2))');
   v3 = fliplr(v(N(1)+N(2)+1:end)');
   plot3(x3,t*ones(size(x3)),v3,'r')
   hold on
   plot3(x2,t*ones(size(x2)),[v3(end) v2],'r')
   plot3(x1,t*ones(size(x1)),[v3(end) v1],'k')
end

Nt = ceil(stim.Tfin/dt);

vrec = zeros(length(eloc)+1,Nt);

Iapp = stim.Gsyn*(dt/2)/A3/Cm;
if eloc == Hlen
   Iapp = Iapp*A3/As;
end

c0 = Iapp.*(t./stim.tau).*exp(1-t./stim.tau).*(t>stim.t1);
t = dt;
c1 = Iapp.*(t./stim.tau).*exp(1-t./stim.tau).*(t>stim.t1);

r = zeros(Hlen,1);
r(eloc) = (stim.vsyn.*(c0 + c1))';

for j=2:Nt,

    B(bloc) = dBe + c1;

    v = B\r;

    vrec(:,j) = v([eloc Hlen]);

    if mod(j,pinc) == 0
         v1 = fliplr(v(1:N(1))');
         v2 = fliplr(v(N(1)+1:N(1)+N(2))');
         v3 = fliplr(v(N(1)+N(2)+1:end)');
         plot3(x3,t*ones(size(x3)),v3,'r')
         plot3(x2,t*ones(size(x2)),[v3(end) v2],'r')
         plot3(x1,t*ones(size(x1)),[v3(end) v1],'k')
    end

    t = j*dt;
    c0 = c1;
    c1 = Iapp.*((t-stim.t1)./stim.tau).*exp(1-(t-stim.t1)./stim.tau).*(t>stim.t1);
    r = 2*v - r;
    r(eloc) = r(eloc) + (stim.vsyn.*(c0 + c1))';

end

t = linspace(0,stim.Tfin,Nt);

if pinc
   xlabel('x (cm)','fontsize',14)
   ylabel('t (ms)','fontsize',14)
   zlabel('v (mV)','fontsize',14)
   axis tight
   hold off

   figure(2)
   plot(t,vrec(1,:),'k',t,vrec(2,:),'r')
   xlabel('t (ms)','fontsize',14)
   ylabel('(mV)','fontsize',14)
   box off
   
end
