Index: /issm/trunk/src/m/utils/Exp/downstreamflowlines.m
===================================================================
--- /issm/trunk/src/m/utils/Exp/downstreamflowlines.m	(revision 9360)
+++ /issm/trunk/src/m/utils/Exp/downstreamflowlines.m	(revision 9360)
@@ -0,0 +1,110 @@
+function flowpath=downstreamflowlines(index,x,y,u,v,x0,y0,varargin)
+%DOWNSTREAMFLOWLINES - compute flowlines from a given set of seed points
+%
+%   Usage:
+%      flowpath=downstreamflowlines(index,x,y,u,v,x0,y0)
+%
+%   the velocity field is given by the couple (u,v) and the coordinates
+%   of the seed points are (x0,y0). One can use one or several seed 
+%   points
+%
+%   Example:
+%      flowpath=downstreamflowlines(md.elements,md.x,md.y,md.vx,md.vy,x0,y0)
+
+%check input size
+if nargin>9 | nargin<7,
+	help flowlines
+	error('flowlines error message: bad usage');
+end
+
+%check input
+if (length(x)~=length(y) | length(x)~=length(u) | length(x)~=length(v)),
+	error('flowlines error message: x,y,u and v must have the same length');
+end
+if length(x)<3,
+	error('flowlines error message: at least one element is required');
+end
+if length(x0)~=length(y0),
+	error('flowlines error message: x0 and y0 do not have the same length');
+end
+
+%get maxiter and precision
+if nargin==9
+	maxiter=varargin{1};%maximum number of iterations
+	precision=varargin{2}; %division of each segment (higer precision increases number of segments)
+else
+	maxiter=200; %maximum number of iterations
+	precision=1; %division of each segment (higer precision increases number of segments)
+end
+
+%check seed points
+tria=TriaSearch(index,x,y,x0,y0);
+pos=find(isnan(tria));
+x0(pos)=[];
+y0(pos)=[];
+
+%initialize other variables
+N=length(x0);
+X=x0; Y=y0;
+flowpath=struct('x',cell(N,1),'y',cell(N,1),'name','','density',1);
+for i=1:N,
+	flowpath(i).x=x0(i);
+	flowpath(i).y=y0(i);
+end
+done=zeros(N,1);
+
+%get avegared length of each element
+length_tria=1/3*(sqrt( (x(index(:,1))-x(index(:,2))).^2+(y(index(:,1))-y(index(:,2))).^2 )+...
+	sqrt((x(index(:,1))-x(index(:,3))).^2+(y(index(:,1))-y(index(:,3))).^2 )+...
+	sqrt((x(index(:,2))-x(index(:,3))).^2+(y(index(:,2))-y(index(:,3))).^2 ));
+
+%take velocity for each element
+u=u(index)*[1;1;1]/3;
+v=v(index)*[1;1;1]/3;
+
+%initialization:
+counter=1;
+
+while any(~done) 
+
+	%find current triangle
+	queue=find(~done);
+	tria=TriaSearch(index,x,y,X(queue),Y(queue));
+
+	%check that the point is actually inside a triangle of the mesh
+	listnan=find(isnan(tria));
+	for i=1:length(listnan)
+		%remove the last point
+		flowpath(queue(listnan(i))).x(end)=[];
+		flowpath(queue(listnan(i))).y(end)=[];
+		done(queue(listnan(i)))=1;
+	end
+	tria(listnan)=[]; 
+	queue(listnan)=[];
+
+	if isempty(tria),
+		break;
+	end
+
+	%velocity of the current triangle and norm it
+	ut=u(tria); vt=v(tria); normv=sqrt(ut.^2+vt.^2);
+	ut=ut./normv;vt=vt./normv;
+
+	%check counter
+	if counter>maxiter
+		disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going forward'])
+		break
+	end
+	counter=counter+1
+
+	%remove stagnant point
+	done(queue(find(ut==0 & vt==0)))=1;
+
+	%build next point
+	for i=1:length(queue)
+		X(queue(i))=flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision;
+		Y(queue(i))=flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision;
+		flowpath(queue(i)).x=[flowpath(queue(i)).x;flowpath(queue(i)).x(end)+ut(i)*length_tria(tria(i))/precision];
+		flowpath(queue(i)).y=[flowpath(queue(i)).y;flowpath(queue(i)).y(end)+vt(i)*length_tria(tria(i))/precision];
+	end
+end