clear;
global T0;
Vmax = 21;
Vmin = -9;
Tplot = 6;

DT1 = 0.2; %0.2
DT2 = 0.4; %0.3
DT3 = 0.7; %0.4
v0 = 1; dt = 0.001; dv = 7; vr = 0;

FPT = @(v) pi/2 - atan(v);
Vsol = @(t,v0) tan(t + atan(v0));
Phi = @(v) (atan(v)-atan(vr)) / FPT(vr);
Pot = @(phi) tan(phi*FPT(vr) + atan(vr));
Delta = @(phi, dv) (atan(Pot(phi)) - atan(Pot(phi)-dv)) / FPT(vr);

T0 = FPT(vr);
AL = 7; % ArrowLength
phi0 = Phi(-1e8);

fs = 14;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(413);
set(gca,'fontsize',fs);

spike1 = FPT(v0);
Time = 0:dt:spike1;
Volt = Vsol(Time,v0);
plot(Time,Volt,'k-','linewidth',2); hold on;
plot([spike1 spike1], [vr Vmax], 'k-','linewidth',2);

T = [0:dt:T0] + Time(end);  % Unperturbed spike
Volt = Vsol([0:dt:T0],vr);
plot(T,Volt,'k-'); 
plot([T(end) T(end)], [vr Vmax], 'k-');

Time = [0:dt:DT1] + Time(end);
Volt = Vsol([0:dt:DT1],vr);
plot(Time,Volt,'k-','linewidth',2); 
plot([Time(end) Time(end)], [min(Volt) 5], 'k:');

plot([Time(end) Time(end)], [Volt(end) Volt(end)-dv], 'k-','linewidth',2);
Time = [0:dt:DT2] + Time(end);
Volt = Vsol([0:dt:DT2], Volt(end)-dv);
plot(Time,Volt,'k-','linewidth',2); 
%plot([Time(end) Time(end)], [min(Volt) 5], 'k:');

plot([Time(end) Time(end)], [Volt(end) Volt(end)-dv], 'k-','linewidth',2);
Time = [0:dt:DT3] + Time(end);
Volt = Vsol([0:dt:DT3], Volt(end)-dv);
plot(Time,Volt,'k-','linewidth',2); 
%plot([Time(end) Time(end)], [min(Volt) 5], 'k:');

spike2 = FPT(Volt(end)-dv);
Actual = (spike2 + DT1 + DT2 + DT3) / T0
plot([Time(end) Time(end)], [Volt(end) Volt(end)-dv], 'k-','linewidth',2);
Time = [0:dt:spike2] + Time(end);
Volt = Vsol([0:dt:spike2], Volt(end)-dv);
plot(Time,Volt,'k-','linewidth',2);

Tail = Tplot - spike1 - spike2;
plot([Time(end) Time(end)], [Vmax vr], 'k-','linewidth',2);
Time = [0:dt:Tail] + Time(end);
Volt = Vsol([0:dt:Tail],vr);
plot(Time,Volt,'k-','linewidth',2); 

phi1 = DT1/T0;
Delta1 = Delta(phi1, dv);

phi2 = phi1 - Delta1 + DT2/T0;
Delta2 = Delta(phi2, dv);

phi3 = phi2 - Delta2 + DT3/T0;
Delta3 = Delta(phi3, dv);

Predicted = Delta1 + Delta2 + Delta3

t0 = spike1; ph0=1-t0/T0;
t1 = t0+DT1;
t2 = t0+DT1+DT2;
t3 = t0+DT1+DT2+DT3;
t4 = t0+DT1+DT2+DT3+spike2;
t5 = t0+DT1+DT2+DT3+spike2+Tail;

tt = [0, t0, t0, t1, t1, t2, t2, t3, t3, t4, t4, t5];
phase = [ph0, 1, 0, phi1, phi1-Delta1, phi2, phi2-Delta2, phi3, phi3-Delta3, 1, 0, Tail/T0]; 

tt0 = [t0, t0+T0, t0+T0];
pp0 = [0, 1, 0]; 

%arrow([t0 10], [t0+T0 10], 'Ends', 'both', 'Length', 9);
text(t0 + T0/2 - 0.2, 13, 'T_0');
%arrow([t0 15], [t4 15], 'Ends', 'both', 'Length', 12);
text((t0 + t4)/2 - 0.2,  18, 'T_3');

ylabel('V(t)');
plot([0 Tplot], [0 0], 'k:');
axis([0 Tplot Vmin Vmax] );
set(gca,'box','off');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(414);
set(gca,'fontsize',fs);

pfs = 13;
plot(tt,phase,'k-','linewidth',2); hold on;
plot(tt0,pp0,'k:','linewidth',1);
hold on;
axis([0 Tplot min(phase)-0.1 1.2]);
plot([0 Tplot], [1 1], 'k:');
plot([0 Tplot], [0 0], 'k:');

plot(t1, phi1, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);
plot(t2, phi2, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);
plot(t3, phi3, 'ok', 'markersize', 5, 'markerfacecolor', [0 0 0]);

h=text(t1+0.1, phi1+0.1, '\phi_1', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_1$');
h=text(t2+0.1, phi2+0.1, '\phi_2', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_2$');
h=text(t3+0.1, phi3+0.1, '\phi_3', 'fontsize', pfs); 
%set(h,'Interpreter','latex','string','$\varphi_3$');

plot([t1 t1],[min(phase) 0.75],':k');
plot([t2 t2],[min(phase) 0.75],':k');
plot([t3 t3],[min(phase) 0.75],':k');

arrow([t1 0.7], [t2 0.7], 'Ends', 'both', 'Length', AL);
arrow([t2 0.7], [t3 0.7], 'Ends', 'both', 'Length', AL);
text(0.5*(t1 + t2)-0.1, 0.85, '\Deltat_{1}', 'fontsize', fs);
text(0.5*(t2 + t3)-0.1, 0.85, '\Deltat_{2}', 'fontsize', fs);

xlabel('time, ms'); 
h = ylabel('\phi(t)','fontsize',pfs);
set(gca,'box','off');
%set(h,'Interpreter','latex','string','$\varphi(t)$');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(413);
arrow([t0 10], [t0+T0 10], 'Ends', 'both', 'Length', AL);
arrow([t0 16], [t4 16], 'Ends', 'both', 'Length', AL);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,1,[1 2]);
set(gca,'fontsize',fs);

phi = phi0:0.01:0;
plot(phi, Delta(phi,dv), 'b-', 'linewidth', 3 ); hold on;

phi = 0:0.01:0.99;
plot(phi, Delta(phi,dv), 'k-', 'linewidth', 1 );

plot([0 0], [0 max(Delta(phi,dv))], 'color', [0.6 0.6 0.6]);

axis([phi0 1 0 max(Delta(phi,dv))]);

xlabel('\phi');
ylabel('\Delta(\phi)');
set(gca,'box','off');

plot([phi1 phi1], [0 Delta(phi1,dv)], 'k:');
plot([phi0 phi1], [Delta(phi1,dv) Delta(phi1,dv)], 'k:');
plot(phi1, Delta(phi1, dv), 'o', 'markerfacecolor', [0 0 0]);
text(phi1-0.03, Delta(phi1, dv)+0.05, '\phi_1', 'fontsize', fs);

plot([phi2 phi2], [0 Delta(phi2,dv)], 'k:');
plot([phi0 phi2], [Delta(phi2,dv) Delta(phi2,dv)], 'k:');
plot(phi2, Delta(phi2, dv), 'o', 'markerfacecolor', [0 0 0]);
text(phi2-0.03, Delta(phi2, dv)+0.05, '\phi_2', 'fontsize', fs);

plot([phi3 phi3], [0 Delta(phi3,dv)], 'k:');
plot([phi0 phi3], [Delta(phi3,dv) Delta(phi3,dv)], 'k:');
plot(phi3, Delta(phi3, dv), 'o', 'markerfacecolor', [0 0 0]);
text(phi3-0.03, Delta(phi3, dv)+0.05, '\phi_3', 'fontsize', fs);

