function saccade_oculomotor_plant_fig
%SACCADE_OCULOMOTOR_PLANT_FIG
%   Visualise the oculomotor plant as a second-order linear system:
%   - Show how the two time constants (T1, T2) shape the impulse response
%   - Show the plant's response to different inputs:
%       1) pulse
%       2) step
%       3) pulse–step
%
%   Intended as optional illustration code for the Neurobiophysics course.

%% General settings
close all;  % close any existing figures
close all;              % close old figures
clearvars;              % clear workspace (local to this function)
clc;                  % optional: clear command window

% Check whether my custom functions are available
hasCbrewer  = ~isempty(which('cbrewer'));
hasNicegraph = ~isempty(which('nicegraph'));
hasSavegraph = ~isempty(which('savegraph'));
hasBftext = ~isempty(which('bf_text'));

% if hasCbrewer
%     % Use custom cbrewer colormap if available
% 	col = cbrewer('qual','Dark2',3);
% else
%     % Fallback to built-in colormap
%     col = lines(3);
% end

%%
dt   = 0.001;          % time step [s]
tmax = 0.8;            % total time [s]
t    = 0:dt:tmax;      % time vector

% Default time constants of the oculomotor plant (in seconds)
T1 = 0.150;            % long time constant (~elasticity)
T2 = 0.012;            % short time constant (~viscosity)

%% 1. Impulse response and effect of T1, T2
% The plant is modeled as a second-order overdamped linear system.
% Its impulse response h(t) is:
%   h(t) = 1/(T1 - T2) * [exp(-t/T1) - exp(-t/T2)]

hconv		= impulse_response(t, T1, T2);

% For illustration: vary T1 and T2 to see their effects
h_T1		= impulse_response(t, T1, 0);   % remove T2
h_T2		= impulse_response(t, 0, T2);   % remove T1

%% 2. Define three input commands: pulse, step, pulse–step
t_on		= 0.10;		% movement onset time [s]
pulse_dur	= 0.050;	% pulse duration [s]

% Step input: tonic command from t_on onwards
u_step		= double(t >= t_on);

% Pulse input: brief burst starting at t_on
A_pulse		= 2.5;           % amplitude of pulse (arbitrary units)
u_pulse		= A_pulse * double(t >= t_on & t < t_on + pulse_dur);

% Pulse–step input: strong pulse on top of a tonic step
A_step		= 1;            % amplitude of step
u_pulsestep = u_pulse + A_step * u_step;

%% 3. Compute plant outputs by convolution (linear system)
% For a linear time-invariant system:
%   y(t) = (u * h)(t) = convolution of input u with impulse response h

U		= dt * conv(u_step,      hconv, 'full'); U = U(1:numel(u_step));
h		= dt * conv(u_pulse,     hconv, 'full'); h = h(1:numel(u_step));
hU		= dt * conv(u_pulsestep, hconv, 'full'); hU = hU(1:numel(u_step));

%% 4. Plot results

% time in msec
t		= t*1000;
t_on	= t_on*1000;
tmax	= tmax*1000;

figure(1);
clf;

% (1) Impulse responses: effect of T1 and T2
subplot(1,3,1);
plot(t, hconv,'-','LineWidth', 2); hold on;
plot(t, h_T1,'-','LineWidth', 1.5);
plot(t, h_T2,'-','LineWidth', 1.5);
xlabel('time (ms)');
ylabel('h(t)');
title('impulse response');
xlim([0 tmax]-t_on)
set(gca,'YTick',[]);
ylim([-1.1 1.1]);

legend( ...
    sprintf('T_1=%.0f ms, T_2=%.0f ms', 1e3*T1, 1e3*T2), ...
    sprintf('only T_1'), ...
    sprintf('only T_2'), ...
    'Location', 'southeast','FontSize',14);

% (2) Input commands: pulse, step, pulse–step
subplot(1,3,2);
plot(t-t_on,u_pulse,'LineWidth', 2); hold on;
plot(t-t_on,u_step+3,'LineWidth', 2);
plot(t-t_on,u_pulsestep+5,'LineWidth', 2);
xlabel('time re input (s)');
ylabel('command (a.u.)');

set(gca,'YTick',[]);

title('neural command signal');
xlim([0 tmax]-t_on);
ylim([-1 9]);
legend('pulse', 'step', 'pulse+step', 'Location', 'northeast','FontSize',14);

% (3) Plant outputs: position responses
subplot(1,3,3);
plot(t-t_on,h,'LineWidth', 2); hold on;
plot(t-t_on,U,'LineWidth', 2);
plot(t-t_on,hU,'LineWidth', 2);
xlabel('time re input (s)');
ylabel('eye position (a.u.)');
xlim([0 tmax]-t_on)
ylim([-0.05 0.23]);
set(gca,'YTick',[]);

title('plant response');
legend('pulse \rightarrow fast, no hold', ...
       'step \rightarrow slow, with hold', ...
       'pulse+step \rightarrow fast, with hold', ...
       'Location', 'northeast','FontSize',14);

% % - (4) Text panel summarizing key points -
% subplot(2,2,4);
% axis off;
% text(0, 1, 'Key ideas:', 'FontWeight', 'bold', 'FontSize', 12);
% text(0, 0.8, sprintf('- The plant is linear and overdamped (no oscillations).'));
% text(0, 0.6, sprintf('- T_1 (long) sets the slow tail (elasticity).'));
% text(0, 0.45,sprintf('- T_2 (short) sets the fast initial rise (viscosity).'));
% text(0, 0.25,sprintf('- Step input: slow movement, final position held.'));
% text(0, 0.10,sprintf('- Pulse input: fast movement, but no fixation (eye drifts back).'));
% text(0, -0.05,sprintf('- Pulse+step: fast movement AND stable fixation (saccade).'));
% xlim([-0.05 1]);
% ylim([-0.2 1.1]);

for ii = 1:3
	subplot(1,3,ii)
	if hasNicegraph
		nicegraph;
		set(gca, 'FontSize', 21);
	else
		% Minimal built-in equivalent
		box off;
		set(gca, 'TickDir', 'out', 'FontSize', 21);
	end
if hasBftext
	bf_text(0.05,0.95,char(96+ii),'FontSize',21)
end

end

%% 5. Optionally save figure
outDir = fullfile('..', '..', 'images');
if ~exist(outDir, 'dir')
    mkdir(outDir);
end
fname = fullfile(outDir, 'saccade_oculomotor_plant_fig.png');

if hasSavegraph
    % Use custom helper
    savegraph(fname,'png');
else
    % Built-in export (R2020a+; otherwise use print/saveas)
    exportgraphics(gcf, fname, 'Resolution', 300);
end

end

% ===== Helper function: impulse response of second-order overdamped system =====
function h = impulse_response(t, T1, T2)
%IMPULSE_RESPONSE  Impulse response of second-order overdamped plant.
%
%   h(t) = 1/(T1 - T2) * [exp(-t/T1) - exp(-t/T2)];
%
%   T1 > T2 > 0 are time constants (in seconds)

if T1 ==0

h = -exp(-t./T2);

elseif T2 == 0
h = exp(-t./T1);

else
h = (exp(-t./T1) - exp(-t./T2));
	
end

h(t < 0) = 0;  % ensure causality (not strictly needed here)

end