function [xi,yi]=intersection(x1,y1,x2,y2,errorY)
% function [xi,yi]=intersection(x1,y1,x2,y2,errorY)
% Find the intersection between line 1 defined by
% x1, y1 and line 2 defined by x2, y2.
% errorY is the maximum allowed difference between
% the two lines at the reported intersection point.
% If there are multiple intersections, the
% intersection reported is not easily described.
% One line can be vertical, in which case errorY is
% used as error in X.
xsmall=.0001;
line1vertical=0;
line2vertical=0;
n=20;
x1max=max(x1);
x2max=max(x2);
xmax=min(x1max,x2max);
xmaxdata=xmax;
x1min=min(x1);
x2min=min(x2);
diffx1=x1max-x1min;
diffx2=x2max-x2min;
if abs(diffx1)<xsmall; line1vertical=1; end
if abs(diffx2)<xsmall; line2vertical=1; end
xmin=max(x1min,x2min);
xmindata=xmin;
deltaX=xmax-xmin;
dx=deltaX/n;

if line1vertical | line2vertical
	% One of the lines is vertical
	if line1vertical
		% Line 1 is vertical
		% disp('line 1 is vertical')
		if x2max<x1max | x2min>x1max
			disp('line 1 is vertical and the lines do not intersect')
			xi=[];
			yi=[];
		else
			xi=x1max;
			yi=interp1(x2,y2,xi,'linear');
		end
	else
		% Line 2 is vertical
		% disp('line 2 is vertical')
		if x1max<x2max | x1min>x2max
			disp('line 2 is vertical and the lines do not intersect')
			xi=[];
			yi=[];
		else
			xi=x2max;
			yi=interp1(x1,y1,xi,'linear');
		end
	end
else; % One line is not vertical
	xtest=xmin:dx:xmax;
	delta=interp1(x1,y1,xtest,'linear')-interp1(x2,y2,xtest,'linear');
	deltamin=min(abs(delta));
	i=find(abs(delta)==deltamin);
	xi=xtest(i);
	yi=interp1(x1,y1,xi,'linear');
	errortest=abs(deltamin);
	Nabove=find(delta>0);
	Nbelow=find(delta<0);

	if length(Nabove)==length(delta)
		disp(' Line 1 is entirely above line 2. There is no intersection.')
		xi=NaN; yi=NaN;
	elseif length(Nbelow)==length(delta)
		disp(' Line 1 is entirely below line 2. There is no intersection.')	
		xi=NaN; yi=NaN;
	elseif i==1 & errortest<=errorY
		disp(' The intersection is at the start of one of the lines.')
	elseif  i==length(delta) & errortest<=errorY
		disp(' The intersection is at the end of one of the lines.')
	else	
		%**
		counter=0;
		while errortest>=errorY
			counter=counter+1;
			if i==1
				imin=1;
				imax=2;
				timetostop=0;
			elseif i==length(xtest)
				imin=i-1;
				imax=i;
			else
				imin=i-1;
				imax=i+1;
			end	
			xmax=xtest(imax);
			xmin=xtest(imin);
			dx=(xmax-xmin)/n;
			xtest=xmin:dx:xmax;
			delta=interp1(x1,y1,xtest,'linear')-interp1(x2,y2,xtest,'linear');
			deltamin=min(abs(delta));
			i=find(abs(delta)==deltamin);
			xi=xtest(i);
			errortest=abs(deltamin);
			yi=interp1(x1,y1,xi,'linear');
			if counter>=10
				disp('Convergence failed in intersection.m')
				xi=NaN; yi=NaN;				
				break
			end
		end
	end; % end of second if
end; % end of first if
