function saccade_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.23;          % 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

%% Make figure
figure(1);
clf;
hold on;


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

% Plot saccade (output)
plot(t, y, 'Color', col(2,:), '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.05, 'target / step', ...
    'HorizontalAlignment', 'left', ...
    'FontSize',            20, ...
    'Color',               'k');

text(0.40, 0.50, {'saccade', 'step response'}, ...
    'HorizontalAlignment', 'left', ...
    'FontSize',            20, ...
    'Color',               col(2,:));

% 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_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