Index: sm/trunk-jpl/src/m/exp/downstreamflowlines.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/downstreamflowlines.m	(revision 20534)
+++ 	(revision )
@@ -1,111 +1,0 @@
-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.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.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
-
-%Create triangulation once for all and check seed points
-trep = triangulation(index,x,y);
-tria = pointLocation(trep,[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 = pointLocation(trep,[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
Index: /issm/trunk-jpl/src/m/exp/flowlines.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/flowlines.m	(revision 20534)
+++ /issm/trunk-jpl/src/m/exp/flowlines.m	(revision 20535)
@@ -3,5 +3,5 @@
 %
 %   Usage:
-%      flowpath=flowlines(index,x,y,u,v,x0,y0)
+%      flowpath=flowlines(index,x,y,u,v,x0,y0,varargin)
 %
 %   the velocity field is given by the couple (u,v) and the coordinates
@@ -11,4 +11,10 @@
 %   Example:
 %      flowpath=flowlines(md.mesh.elements,md.mesh.x,md.mesh.y,md.vx,md.initialization.vy,x0,y0)
+%
+%   Options:
+%      - 'maxiter':   how many steps upstream and downstream of the seed points (default: 200)
+%      - 'precision': division of each segment (higer precision increases number of segments, default: 1)
+%      - 'downstream':flow line upstream of the seed points (default: 1)
+%      - 'upstream':  flow line upstream of the seed points (default: 1)
 
 %check input size
@@ -29,11 +35,8 @@
 end
 
-%get maxiter and precision
-if nargin==8
-	maxiter=varargin{1};
-else
-	maxiter=200; %maximum number of iterations
-end
-precision=1; %division of each segment (higer precision increases number of segments)
+%process options
+options   = pairoptions(varargin{:});
+maxiter   = getfieldvalue(options,'maxiter',200);
+precision = getfieldvalue(options,'precision',1);
 
 %Create triangulation once for all and check seed points
@@ -63,102 +66,106 @@
 v=v(index)*[1;1;1]/3;
 
-%initialization:
-counter=1;
+if downstream,
+	%initialization:
+	counter=1;
 
-while any(~done) 
+	while any(~done) 
 
-	%find current triangle
-	queue=find(~done);
-	tria = pointLocation(trep,[X(queue),Y(queue)]);
+		%find current triangle
+		queue=find(~done);
+		tria = pointLocation(trep,[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)=[];
+		%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
+		if isempty(tria),
+			break;
+		end
 
-	%velocity of the current triangle and norm it
-	ut=u(tria); vt=v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));
-	ut=ut./normv;vt=vt./normv;
+		%velocity of the current triangle and norm it
+		ut=u(tria); vt=v(tria); normv=max(eps,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;
+		%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;
+		%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];
+		%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
 end
 
 %same process but reverse (vel=-vel) to have a vcomplete flow line
-queue=[];
-counter=1;
-X=x0; Y=y0;
-done=zeros(N,1);
+if upstream,
+	queue=[];
+	counter=1;
+	X=x0; Y=y0;
+	done=zeros(N,1);
 
-while any(~done) 
+	while any(~done) 
 
-	%find current triangle
-	queue=find(~done);
-	tria = pointLocation(trep,[X(queue),Y(queue)]);
+		%find current triangle
+		queue=find(~done);
+		tria = pointLocation(trep,[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(1)=[];
-		flowpath(queue(listnan(i))).y(1)=[];
-		done(queue(listnan(i)))=1;
-	end
-	tria(listnan)=[]; 
-	queue(listnan)=[];
+		%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(1)=[];
+			flowpath(queue(listnan(i))).y(1)=[];
+			done(queue(listnan(i)))=1;
+		end
+		tria(listnan)=[]; 
+		queue(listnan)=[];
 
-	if isempty(tria),
-		break;
-	end
+		if isempty(tria),
+			break;
+		end
 
-	%velocity of the current triangle and norm it
-	ut=-u(tria); vt=-v(tria); normv=max(eps,sqrt(ut.^2+vt.^2));
-	ut=ut./normv;vt=vt./normv;
+		%velocity of the current triangle and norm it
+		ut=-u(tria); vt=-v(tria); normv=max(eps,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 backward'])
-		break
-	end
-	counter=counter+1;
+		%check counter
+		if counter>maxiter
+			disp(['Maximum number of iterations (' num2str(maxiter) ') reached while going backward'])
+			break
+		end
+		counter=counter+1;
 
-	%remove stagnant point
-	done(queue(find(ut==0 & vt==0)))=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(1)+ut(i)*length_tria(tria(i))/precision;
-		Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision;
-		flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x];
-		flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y];
+		%build next point
+		for i=1:length(queue)
+			X(queue(i))=flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision;
+			Y(queue(i))=flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision;
+			flowpath(queue(i)).x=[flowpath(queue(i)).x(1)+ut(i)*length_tria(tria(i))/precision; flowpath(queue(i)).x];
+			flowpath(queue(i)).y=[flowpath(queue(i)).y(1)+vt(i)*length_tria(tria(i))/precision; flowpath(queue(i)).y];
+		end
 	end
 end
 
-%EXP compatibility
+%EXP compatibility (add name)
 for i=1:length(queue)
 	flowpath(queue(i)).name=['flowline' num2str(i)];