Index: /issm/trunk-jpl/src/m/mesh/meshband.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshband.m	(revision 22938)
+++ /issm/trunk-jpl/src/m/mesh/meshband.m	(revision 22938)
@@ -0,0 +1,92 @@
+function band=meshband(mh,outerdomain,varargin)
+
+	%process options: 
+	options=pairoptions(varargin{:});
+
+	%some checks on the mesh: 
+	if (isempty(mh.x) | isempty(mh.y)  ) 
+		error('meshband error message: mesh has one of the following empty: ''x'',''y''');
+	end
+
+	%give ourselves a unique temporary directory: 
+	temproot=tempname; mkdir(temproot);
+
+	%align external segments on internal mesh: 
+	mh.segments=alignsegments(mh.segments);
+
+	%figure out domain outline: 
+	meshtodomain(mh,[temproot '/innerdomain.exp']);
+
+	%create domain outine from inner and outer domain
+	inner=expread([temproot '/innerdomain.exp']);
+	if inner.closed==0,
+		inner.nods=inner.nods+1;
+		inner.x(end+1)=inner.x(1);
+		inner.y(end+1)=inner.y(1);
+		inner.closed=1;
+		if getfieldvalue(options,'invert',0),
+			inner.x=flipud(inner.x);
+			inner.y=flipud(inner.y);
+		end
+	end
+
+	[path,name,ext]=fileparts(outerdomain);
+	if strcmpi(ext,'.shp'),
+		outer=shpread(outerdomain);
+	else
+		outer=expread(outerdomain);
+	end
+	if outer.closed==0,
+		outer.nods=outer.nods+1;
+		outer.x(end+1)=outer.x(1);
+		outer.y(end+1)=outer.y(1);
+		outer.closed=1;
+	end
+
+	domain(1).x=outer.x; 
+	domain(1).y=outer.y; 
+	domain(1).density=outer.density; 
+	domain(1).nods=outer.nods; 
+	domain(1).closed=outer.closed;
+	
+	
+	domain(2).x=inner.x; 
+	domain(2).y=inner.y; 
+	domain(2).density=inner.density; 
+	domain(2).nods=inner.nods; 
+	domain(2).closed=inner.closed;
+
+
+	expwrite(domain,[temproot '/Band.exp']);
+	
+	if getfieldvalue(options,'plot',0),
+		figure(1),clf,hold on,axis image;
+		expdisp([temproot '/Band.exp'],'linestyle','k-*');
+	end
+
+	%mesh: 
+	md=bamg(model(),'domain',[temproot '/Band.exp'],'MaxCornerAngle',1e-15,'gradation',10000); band=md.mesh; clear md;
+
+	if getfieldvalue(options,'plot',0),
+		figure(2),clf,trisurf(band.elements,band.x,band.y,band.x),view(2),shading faceted;
+		hold on,expdisp([temproot '/Band.exp'],'linestyle','k-*');
+	end
+
+	%check that the domain vertices = the number of segments: 
+	if abs(length(band.segments) - (length(domain(1).x)+length(domain(2).x)-2))>=2,
+		disp(sprintf('band mesh not consistent: %i!=%i+%i\n',length(band.segments), length(domain(1).x), length(domain(2).x)));
+
+		figure(3),clf,expdisp([temproot '/Band.exp'],'linestyle','r*');
+		hold on; 
+		for i=1:length(band.segments),
+			i1=band.segments(i,1); i2=band.segments(i,2);
+			plot([band.x(i1) band.x(i2)],[band.y(i1) band.y(i2)],'k*-');
+			plot([band.x(i1)+ band.x(i2)]/2,[band.y(i1)+ band.y(i2)]/2,'g*');
+		end
+		%error('band mesh not consistent');
+	end
+
+	%erase temporary directory: 
+	system(['rm -rf ' temproot]);
+
+end
Index: /issm/trunk-jpl/src/m/mesh/meshintersect3d.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshintersect3d.m	(revision 22938)
+++ /issm/trunk-jpl/src/m/mesh/meshintersect3d.m	(revision 22938)
@@ -0,0 +1,33 @@
+function indices=meshintersect3d(x,y,z,xs,ys,zs,varargin)
+%MESHINTERSECT - return indices (into x,y and z) of common values between (x,y,z) and (xs,ys,zs). 
+%  i.e: x(index)=xs; y(index)=ys;
+%
+
+
+	%process options: 
+	options=pairoptions(varargin{:});
+
+	%retrieve tolerance: 
+	tolerance=getfieldvalue(options,'tolerance',10);
+
+	%go through lats,longs and find within tolerance, the index of the corresponding value in lat,long: 
+	indices=zeros(length(xs),1);
+	
+	for i=1:length(xs),
+		distance=sqrt((x-xs(i)).^2+(y-ys(i)).^2+(z-zs(i)).^2);
+		s=find(distance<tolerance);
+		if length(s)>1,
+			for j=1:length(s),
+				hold on;plot3(x(s(j)),y(s(j)),z(s(j)),'c.','MarkerSize',40)
+			end
+			distance(s)
+			error(sprintf('one or more vertices on the global mesh were duplicated (offset %i)',i));
+		elseif isempty(s),
+			plot(distance);
+			min(distance);
+			i
+			error('cannot find concurrent vertics!');
+		else
+			indices(i)=s;
+		end
+	end