function saccade_slow_step_fig
%SACCADE_STEP_FIG  Generate illustration of a saccade as a step response.
%
% This script plots:
%   - the target step (input)
%   - the saccadic eye movement (output)
% in the same axis, and saves the figure as
%   ../../images/saccade_step_fig.png
%
% Intended for use in the Neurobiophysics course (optional download).

%% Initialisation
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'));

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


%% Time axis
tEnd = 1.0;             % total duration [s]
Fs   = 1000;            % sampling rate [Hz]
t    = linspace(0, tEnd, tEnd*Fs + 1);  % time vector

%% Target step (input)
tStep = 0.20;           % time of target jump [s]
x     = double(t >= tStep);  % step from 0 to 1 at t = tStep

%% Saccade simulation (output)
% We model the saccadic eye position as the time integral of a brief
% velocity pulse. The pulse itself is chosen as an asymmetric "gamma-like"
% function; its exact shape is not important, only that it is brief and
% positive and peaks soon after saccade onset and is skewed.

tOnset = 0.20;          % saccade onset (slightly after target step) [s]
tau    = 0.03;         % time scale of the pulse [s]

% Define a velocity pulse v(t)
v = ((t - tOnset) ./ tau).^3 .* exp(-(t - tOnset) ./ tau);  % gamma-like
v(t < tOnset) = 0;      % no velocity before onset

% Integrate velocity to get position, and normalize to step size 1
y = cumtrapz(t, v);
% y = v;
y = y ./ max(y);        % normalize so that final position is ~1

%% Gaussian velocity pulse (slow saccade)
% Parameters controlling the Gaussian
tOnset = 0.23;        % saccade onset [s]
mu     = tOnset + 0.4;   % center of the Gaussian pulse [s]
sigma  = 0.2;           % spread (larger = slower saccade)
A      = 300;             % amplitude of velocity pulse [deg/s]

% Gaussian velocity
v = exp(-0.5 * ((t - mu) ./ sigma).^2);

% Zero velocity before onset
v(t < tOnset) = 0;

% Integrate velocity to obtain position
y_gauss = cumtrapz(t, v);

% Normalize to final amplitude 1 (optional)
y_gauss = y_gauss ./ max(y_gauss);
%% Make figure
figure(1);
clf;
hold on;


% Plot target step (input)
plot(t, x, 'Color', 'k', 'LineWidth', 2);

% Plot saccade (output)
plot(t, y, 'k:','Color', col(2,:), 'LineWidth', 1.5);
plot(t, y_gauss, 'Color', col(3,:), 'LineWidth', 2);

% Axes formatting
% xlim([0 tEnd]);
% ylim([-0.15 1.2]);

ylim([-0.5 1.5]);
xlim([-0.1 1.1])
if hasNicegraph
    nicegraph;
	axis normal;
    set(gca, 'FontSize', 24);
	axis off;
else
    % Minimal built-in equivalent
    box off;
	axis off;
    set(gca, 'TickDir', 'out', 'FontSize', 24);
end

% Manual "axes" labels (for didactic clarity)
text(0, -0.12, 'time', ...
    'HorizontalAlignment', 'left', ...
    'VerticalAlignment',   'top', ...
    'FontSize',            18, ...
    'Color',               'k');

text(-0.05, 0, 'position', ...
    'HorizontalAlignment', 'left', ...
    'VerticalAlignment',   'bottom', ...
    'FontSize',            18, ...
    'Color',               'k', ...
    'Rotation',            90);

% Annotations / labels
text(0.22, 1.1, {'step command','step-like motoneuron activity'}, ...
    'HorizontalAlignment', 'left', ...
    'FontSize',            20, ...
    'Color',               'k');

text(0.35, 0.50, {'actual saccade'}, ...
    'HorizontalAlignment', 'left', ...
    'FontSize',            20, ...
    'Color',               col(2,:));
annotation(gcf,'arrow',[0.33 0.305]+0.09,...
	[0.5 0.5]+0.015);


text(0.68, 0.50, {'step response plant'}, ...
    'HorizontalAlignment', 'left', ...
    'FontSize',            20, ...
    'Color',               col(3,:));
annotation(gcf,'arrow',[0.33 0.305]+0.30,...
	[0.5 0.5]+0.015);
% Optional arrows to indicate step input (coordinates are normalized)

annotation(gcf,'arrow',[0.17 0.3],...
	[0.27 0.27]);
annotation(gcf,'arrow',[0.17 0.17],...
	[0.27 0.45]);
%% Save figure
fname = fullfile('..','..','images','saccade_slow_step_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