function saccade
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% simple saccade generator model adapted from Scudder 1988 and 
% Jurgens et al. 1981
% (see lecture notes for model schematic)
%
% 2 ways to implement the model are presented:
% - using the Matlab control systems toolbox
% - using numerical integration in the time domain
%
% written by G. Blohm (2015)
%--------------------------------------------------------------------------

clear

%% initialize 
dt = 0.001;
t = 0:dt:.5; % time array
deltaE = zeros(size(t));
deltaE(200:end) = 20; % position error (for model)
deltaE0 = zeros(size(t));
deltaE0(100:end) = 20; % position error (for plots) - used to implement saccade latency

T1 = 0.175; % visco-elastic time constant
T2 = 0.013; % inertial time constant
gPS = 600; % pulse generator gain

%% CONTROL SYSTEMS TOOLBOX CODE
% %% definition of model
% numRI = 1; denRI = [1 0]; % resettable integrator
% numPG = [gPS]; denPG = 1; % pulse generator (simplest form)
% [numCOM, denCOM] = feedback(numPG, denPG, numRI, denRI); % desired eye velocity command
% 
% [numPS, denPS] = parallel(1, [1 0], T1, 1); % pulse-step generation
% numE = 1; denE = [(T1*T2) (T1+T2) 1]; % eye plant
% [numPM, denPM] = series(numPS, denPS, numE, denE); % pre-motor circuit
% 
% [numS, denS] = series(numCOM, denCOM, numPM, denPM); % full model
% 
% %% simulation
% y = lsim(numS, denS, deltaE, t); % numerical integration of linear system

%% NUMERICAL INTEGRATION CODE
y(1:2) = 0; x2(1:2) = 0; Eint(1:2) = 0; x1(1:2) = 0; dEstar(1:2) = 0; Edot(1:2) = 0; Eint(1:2) = 0;
dEstar = zeros(length(t)); 
gRI = 1; % resettable integrator gain
for i = 1:length(t)-1,
    e(i) = deltaE(i) - dEstar(i); % internally estimated saccade error
    dEstar(i+1) = dEstar(i) + dt*gRI*Edot(i); % resettable integrator
    Edot(i+1) = gPS*e(i); % pulse generator (default gain = 600)
    Eint(i+1) = Eint(i) + dt*Edot(i); % neural integrator
    x1(i) = (T1*Edot(i) + Eint(i)); % pulse-step addition
    x2(i+1) = x2(i) + dt/T2*(-x2(i) + x1(i)); % second ocular time constant (inertia)
    y(i+1) = y(i) + dt/T1*(-y(i) + x2(i)); % output after second ocular time constant (visco-elasticity)
end

%% plot results
figure
hold on;
plot(t,deltaE0,'k:')
plot(t(1:end-1),y(1:end-1),'r')
ylabel('angular position (deg)')
