Index: /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24845)
+++ /issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m	(revision 24846)
@@ -142,4 +142,7 @@
 			
 			for gid=1:length(self.glacier_connectivity),
+				if mod(gid,1000)==0,
+					disp(num2str(gid/length(self.glacier_connectivity)*100));
+				end
 				el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));
 				x1=latel(1); x2=latel(2); x3=latel(3); 
@@ -151,54 +154,19 @@
 				
 				if ~isintriangle(x0,x1,x2,x3,y0,y1,y2,y3),
-					%plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-'); text(x1,y1,'1'); text(x2,y2,'2'); text(x3,y3,'3'); plot(x0,y0,'r*');
-					%error(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
-					disp(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
-					counter=counter+1
-				end
-			end
-
-		end % }}}
-		function connectivity_check2(self,md) % {{{
-			disp(sprintf('   Checking glacier connectivity:')); 
-			figure(1),clf,hold on;
-			
-			%some quick checks:  % {{{
-			if ~isempty(find(self.glacier_connectivity>md.mesh.numberofelements)),
-				error('glacier_connectivity table is pointing to elements that are > md.mesh.numberofelements!');
-			end
-			if ~isempty(find(self.glacier_connectivity<0)),
-				error('glacier_connectivity table is pointing to elements that are < 0!');
-			end  %}}}
-			
-			for gid=1:length(self.glacier_connectivity),
-				el=self.glacier_connectivity(gid); latel=md.mesh.lat(md.mesh.elements(el,:)); longel=md.mesh.long(md.mesh.elements(el,:));
-				[lid,rid]=self.gidtolid(gid);
-				lat=self.regions(rid).CenLat(lid);
-				long=self.regions(rid).CenLon(lid);
-				
-				proj=laea(lat,long);
-				[longel;long]
-				[x,y]=gdaltransform([longel;long],[latel;lat],'EPSG:4326',proj);
-				x1=x(1); y1=y(1); x2=x(2); y2=y(2); x3=x(3); y3=y(3); x0=x(4); y0=y(4);
-				
-				%test if x0,y0 is in the 1,2,3 triangle: 
-				epsilon=10^-10;
-				delta=x2*y3-y2*x3-x1*y3+y1*x3+x1*y2-y1*x2;
-
-				%first area coordinate
-				area_1=((y3-y2)*(x0-x3)-(x3-x2)*(y0-y3))./delta
-				%second area coordinate
-				area_2=((x3-x1)*(y0-y3)-(y3-y1)*(x0-x3))./delta
-				%third area coordinate
-				area_3=1-area_1-area_2
-
-				isintriangle=((area_1>=0-epsilon) & (area_2>=0-epsilon) & (area_3>=0-epsilon));
-
-				plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-');
-				plot(x0,y0,'r*');
-				if ~isintriangle,
-					error();
-				end
-
+					
+					%check in local coordinate system: 
+					proj=laea(x0,y0);
+					[xp,yp]=gdaltransform([y0,y1,y2,y3],[x0,x1,x2,x3],'EPSG:4326',proj);
+					x0=xp(1); x1=xp(2); x2=xp(3); x3=xp(4);
+					y0=yp(1); y1=yp(2); y2=yp(3); y3=yp(4);
+				
+					if ~isintriangle(x0,x1,x2,x3,y0,y1,y2,y3),
+					
+						plot([x1,x2,x3,x1],[y1,y2,y3,y1],'k*-'); text(x1,y1,'1'); text(x2,y2,'2'); text(x3,y3,'3'); plot(x0,y0,'r*');
+						error(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
+						disp(sprintf('Vertex %i in region %i and lid %i is not within element %i',gid,rid,lid,el));
+						counter=counter+1
+					end
+				end
 			end