function [tNew,yNew,breakIndicesNew] = reduceTimePoints(t,y,breakIndices,desiredPoints)
%REDUCETIMEPOINTS remove time points from a time simulation of a switched
%dynamical system
%
% Input arguments:
%   t:              one dimensional vector of time points
%   y:              output of the dynamical system at the time points t
%   breakIndices:   indices of t where the switching happens
%   desiredPoints:  approximately the number of points after removal
%                   if desiredPoints == 0, then no reduction is performed
%
% Output arguments:
%   tNew:           reduced vector of time points
%   yNew:           reduced output of the dynamical system at time points
%                   tNew
%   breakIndicesNew:the indices of tNew where the switching happens
%
% Authors:
%   Philipp Schulze, pschulze@math.tu-berlin.de
%   Benjamin Unger, unger@math.tu-berlin.de
%
% Version:
%   1.0 (December 22, 2017)
% 
if and(desiredPoints>0,length(t)>desiredPoints)
	switchLengths   = [breakIndices(1) diff(breakIndices)] ;
	fak             = floor(length(t)/desiredPoints) ;
	nSwitchLengths  = floor((switchLengths-1)/fak)+1 ;
    nPoints         = sum(nSwitchLengths) ;
    tNew            = NaN(nPoints,1) ;
    if (size(y,1)>size(y,2)) 
        y = transpose(y);
    end
    yNew            = NaN(size(y,1),nPoints) ;
    breakIndicesNew = NaN(size(breakIndices)) ;
    startInd        = 0;
    startIndNew     = 0;
    
    for i=1:length(breakIndices)
        breakIndicesNew(i) = sum(nSwitchLengths(1:i)) ;
        tNew(startIndNew+1:breakIndicesNew(i)) = t(startInd+1:fak:breakIndices(i)) ;
        yNew(:,startIndNew+1:breakIndicesNew(i)) = y(:,startInd+1:fak:breakIndices(i)) ;
        startInd    = breakIndices(i) ;
        startIndNew = breakIndicesNew(i) ;
    end
    if(abs((tNew(end)-t(end))/t(end))>0.001)
        tNew = [tNew;t(end)] ;
        yNew = [yNew y(end)] ;
        breakIndicesNew(end) = breakIndicesNew(end)+1 ;
    end
else
	tNew = t ;
	yNew = y ;
	breakIndicesNew = breakIndices ;
end

