Index: /issm/trunk-jpl/src/m/plot/plotchannels.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotchannels.m	(revision 25778)
+++ /issm/trunk-jpl/src/m/plot/plotchannels.m	(revision 25779)
@@ -4,4 +4,5 @@
 %   Usage:
 %      plotchannels(md,channelarea,options)
+%      plotchannels(md,channeldischarge,options)
 %
 %   List of supported options
@@ -10,13 +11,22 @@
 %      - 'colormap' colormap used (default is 'gray')
 %      - 'linewidth' linewidth (default is 2)
+%      - 'quiver' only used for discharge(default is 2)
 %
 %   Example:
 %      plotchannels(md,md.results.TransientSolution(end).ChannelArea,'min',15,'max',36);
+%      plotchannels(md,md.results.TransientSolution(end).ChannelDischarge,'min',15,'max',36,'quiver',1);
 
 %Process options
 options = pairoptions(varargin{:});
 
+%is quiver?
+isquiver = (getfieldvalue(options,'quiver',0)==1);
+
 %define level
-level = channelarea;
+if isquiver
+	level = abs(channelarea);
+else
+	level = channelarea;
+end
 
 %Some processing
@@ -38,67 +48,88 @@
 
 %Reconstruct edges
-if 0,
-	edges = [md.mesh.edges(:,1:2) md.mesh.edges(:,1)]';
-else
-	tic
-	%Maximum number of edges
-	maxnbf = 3*md.mesh.numberofelements;
-	%Initialize intermediaries
-	edges = zeros(maxnbf,2);
-	%Chaining algorithm
-	head_minv = -1*ones(md.mesh.numberofvertices,1);
-	next_face = zeros(maxnbf,1);
-	%Initialize number of faces
-	nbf = 0;
-	for i=1:md.mesh.numberofelements
-		for j=1:3
-         %Get the two indices of the edge number j of the ith triangle
-         v1 = md.mesh.elements(i,j);
-			if(j==3)
-				v2 = md.mesh.elements(i,1);
-			else
-				v2 = md.mesh.elements(i,j+1);
+% {{{
+tic
+%Maximum number of edges
+maxnbf = 3*md.mesh.numberofelements;
+%Initialize intermediaries
+edges = zeros(maxnbf,3);
+exchange = zeros(maxnbf,1);
+%Chaining algorithm
+head_minv = -1*ones(md.mesh.numberofvertices,1);
+next_face = zeros(maxnbf,1);
+%Initialize number of faces
+nbf = 0;
+for i=1:md.mesh.numberofelements
+	for j=1:3
+		%Get the two indices of the edge number j of the ith triangle
+		v1 = md.mesh.elements(i,j);
+		if(j==3)
+			v2 = md.mesh.elements(i,1);
+		else
+			v2 = md.mesh.elements(i,j+1);
+		end
+		%sort
+		if(v2<v1)
+			v3=v2; v2=v1; v1=v3;
+		end
+		%This edge a priori has not been processed yet
+		exists = false;
+		%Go through all processed faces connected to v1 and check whether we have seen this edge yet
+		e=head_minv(v1);
+		while(e~=-1)
+			if(edges(e,2)==v2)
+				exists = true;
+				break;
 			end
-			%sort
-         if(v2<v1)
-            v3=v2; v2=v1; v1=v3;
-			end
-         %This edge a priori has not been processed yet
-         exists = false;
-			%Go through all processed faces connected to v1 and check whether we have seen this edge yet
-			e=head_minv(v1);
-			while(e~=-1)
-            if(edges(e,2)==v2)
-               exists = true;
-               break;
-				end
-				e=next_face(e);
-			end
-			%If this edge is new, add it to the lists
-         if(~exists)
-            %Update edges
-            edges(nbf+1,1) = v1;
-            edges(nbf+1,2) = v2;
-            %Update chain
-            next_face(nbf+1) = head_minv(v1);
-            head_minv(v1)    = nbf+1;
-            %Increase number of faces
-            nbf=nbf+1;
-			end
+			e=next_face(e);
+		end
+		%If this edge is new, add it to the lists
+		if(~exists)
+			%Update edges
+			edges(nbf+1,1) = v1; %vertex 1
+			edges(nbf+1,2) = v2; %vertex 2
+			edges(nbf+1,3) = i;  %element 1 (ignore second one)
+			if(v1~=md.mesh.elements(i,j)) exchange(nbf+1)=1; end
+			%Update chain
+			next_face(nbf+1) = head_minv(v1);
+			head_minv(v1)    = nbf+1;
+			%Increase number of faces
+			nbf=nbf+1;
 		end
 	end
-	edges = edges(1:nbf,:);
-	toc
 end
 
+edges = edges(1:nbf,:);
+pos = find(exchange);
+v3 = edges(pos,1);
+edges(pos,1) = edges(pos,2);
+edges(pos,2) = v3;
+toc
+% }}}
+
 %Change edges formatting so that plot looks ok
-edges = [edges(:,1:2) edges(:,1)]';
+myedges = [edges(:,1:2) edges(:,1)]';
 
 %Loop over all levels and plot
 for i=1:numcolors
 	pos=find(colorind==i);
-	hprime=plot(md.mesh.x(edges(:,pos)),md.mesh.y(edges(:,pos)),...
-		'-','Color',palette(i,:),'LineWidth',2);
+	hprime=plot(md.mesh.x(myedges(:,pos)),md.mesh.y(myedges(:,pos)),'-','Color',palette(i,:),'LineWidth',2);
 	if i==1; hold on; end
+
+	if isquiver
+		x1 = md.mesh.x(edges(pos,1)); x2 = md.mesh.x(edges(pos,2));
+		y1 = md.mesh.y(edges(pos,1)); y2 = md.mesh.y(edges(pos,2));
+		len = sqrt( (x2-x1).^2 + (y2-y1).^2);
+
+		xq = mean([x1 x2],2);
+		yq = mean([y1 y2],2);
+		tx = sign(channelarea(pos)).*(x2-x1)./len .* len/10;
+		ty = sign(channelarea(pos)).*(y2-y1)./len .* len/10;
+		px = -ty;
+		py = tx;
+
+		%hprime=quiver(xq,yq,tx,ty,'color',palette(i,:),'showarrowhead','on','autoscale','off');
+		num = length(pos);
+		patch('Faces',reshape(1:3*num,num,3),'Vertices',[[xq;xq+(-tx+px);xq+(-tx-px)],[yq;yq+(-ty+py);yq+(-ty-py)]],'FaceColor',palette(i,:),'EdgeColor','none');
+	end
 end
-hold off