function [square_wave_jerks square_wave_jerks2] = swjs( sacc_starts, sacc_ends, sacc_magnitudes, sacc_directions, blinks, parameters )
%-------------------------------------------------------------------
%
%  FUNCTION swjs.m
%
%   Detects all the SWJs in a sequence of consecutive saccades
%
%   Distinctive Features of Saccadic Intrusions and Microsaccades in
%   Progressive Supranuclear Palsy. Otero-Millan, Leigh, Serra, Troncoso,
%   Macknik, Martinez-Conde. Journal of Neuroscience (under review).
%   Patent pending.
%
%-------------------------------------------------------------------
%
%  INPUT:
%
%  sacc_starts          start indices of the saccades
%  sacc_ends            end indices of the saccades
%  sacc_magnitudes      magnitude of the saccades
%  sacc_directions      direction of the saccades
%  blinks               binary vector indicating for each sample if it
%                       belons to a blink (1) or not (0)
%  paramters            parameters necessary to calculate the swj index
%     parameters.DM     Difference in direction mixed gaussian means [270 270]
%     parameters.DS     Difference in direction mixed gaussian stds [30 7]
%     parameters.DR     Difference in direction mixed gaussian ratios [0.4 0.6]
%     
%     parameters.RM     Magnitude similarity index mixed gaussian means [0 0]
%     parameters.RS     Magnitude similarity index mixed gaussian stds [0.39 0.16]
%     parameters.RR     Magnitude similarity index mixed gaussian ratios [0.4 0.6]
%     parameters.ISIP   Intersaccadic interval ex-gaussian mean, std and tau [125 60 180]
%     parameters.swjindexth     Threshold for the swj index [0.0014]
%
%  OUTPUT:
%
%  square_wave_jerks   	indices of the first saccade of each SWJ
%  square_wave_jerks2   indices of the second saccade of each SWJ
%
%---------------------------------------------------------------------

% find possible square wave jerks

square_wave_jerks = [];
swjindexes = [];

i=0;
while( i < length(sacc_starts)-1 )
    i = i+1;
    [isswj swjindex] = is_possible_swj( i, sacc_starts, sacc_ends, sacc_magnitudes, sacc_directions, blinks, parameters );
    if ( isswj )
        square_wave_jerks(end+1) = i;
        swjindexes(end+1) = swjindex;
    end
end

% find disjoint SWJs
swj = square_wave_jerks;
square_wave_jerks = [];


i=1;
while (i< length(swj) )
    
    % amount of consecutive square wave jerks
    s = sum( swj(i:end) == swj(i)-1+ [1:length(swj(i:end))] );
    
    % if it is an even number I can chose how to get the swjs
    if ( mod(s,2) == 0 )
        %I select the ones that are smaller in time
        %if ( sum(usacc_ends(swj(i:2:i+s-2)+1)   -   usacc_starts(swj(i:2:i+s-2) )) > sum(usacc_ends(swj(i+1:2:i+s-1)+1)   -   usacc_starts(swj(i+1:2:i+s-1) )))
        if sum((swjindexes(i:2:i+s-2))) > sum((swjindexes(i+1:2:i+s-1)))
            square_wave_jerks = [square_wave_jerks swj(i+1:2:i+s-1) ];
        else
            square_wave_jerks = [square_wave_jerks swj(i:2:i+s-2) ];
        end
        i = i+s;
        % if it is an odd number I do not have election
    elseif ( s==1)
        square_wave_jerks = [square_wave_jerks swj(i:2:i+s-1) ];
        i = i+s;
    else
        if sum((swjindexes(i:2:i+s-1))) > sum((swjindexes(i+1:2:i+s-2)))
            square_wave_jerks =  [square_wave_jerks swj(i:2:i+s-1) ];
            i = i+s;
        else
            square_wave_jerks =  [square_wave_jerks swj(i+1:2:i+s-2) ];
            i = i+s;
        end
    end
end
square_wave_jerks = square_wave_jerks';
square_wave_jerks2 = square_wave_jerks +1;


function [result swjindex] = is_possible_swj( sacc_index, sacc_starts, sacc_ends, sacc_magnitudes, sacc_directions, blinks, parameters)
% result = is_possible_swj( sacc_index, sacc_starts, sacc_ends, sacc_magnitudes, sacc_directions, blinks)
%
% given the index of a microsaccade decide if that microsaccade and the next could form an square
% wave jerk
%


result = false;
swjindex = 0;

%% no blinks or saccades in between the two microsaccades
if ( sum( blinks( sacc_starts(sacc_index):sacc_ends(sacc_index+1) ) ) > 0 )
	return
end

relmag = (usacc_magnitudes(usacc_index+1) - usacc_magnitudes(usacc_index)) ./ (usacc_magnitudes(usacc_index+1) + usacc_magnitudes(usacc_index));
dirdif = mod(acos(cos((usacc_directions(usacc_index) - usacc_directions(usacc_index+1))*pi/180))*180/pi .* sign((usacc_directions(usacc_index) - usacc_directions(usacc_index+1)))+90,360);
isi = (usacc_starts(usacc_index+1) - usacc_ends(usacc_index)+1)*2;

if ( ~exist('parameters' ) )
    DM = [270 270];
    DS = [30 7];
    DR = [0.4 0.6];
    
    RM = [0 0];
    RS = [ 0.39 0.16];
    RR = [0.4 0.6];
    
    ISIP = [125 60 180];
else
    DM = parameters.DM;
    DS =  parameters.DS;
    DR =  parameters.DR;
    
    RM =  parameters.RM;
    RS =  parameters.RS;
    RR =  parameters.RR;
    
    ISIP =  parameters.ISIP;
end

f1 = 1-normcdf( abs(dirdif-DM(1)), 0, DS(1))*DR(1)-normcdf( abs(dirdif-DM(2)), 0, DS(2))*DR(2);
f2 = 1-normcdf( abs(relmag-RM(1)), 0, RS(1))*RR(1)-normcdf( abs(relmag-RM(2)), 0, RS(2))*RR(2);
f3 = ( double(isi>=200).*(1-exgausscdf( isi,ISIP))+double(isi<200).*(exgausscdf( isi, ISIP)));

swjindex = f1*f2*f3;

if ( ~exist('parameters' ) )
    result = swjindex > 0.0014;
else
    result = swjindex > parameters.swjindexth;
end
    


