% Code to extract firing threholds for VelPos experiment
% Note Eye velocity is in °/frame
tic
SavePath='ChangeAccordingly';
SaveStructP=cell(1,889);
SaveStructV=cell(1,889);

for k=1:size(meshfitting.mesh,3)
    meshdata=meshfitting.mesh(:,:,k);
    exportmatdata=meshfitting.exportmat(k,:);
    rawdata=meshfitting.RawData(:,:,k);

    
    %% Position 
    
    
    PosTuning=nanmean(meshdata(7:8,:));
    %get peak orientation
    findpeak=PosTuning';
    findpeak(:,2)=1:30;
    findpeak(isnan(findpeak(:,1)),:)=[];
    findpeak=sortrows(findpeak,1);
    running=0;
    for pfk=0:2
        if findpeak(end-pfk,2)<mean(findpeak(:,2))
            running=running-1; %remove one when peak is to the left
        else
            running=running+1; %remove one when peak is to the left
        end
    end
    if running>0
        peak=1;
    elseif running<0
        peak=0;
    else
        FORCEDERROR
    end
    clear running findpeak pfk 
    %get firing threshold
    Thres=NaN;
    ThresIndex=NaN;
    
    if peak==0 %check the firing threshold dependet on the side of activity
        iterbasesta=isnan(PosTuning);
        iterbasesta2=find(iterbasesta==0,1,'last');
        pvaluesIter=nan(size(PosTuning,2),1); 
        hvaluesIter=nan(size(PosTuning,2),1);
        for g=1:(iterbasesta2-2-(find(iterbasesta==0,1,'first'))) %highest x value with big enough sample size (iterbasesta2))  - 2 (from initial Baseline) - lowest xvalue with enough sample size = number of times I can run ranksum
            iterbase=vertcat(reshape(cell2mat(rawdata(7,iterbasesta2-1-g:iterbasesta2)),[],1),... % highest x value with big enough sample size -1 -g (starts with 1)= (for g=1) 3 samples as baseline + increasing with loop 
                reshape(cell2mat(rawdata(8,iterbasesta2-1-g:iterbasesta2)),[],1));
        try
            [pvaluesIter(iterbasesta2-2-g,1) hvaluesIter(iterbasesta2-2-g,1)]=...
                ranksum(vertcat(reshape(cell2mat(rawdata(7,iterbasesta2-2-g)),[],1),reshape(cell2mat(rawdata(8,iterbasesta2-2-g)),[],1)),...
                iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline  
        catch
            pvaluesIter(iterbasesta2-2-g,1)=NaN;
            hvaluesIter(iterbasesta2-2-g,1)=NaN;
        end
        end
        
    elseif peak==1
        iterbasesta=isnan(PosTuning);
        iterbasesta2=find(iterbasesta==0,1,'first');
        pvaluesIter=nan(size(PosTuning,2),1); 
        hvaluesIter=nan(size(PosTuning,2),1);
        for g=1:(find(iterbasesta==0,1,'last')-iterbasesta2-2) %highest x value with big enough sample size  - 2 (from initial Baseline) - %lowest xvalue with enough sample size = number of times I can run ranksum
            iterbase=vertcat(reshape(cell2mat(rawdata(7,iterbasesta2:iterbasesta2+1+g)),[],1),...
                reshape(cell2mat(rawdata(8,iterbasesta2:iterbasesta2+1+g)),[],1));
            
            try
            [pvaluesIter(iterbasesta2+2+g,1) hvaluesIter(iterbasesta2+2+g,1)]=...
            ranksum(vertcat(reshape(cell2mat(rawdata(7,iterbasesta2+2+g)),[],1),reshape(cell2mat(rawdata(8,iterbasesta2+2+g)),[],1)),...
            iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline  
            catch
                pvaluesIter(iterbasesta2+2+g,1)=NaN;
                hvaluesIter(iterbasesta2+2+g,1)=NaN;
            end
   
        end
    else
        Forcederror
    end
    
    for kk=1:size(hvaluesIter,1) % set 0 as NaN's so that only significant values := 1
        if hvaluesIter(kk,1)==0
            hvaluesIter(kk,1)=NaN;
        end
        if hvaluesIter(kk,1)==1
            hvaluesIter(kk,1)=PosTuning(1,kk);
        end
    end
    clear g iterbase iterbasesta iterbasesta2 kk 
    
    
    %find firing threhsold
    if peak==1 %peak right     
    indexvector=find(~isnan(hvaluesIter));
    threshrunning=0;
    for ivk=1:size(indexvector,1)
     if ivk~=size(indexvector,1) %significant if at least one or both of the following points are also significant
        if (ismember(indexvector(ivk)+1,indexvector) ||...
            ismember(indexvector(ivk)+2,indexvector)) &&...
            threshrunning==0
        indexvectorEyePos=-14.5:1:14.5;
        Thres=indexvectorEyePos(indexvector(ivk)); %get threshold
        threshrunning=threshrunning+1;
        ThresIndex=indexvector(ivk);
        end
     elseif ivk==size(indexvector,1)  %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
         if threshrunning==0 && length(indexvector)>=3
            indexvectorEyePos=-14.5:1:14.5;
            Thres=indexvectorEyePos(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndex=indexvector(ivk);
         end
     else
         FORCEDERROR
     end
    end
    elseif peak==0 %peak left
    indexvector=find(~isnan(hvaluesIter));
    threshrunning=0;
    for ivk=size(indexvector,1):-1:1
     if ivk~=1
        if (ismember(indexvector(ivk)-1,indexvector) ||... %significant if at least one or both of the following points are also significant
            ismember(indexvector(ivk)-2,indexvector)) &&...
            threshrunning==0
            indexvectorEyePos=-14.5:1:14.5;
            Thres=indexvectorEyePos(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndex=indexvector(ivk);
        end
     elseif ivk==1
        if threshrunning==0 && length(indexvector)>=3 %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
            indexvectorEyePos=-14.5:1:14.5;
            Thres=indexvectorEyePos(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndex=indexvector(ivk);
        end
     else 
         FORCEDERROR
     end
    end        
    else
        FORCEDERROR
    end
    clear indexvector indexvectorEyePos ivk threshrunning
  

    %% Now for bootstrap
    for bootsk=1:1000
    %Boots Position
    ThresBo=NaN;
    ThresIndexBo=NaN;
    if peak==0 %peak left
        iterbasestaBo=isnan(PosTuning);
        iterbasesta2Bo=find(iterbasestaBo==0,1,'last');
        pvaluesIterBo=nan(size(PosTuning,2),1); 
        hvaluesIterBo=nan(size(PosTuning,2),1);
        for gBo=1:(iterbasesta2Bo-2-(find(iterbasestaBo==0,1,'first')))
            % get baseline
            baselineBo=getBaselineBoPosNeg(iterbasesta2Bo-1-gBo:iterbasesta2Bo,rawdata);
            clear newbin
            newbinVal=GetRawDataPosNeg(iterbasesta2Bo-2-gBo,rawdata);
            try
            [pvaluesIterBo(iterbasesta2Bo-2-gBo,1) hvaluesIterBo(iterbasesta2Bo-2-gBo,1)]=...
                ranksum(newbinVal,baselineBo, 'tail', 'right','alpha', 0.05/gBo);
            catch
                pvaluesIterBo(iterbasesta2Bo-2-gBo,1)=NaN;
                hvaluesIterBo(iterbasesta2Bo-2-gBo,1)=NaN;
            end
        end % end gBo
        
    elseif peak==1 %peak right
        iterbasestaBo=isnan(PosTuning);
        iterbasesta2Bo=find(iterbasestaBo==0,1,'first');
        pvaluesIterBo=nan(size(PosTuning,2),1); 
        hvaluesIterBo=nan(size(PosTuning,2),1);
        for gBo=1:(find(iterbasestaBo==0,1,'last')-iterbasesta2Bo-2) 
            baselineBo=getBaselineBoPosNeg(iterbasesta2Bo:iterbasesta2Bo+1+gBo,rawdata); %Pos and Neg peak should not make difference
            clear newbin
            newbinVal=GetRawDataPosNeg(iterbasesta2Bo+2+gBo,rawdata); %Pos and Neg peak should not make difference
            try 
                [pvaluesIterBo(iterbasesta2Bo+2+gBo,1) hvaluesIterBo(iterbasesta2Bo+2+gBo,1)]=...
                    ranksum(newbinVal,baselineBo, 'tail', 'right','alpha', 0.05/gBo);
            catch 
                pvaluesIterBo(iterbasesta2Bo+2+gBo,1)=NaN;
                hvaluesIterBo(iterbasesta2Bo+2+gBo,1)=NaN;
            end
        end
        
    else
        FORCEDERROR
    end
      
    
    
    for kkBo=1:size(hvaluesIterBo,1) % set 0 as NaN's so that only significant values := 1
        if hvaluesIterBo(kkBo,1)==0
            hvaluesIterBo(kkBo,1)=NaN;
        end
        if hvaluesIterBo(kkBo,1)==1
            hvaluesIterBo(kkBo,1)=PosTuning(1,kkBo);
        end
    end
    clear gBo iterbaseBo iterbasestaBo iterbasesta2Bo kkBo 
        
     %find firing threhsold
    if peak==1 %peak right     
    indexvectorBo=find(~isnan(hvaluesIterBo));
    threshrunningBo=0;
    for ivkBo=1:size(indexvectorBo,1)
     if ivkBo~=size(indexvectorBo,1) %significant if at least one or both of the following points are also significant
        if (ismember(indexvectorBo(ivkBo)+1,indexvectorBo) ||...
            ismember(indexvectorBo(ivkBo)+2,indexvectorBo)) &&...
            threshrunningBo==0
        indexvectorEyePos=-14.5:1:14.5;
        ThresBo=indexvectorEyePos(indexvectorBo(ivkBo)); %get threshold
        threshrunningBo=threshrunningBo+1;
        ThresIndexBo=indexvectorBo(ivkBo);
        end
     elseif ivkBo==size(indexvectorBo,1)  %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
         if threshrunningBo==0 && length(indexvectorBo)>=3
            indexvectorEyePos=-14.5:1:14.5;
            ThresBo=indexvectorEyePos(indexvectorBo(ivkBo));
            threshrunningBo=threshrunningBo+1;
            ThresIndexBo=indexvectorBo(ivkBo);
         end
     else
         FORCEDERROR
     end
    end
    elseif peak==0 %peak left
    indexvectorBo=find(~isnan(hvaluesIterBo));
    threshrunningBo=0;
    for ivkBo=size(indexvectorBo,1):-1:1
     if ivkBo~=1
        if (ismember(indexvectorBo(ivkBo)-1,indexvectorBo) ||... %significant if at least one or both of the following points are also significant
            ismember(indexvectorBo(ivkBo)-2,indexvectorBo)) &&...
            threshrunningBo==0
            indexvectorEyePos=-14.5:1:14.5;
            ThresBo=indexvectorEyePos(indexvectorBo(ivkBo));
            threshrunningBo=threshrunningBo+1;
            ThresIndexBo=indexvectorBo(ivkBo);
        end
     elseif ivkBo==1
        if threshrunningBo==0 && length(indexvectorBo)>=3 %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
            indexvectorEyePos=-14.5:1:14.5;
            ThresBo=indexvectorEyePos(indexvectorBo(ivkBo));
            threshrunningBo=threshrunningBo+1;
            ThresIndexBo=indexvectorBo(ivkBo);
        end
     else 
         FORCEDERROR
     end
    end        
    else
        FORCEDERROR
    end
    
    clear indexvectorBo indexvectorEyePosBo ivkBo threshrunningBo
    
  
    SaveStructP{1,k}(1,bootsk)=ThresBo;
    
    end %bootsend

    % end boots Pos
    %% velocity
   
    VelTuning=nanmean(meshdata(:,15:16),2);
    %get peak orientation
    findpeakVel=VelTuning;
    findpeakVel(:,2)=1:14;
    findpeakVel(isnan(findpeakVel(:,1)),:)=[];
    findpeakVel=sortrows(findpeakVel,1);
    running=0;
    for pfk=0:2
        if findpeakVel(end-pfk,2)<mean(findpeakVel(:,2))
            running=running-1; %remove one when peak is to the left
        else
            running=running+1; %add one when peak is to the right
        end
    end
    if running>0
        peakVel=1;
    elseif running<0
        peakVel=0;
    else
        FORCEDERROR
    end
    clear running findpeakVel pfk 
    
    
    %get firing threshold
    ThresVel=NaN;
    ThresIndexVel=NaN;
    if peakVel==0
        iterbasesta=isnan(VelTuning);
        iterbasesta2=find(iterbasesta==0,1,'last');
        pvaluesIterVel=nan(size(VelTuning,1),1); 
        hvaluesIterVel=nan(size(VelTuning,1),1);
        for g=1:(iterbasesta2-2-(find(iterbasesta==0,1,'first'))) %highest x value with big enough sample size (iterbasesta2))  - 2 (from initial Baseline) - lowest xvalue with enough sample size = number of times I can run ranksum
            iterbase=vertcat(reshape(cell2mat(rawdata(iterbasesta2-1-g:iterbasesta2,15)'),[],1),... % highest x value with big enough sample size -1 -g (starts with 1)= (for g=1) 3 samples as baseline + increasing with loop 
                reshape(cell2mat(rawdata(iterbasesta2-1-g:iterbasesta2,16)'),[],1));
        try
            [pvaluesIterVel(iterbasesta2-2-g,1) hvaluesIterVel(iterbasesta2-2-g,1)]=...
                ranksum(vertcat(reshape(cell2mat(rawdata(iterbasesta2-2-g,15)'),[],1),reshape(cell2mat(rawdata(iterbasesta2-2-g,16)'),[],1)),...
                iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline  
        catch
            pvaluesIterVel(iterbasesta2-2-g,1)=NaN;
            hvaluesIterVel(iterbasesta2-2-g,1)=NaN;
        end
        end
        
    elseif peakVel==1
        iterbasesta=isnan(VelTuning);
        iterbasesta2=find(iterbasesta==0,1,'first');
        pvaluesIterVel=nan(size(VelTuning,1),1); 
        hvaluesIterVel=nan(size(VelTuning,1),1);
        for g=1:(find(iterbasesta==0,1,'last')-iterbasesta2-2) %highest x value with big enough sample size  - 2 (from initial Baseline) - %lowest xvalue with enough sample size = number of times I can run ranksum
            iterbase=vertcat(reshape(cell2mat(rawdata(iterbasesta2:iterbasesta2+1+g,15)'),[],1),...
                reshape(cell2mat(rawdata(iterbasesta2:iterbasesta2+1+g,16)'),[],1));
            
            try
            [pvaluesIterVel(iterbasesta2+2+g,1) hvaluesIterVel(iterbasesta2+2+g,1)]=...
            ranksum(vertcat(reshape(cell2mat(rawdata(iterbasesta2+2+g,15)'),[],1),reshape(cell2mat(rawdata(iterbasesta2+2+g,16)'),[],1)),...
            iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline    
            catch
                pvaluesIterVel(iterbasesta2+2+g,1)=NaN;
                hvaluesIterVel(iterbasesta2+2+g,1)=NaN;
            end
   
        end
    else
        FORCEDERROR
    end    
  
    VelTuning=VelTuning'; %to be consistent with PosTuning 
    
        for kk=1:size(hvaluesIterVel,1) % set 0 as NaN's so that only significant values := 1
        if hvaluesIterVel(kk,1)==0
            hvaluesIterVel(kk,1)=NaN;
        end
        if hvaluesIterVel(kk,1)==1
            hvaluesIterVel(kk,1)=VelTuning(1,kk);
        end
        end
    clear g iterbase iterbasesta iterbasesta2 kk 
         
    %find firing threhsold
    if peakVel==1 %peak right     
    indexvectorVel=find(~isnan(hvaluesIterVel));
    threshrunning=0;
    for ivk=1:size(indexvectorVel,1)
     if ivk~=size(indexvectorVel,1)
        if (ismember(indexvectorVel(ivk)+1,indexvectorVel) ||... %significant if at least one or both of the following points are also significant
            ismember(indexvectorVel(ivk)+2,indexvectorVel)) &&...
            threshrunning==0
        indexvectorEyeVel=-3.25:0.5:3.25;
        ThresVel=indexvectorEyeVel(indexvectorVel(ivk));
        threshrunning=threshrunning+1;
        ThresIndexVel=indexvectorVel(ivk); 
        end
     elseif ivk==size(indexvectorVel,1) 
         if threshrunning==0 && length(indexvectorVel)>=3 %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresVel=indexvectorEyeVel(indexvectorVel(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvectorVel(ivk);
         end
     else
         FORCEDERROR
     end
    end
    elseif peakVel==0 %peak left
    indexvectorVel=find(~isnan(hvaluesIterVel));
    threshrunning=0;
    for ivk=size(indexvectorVel,1):-1:1
     if ivk~=1
        if (ismember(indexvectorVel(ivk)-1,indexvectorVel) ||... %significant if at least one or both of the following points are also significant
            ismember(indexvectorVel(ivk)-2,indexvectorVel)) &&...
            threshrunning==0
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresVel=indexvectorEyeVel(indexvectorVel(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvectorVel(ivk);
        end
     elseif ivk==1
        if threshrunning==0 && length(indexvectorVel)>=3 %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresVel=indexvectorEyeVel(indexvectorVel(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvectorVel(ivk);
        end
     else 
         FORCEDERROR
     end
    end        
    else
        FORCEDERROR
    end

    clear indexvectorVel indexvectorEyeVel ivk threshrunning 

    
    %% Boots Velocity
for kbootsv=1:1000
    ThresBoV=NaN;
    ThresIndexBoV=NaN;
    if peakVel==0 %peak left
        iterbasestaBoV=isnan(VelTuning);
        iterbasesta2BoV=find(iterbasestaBoV==0,1,'last');
        pvaluesIterBoV=nan(size(VelTuning,2),1); 
        hvaluesIterBoV=nan(size(VelTuning,2),1);
        for gBoV=1:(iterbasesta2BoV-2-(find(iterbasestaBoV==0,1,'first')))
%             get baseline
            baselineBoV=getBaselineBoVel(iterbasesta2BoV-1-gBoV:iterbasesta2BoV,rawdata);
            clear newbin
            newbinValV=GetRawDataV(iterbasesta2BoV-2-gBoV,rawdata);
            try
            [pvaluesIterBoV(iterbasesta2BoV-2-gBoV,1) hvaluesIterBoV(iterbasesta2BoV-2-gBoV,1)]=...
                ranksum(newbinValV,baselineBoV, 'tail', 'right','alpha', 0.05/gBoV);
            catch
                pvaluesIterBoV(iterbasesta2BoV-2-gBoV,1)=NaN;
                hvaluesIterBoV(iterbasesta2BoV-2-gBoV,1)=NaN;
            end
        end % end gBo
        
    elseif peakVel==1 %peak right
        iterbasestaBoV=isnan(VelTuning);
        iterbasesta2BoV=find(iterbasestaBoV==0,1,'first');
        pvaluesIterBoV=nan(size(VelTuning,2),1); 
        hvaluesIterBoV=nan(size(VelTuning,2),1);
        for gBoV=1:(find(iterbasestaBoV==0,1,'last')-iterbasesta2BoV-2) 
            baselineBoV=getBaselineBoVel(iterbasesta2BoV:iterbasesta2BoV+1+gBoV,rawdata); %        
            clear newbin
            newbinValV=GetRawDataV(iterbasesta2BoV+2+gBoV,rawdata); %Pos and Neg peak should not make difference
            try 
                [pvaluesIterBoV(iterbasesta2BoV+2+gBoV,1) hvaluesIterBoV(iterbasesta2BoV+2+gBoV,1)]=...
                    ranksum(newbinValV,baselineBoV, 'tail', 'right','alpha', 0.05/gBoV);
            catch 
                pvaluesIterBoV(iterbasesta2BoV+2+gBoV,1)=NaN;
                hvaluesIterBoV(iterbasesta2BoV+2+gBoV,1)=NaN;
            end
        end
        
    else
        FORCEDERROR
    end

    
    for kkBoV=1:size(hvaluesIterBoV,1) % set 0 as NaN's so that only significant values := 1
        if hvaluesIterBoV(kkBoV,1)==0
            hvaluesIterBoV(kkBoV,1)=NaN;
        end
        if hvaluesIterBoV(kkBoV,1)==1
            hvaluesIterBoV(kkBoV,1)=VelTuning(1,kkBoV);
        end
    end
    clear gBoV iterbaseBoV iterbasestaBoV iterbasesta2BoV kkBoV 
        
     %find firing threhsold
    if peakVel==1 %peak right     
    indexvectorBoV=find(~isnan(hvaluesIterBoV));
    threshrunningBoV=0;
    for ivkBoV=1:size(indexvectorBoV,1)
     if ivkBoV~=size(indexvectorBoV,1) %significant if at least one or both of the following points are also significant
        if (ismember(indexvectorBoV(ivkBoV)+1,indexvectorBoV) ||...
            ismember(indexvectorBoV(ivkBoV)+2,indexvectorBoV)) &&...
            threshrunningBoV==0
        indexvectorEyeVel=-3.25:0.5:3.25;
        ThresBoV=indexvectorEyeVel(indexvectorBoV(ivkBoV)); %get threshold
        threshrunningBoV=threshrunningBoV+1;
        ThresIndexBoV=indexvectorBoV(ivkBoV);
        end
     elseif ivkBoV==size(indexvectorBoV,1)  %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
         if threshrunningBoV==0 && length(indexvectorBoV)>=3
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresBoV=indexvectorEyeVel(indexvectorBoV(ivkBoV));
            threshrunningBoV=threshrunningBoV+1;
            ThresIndexBoV=indexvectorBoV(ivkBoV);
         end
     else
         FORCEDERROR
     end
    end
    elseif peakVel==0 %peak left
    indexvectorBoV=find(~isnan(hvaluesIterBoV));
    threshrunningBoV=0;
    for ivkBoV=size(indexvectorBoV,1):-1:1
     if ivkBoV~=1
        if (ismember(indexvectorBoV(ivkBoV)-1,indexvectorBoV) ||... %significant if at least one or both of the following points are also significant
            ismember(indexvectorBoV(ivkBoV)-2,indexvectorBoV)) &&...
            threshrunningBoV==0
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresBoV=indexvectorEyeVel(indexvectorBoV(ivkBoV));
            threshrunningBoV=threshrunningBoV+1;
            ThresIndexBoV=indexvectorBoV(ivkBoV);
        end
     elseif ivkBoV==1
        if threshrunningBoV==0 && length(indexvectorBoV)>=3 %for last point only that one has to be significant if at least 3 points om that tuning curve are significant
            indexvectorEyeVel=-3.25:0.5:3.25;
            ThresBoV=indexvectorEyeVel(indexvectorBoV(ivkBoV));
            threshrunningBoV=threshrunningBoV+1;
            ThresIndexBoV=indexvectorBoV(ivkBoV);
        end
     else 
         FORCEDERROR
     end
    end        
    else
        FORCEDERROR
    end
    
    clear indexvectorBoV indexvectorEyePosBoV ivkBoV threshrunningBoV
    
     SaveStructV{1,k}(1,kbootsv)=ThresBoV;
end
disp('test')

    
    % End Boots Velocity
    %% Plot
    h1=figure(1);
    subplot(2,2,1)
    MyNeuron = meshdata;
    [nr,nc] = size(MyNeuron); %workaround to plot NaN white
    pcolor([MyNeuron nan(nr,1); nan(1,nc+1)]); %add one row and column NaNs to ensure that also the last row/colum is plotted
    shading flat; %removes grid
    set(gca, 'ydir', 'reverse'); %changes origin to lowleft 
    xlabel('Eye position')
    ylabel ('Velocity')
    set(gca,'xtick',[1,9,16,23,31])
    set(gca,'xticklabel',{'-15','-7','0','7','15'})
    set(gca,'ytick',[1,8,15])
    set(gca,'yticklabel',{'-3.5','0','3.5'})
    clear nr nc MyNeuron
    title(strcat('Neuron Number:',num2str(k),'; PVIndex:',num2str(exportmatdata{1,9}(1,16))))
    
    h3=subplot(2,2,3);
    hold on
    eyeposveci=-14.5:1:14.5;
    for gdk=1:30
        clear plotmat
        plotmat=meshdata(:,gdk);
        plotmat(:,2)=eyeposveci(gdk);
        plotmat(:,3)=10;
        scatter(plotmat(:,2),plotmat(:,1),plotmat(:,3),[0.6 0.6 0.6],'filled')
    end
    scatter(-14.5:1:14.5,hvaluesIter,'o','r','filled')
    plot(-14.5:1:14.5,PosTuning(1,:),'k')

    lineyvalH3=get(h3,'YLim');
    if ~isnan(Thres)
        line([Thres Thres],[lineyvalH3(1) lineyvalH3(2)],'LineStyle','--','Color','k')
    end
    xlim([-15,15])
    title('Pos Tuning')
    clear eyeposveci gdk plotmat
    
    
    h4=subplot(2,2,4);
    hold on
    eyevelveci=-3.25:0.5:3.25;
    for gdk=1:14
        clear plotmat
        plotmat=meshdata(gdk,:);
        plotmat=plotmat';
        plotmat(:,2)=eyevelveci(gdk);
        plotmat(:,3)=10;
        scatter(plotmat(:,2),plotmat(:,1),plotmat(:,3),[0.6 0.6 0.6],'filled')
    end
    scatter(-3.25:0.5:3.25,hvaluesIterVel,'o','r','filled')
    plot(-3.25:0.5:3.25,VelTuning,'k')
    lineyvalH4=get(h4,'YLim');
    if ~isnan(ThresVel)
         line([ThresVel ThresVel],[lineyvalH4(1) lineyvalH4(2)],'LineStyle','--','Color','k')
    end
    xlim([-3.5,3.5])
    title('Vel Tuning')
    
    ylim1min=min(lineyvalH3(1),lineyvalH4(1));
    ylim1max=max(lineyvalH3(2),lineyvalH4(2));
    set(h3,'YLim',[ylim1min ylim1max])
    set(h4,'YLim',[ylim1min ylim1max])
    
    
    saveas(h1,fullfile(SavePath,strcat('NeuronNumer_',num2str(k))),'fig')
    saveas(h1,fullfile(SavePath,strcat('NeuronNumer_',num2str(k))),'jpeg')
    clf(h1,'reset')
    
    savematrix(k,1)=exportmatdata{1,9}(1,16);
    savematrix(k,2)=peak;
    savematrix(k,3)=Thres;
    savematrix(k,4)=peakVel;
    savematrix(k,5)=ThresVel;
    
    
    clear lineyvalH3 lineyvalH4 eyeposveci  
    clear meshdata exportmatdata PosTuning h1 h3 h4
    clear hvaluesIter hvaluesIterVel k peak peakVel pvaluesIter pvaluesIterVel rawdata Thres ThresIndexVel ThresVel VelTuning ThresIndex
end
save('savematrix.mat','savematrix')
save('SaveStructP.mat','SaveStructP')
save('SaveStructV.mat','SaveStructV')
disp('batch finished')
toc