function yi= interp_linear(xdata,ydata,xi)
	%% linear interpolation 
	% xdata, ydata - known data points
	% xi - points at which we want interpolated values yi
	% WARNING: no checks is done that xi values are inside xdata range

	%% First we need to sort our input vectors
	% we need monotonous grow of x

	% Matlab sort function has extra return value index
	% with the list of indexes so xdata(index) is sorted
	[x,index]=sort(xdata); 
	% we reuse this index to sort 'y' vector appropriately
	y=ydata(index); 


	%% Second we want to calculate parameters of the connecting lines
	Nintervals=length(xdata)-1;
	% for Ndata points we will have Ndata-1 connecting lines
	% each of them will have 2 parameters slope and offset
	% so we need a parameters matrix with Nintervals x 2 values
	p=zeros(Nintervals,2);
	% p(i, 1) is the slope  for the interval between x(i) and x(i+1)
	% p(i, 2) is the offset for the interval between x(i) and x(i+1)
	% so y = offset*x+slope = p1*x+p2 at this interval
	for i=1:Nintervals
		slope = ( y(i+1)-y(i) ) / ( x(i+1)-x(i) );  % slope
		offset = y(i)-slope*x(i);  % offset at x=0
		p(i,1)=slope;
		p(i,2)=offset; 
	end

	%% It is crucial to have yi preinitialized !
	% It significantly speeds up the algorithm, 
	% since computer does not have to reallocate memory for new data points
	% Try to comment the following line and compare the execution time 
	% for a large length of xi
	yi=0.*xi; % simple shortcut to initialize return vector with zeros 
	          % this also takes care of deciding  xdata type (row or column).

    %% Finally we interpolate
	N=length(xi); 
	for i=1:N % we will go through all points to be interpolated
		% let's find nearest left neighbor for xi
		displacement=(xi(i)-x);
		% such neighbor must have smallest possible positive displacement
		[displ_min_positive, interval_index] = min( displacement(displacement >= 0 ));
		if (interval_index > Nintervals ) 
			% xi coincides with the last left most data point
			% we can reuse the last interval parameters
			interval_index = Nintervals;
		end
		% index tells which  interval to use for the linear approximation
		% I will use 'polyval' matlab function
		% it evaluate value of polynom of degree n
		% at point x=p_1*x^n+p_2*x^(n-1)+ p_n*x +p_{n+1}
		% since line is polynomial of degree 1 we can use it
		% yi(i)= p(interval_index,1) * xi(i) +p(interval_index,2);
		poly_coef=p(interval_index,:);
		yi(i)=polyval(poly_coef,xi(i));
		
	end
	plot(x,y,'o', xi,yi,'-');
end


