function [peak, h9, LEBiStimThresEyeVel, REBiStimThresEyeVel, indexvectorEyeVel,...
    output1S1,output1S2,peakLeMono,peakReMono,LELStimThresEyeVel,RERStimThresEyeVel,output1Left,output1Right...
    ]=figure9subfunctionvel(DFFdec, handles)
% to ensure variable is defined
REBiStimThresEyeVel=nan;
LEBiStimThresEyeVel=nan;
indexvectorEyeVel=-4.75:0.5:4.75;
LELStimThresEyeVel=nan;
RERStimThresEyeVel=nan;
peakLeMono=nan;
peakReMono=nan;

% tuning cure for Left eye monocular stimulation
h9=figure(9);
h91=subplot(3,2,1);
toTuningvel=nan(size(DFFdec(1:300),2),2);
toTuningvel(:,2)=DFFdec(1:300);
toTuningvel(2:end,1)=diff(handles.desRfile(1:300,2));
[output1,output2]=averageEye_V3vel(toTuningvel); 
output1(find(output1(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1(:,1),'k')
hold on
toTuningvel(:,2)=zscore(DFFdec(1:300));
toTuningvel(2:end,1)=diff(handles.desRfile(1:300,2));
[output1S1,output2S1]=averageEye_V3vel(toTuningvel);
output1S1(find(output1S1(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1S1(:,1),'r')
title('Left eye monocular stimulation; zscore red')
set(h91,'XLim',[-5 5])
clear toTuningvel

% tuning cure for Right eye monocular stimulation 
h92=subplot(3,2,2);
toTuningvel=nan(size(DFFdec(301:600),2),2);
toTuningvel(:,2)=DFFdec(301:600);
toTuningvel(2:end,1)=diff(handles.desRfile(301:600,3));
[output1,output2]=averageEye_V3vel(toTuningvel);
output1(find(output1(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1(:,1),'k')
hold on
toTuningvel(:,2)=zscore(DFFdec(301:600));
toTuningvel(2:end,1)=diff(handles.desRfile(301:600,3));
[output1S2,output2S2]=averageEye_V3vel(toTuningvel);
output1S2(find(output1S2(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1S2(:,1),'r')
title('right eye monocular stimulation; zscore red')
set(h92,'XLim',[-5 5])
clear toTuningvel

% tuning cure for Left eye binocular stimulation
h93=subplot(3,2,3);
toTuningvel=nan(size(DFFdec(601:end),2),2);
toTuningvel(:,2)=DFFdec(601:end);
toTuningvel(2:end,1)=diff(handles.desRfile(601:end,2));
[output1,output2]=averageEye_V3vel(toTuningvel);
output1(find(output1(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1(:,1),'k')
hold on
toTuningvel(:,2)=zscore(DFFdec(601:end));
toTuningvel(2:end,1)=diff(handles.desRfile(601:end,2));
[output1Left,outputRawLeft]=averageEye_V3vel(toTuningvel);
output1Left(find(output1Left(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1Left(:,1),'r')
output1Left(:,3)=-4.75:0.5:4.75;
title('Left eye binocular stimulation; zscore red')
set(h93,'XLim',[-5 5])
clear output1 toTuningvel

% find peak 1==left
findpeakLeft=output1Left;
findpeakLeft(:,3)=-4.75:0.5:4.75;
findpeakLeft(isnan(findpeakLeft(:,1)),:)=[];
findpeakLeft=sortrows(findpeakLeft,1);
running=0;
for pfk=0:2
    if findpeakLeft(end-pfk,3)<mean(findpeakLeft(:,3))
        running=running-1; %remove one when peak is to the left
    else
        running=running+1; %add 1 when peak is to the right
    end
end

% tuning cure for Right eye binocular stimulation
h94=subplot(3,2,4);
toTuningvel=nan(size(DFFdec(601:end),2),2);
toTuningvel(:,2)=DFFdec(601:end);
toTuningvel(2:end,1)=diff(handles.desRfile(601:end,3));
[output1,output2]=averageEye_V3vel(toTuningvel);
output1(find(output1(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1(:,1),'k')
hold on
toTuningvel(:,2)=zscore(DFFdec(601:end));
toTuningvel(2:end,1)=diff(handles.desRfile(601:end,3));
[output1Right,outputRawRight]=averageEye_V3vel(toTuningvel);
output1Right(find(output1Right(:,2)==0),1)=NaN;
plot(-4.75:0.5:4.75,output1Right(:,1),'r')
output1Right(:,3)=-4.75:0.5:4.75;
title('Right eye binocular stimulation; zscore red')
set(h94,'XLim',[-5 5])


% find peak 1==left
findpeakRight=output1Right;
findpeakRight(:,3)=-4.75:0.5:4.75;
findpeakRight(isnan(findpeakRight(:,1)),:)=[];
findpeakRight=sortrows(findpeakRight,1);
running2=0;
for pfk=0:2
    if findpeakRight(end-pfk,3)<mean(findpeakRight(:,3))
        running2=running2-1; %remove one when peak is to the left
    else
        running2=running2+1; %add 1 when peak is to the right
    end
end

% plot activity and left eye velocity
subplot(3,2,5)
vel_temp1=nan(size(handles.desRfile,1),1);
vel_temp1(2:end,1)=diff(handles.desRfile(:,2));
plot(vel_temp1,'k')
title('eye vel for left eye (dec kernelled DFF red)')
hold on
plot(DFFdec*20,'r')
clear vel_tmep1

% plot activity and right eye velocity
subplot(3,2,6)
vel_temp2=nan(size(handles.desRfile,1),1);
vel_temp2(2:end,1)=diff(handles.desRfile(:,3));
plot(vel_temp2,'k')
hold on
title('eye vel for right eye (dec kernelled DFF red)')
plot(DFFdec*20,'r')

%parameters to save figure
set(h9,'PaperPositionMode','auto')
set(h9, 'units','normalized','outerposition',[0 0 1 1]);


% peak 1 := peak is to the right
% peak 0 := peak is to the left
if running>0 && running2>0
    peak=1;
elseif running<0 && running2<0
    peak=0;
else
    peak=NaN;
end


%% calculate statistical firing threshold for subplot 3 (LE BiStim)

 [LEBiStimThresEyeVel]=getFiringThesholdNewVel(peak,output1Left,outputRawLeft,h93);
 
 
%% calculate statistical firing threshold for subplot 4 (RE BiStim)

 [REBiStimThresEyeVel]=getFiringThesholdNewVel(peak,output1Right,outputRawRight,h94);
% end original subfunction
 
 
if ismember(cell2mat(handles.activeCell(1,67)),30:53) %check if neurons is from new monocular exclusive regressors 
    %define peak for monocular neurons
    % for left eye 
    if ismember(cell2mat(handles.activeCell(1,67)),42:45) %LE Vel
        
        output1S1Mono=output1S1;
        output1S1Mono(:,3)=-4.75:0.5:4.75;
        output1S1Mono(isnan(output1S1Mono(:,1)),:)=[];
        output1S1Mono=sortrows(output1S1Mono,1);
        runningLeMono=0;
        for pfk=0:2
            if output1S1Mono(end-pfk,3)<mean(output1S1Mono(:,3))
                runningLeMono=runningLeMono-1;%remove one when peak is to the left
            else
                runningLeMono=runningLeMono+1;%add 1 when peak is to the right
            end
        end
        clear pfk
        
        if runningLeMono>0
            peakLeMono=1;
        elseif runningLeMono<0
            peakLeMono=0;
        end
        
        [LELStimThresEyeVel]=getFiringThesholdNewVel(peakLeMono,output1S1,output2S1,h91);

    elseif ismember(cell2mat(handles.activeCell(1,67)),50:53) %RE Pos
        
        output1S2Mono=output1S2;
        output1S2Mono(:,3)=-4.75:0.5:4.75;
        output1S2Mono(isnan(output1S2Mono(:,1)),:)=[];
        output1S2Mono=sortrows(output1S2Mono,1);
        runningReMono=0;
        for pfk=0:2
            if output1S2Mono(end-pfk,3)<mean(output1S2Mono(:,3))
                runningReMono=runningReMono-1;%remove one when peak is to the left
            else
                runningReMono=runningReMono+1;%add 1 when peak is to the right
            end
        end
        clear pfk
        
        if runningReMono>0
            peakReMono=1;
        elseif runningReMono<0
            peakReMono=0;
        end
        
       [RERStimThresEyeVel]=getFiringThesholdNewVel(peakReMono,output1S2,output2S2,h92);

    end
end
end %function end