Index: issm/trunk/src/m/classes/sealevelmodel.m
===================================================================
--- issm/trunk/src/m/classes/sealevelmodel.m	(revision 24313)
+++ issm/trunk/src/m/classes/sealevelmodel.m	(revision 25836)
@@ -12,5 +12,5 @@
 	properties (SetAccess=public) %Model fields
 		% {{{
-		icecaps          = {}; % list of ice cap models
+		icecaps          = {}; % list of land/ice models, name should  be change longer term.
 		earth            = 0;  % model for the whole earth
 		basins           = {}; % list  of basins, matching icecaps, where shapefile info is held.
@@ -19,25 +19,27 @@
 		settings         = 0;
 		private          = 0;
-		range            = 0;
 		mergedcaps       = 0;
 		transitions      = {};
-		eltransitions      = {};
+		eltransitions    = {};
+		planet           = '';
 		%}}}
 	end
 	methods
 		function slm = sealevelmodel(varargin) % {{{
-
-			if nargin==0, 
-				slm=setdefaultparameters(slm);
-			else 
-				slm=setdefaultparameters(slm);
-
-				options=pairoptions(varargin{:}); 
-			
+			slm=setdefaultparameters(slm);
+
+			if nargin==1,
+
+				options=pairoptions(varargin{:});
+
 				%recover all the icecap models: 
 				slm.icecaps=getfieldvalues(options,'ice_cap',{}); 
 				
 				%recover the earth model:
-				slm.earth = getfieldvalue(options,'earth');
+				slm.earth = getfieldvalue(options,'earth',0);
+
+				%set planet type: 
+				slm.planet=getfieldvalue(options,'planet','earth');
+
 			end
 		end
@@ -89,7 +91,7 @@
 			slm.private           = private();
 			slm.cluster           = generic();
-			slm.range             = {};
 			slm.transitions       = {};
-			slm.eltransitions       = {};
+			slm.eltransitions     = {};
+			slm.planet            = 'earth';
 		end
 		%}}}
@@ -100,11 +102,10 @@
 			disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
 			disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
-			disp(sprintf('%19s: %-22s -- %s','range'   ,['[1x1 ' class(self.range) ']'],'ranges'));
 		end % }}}
 		function self=mergeresults(self) % {{{
-			champs=fields(self.icecaps{1}.results.TransientSolution);
+			champs=fieldnames(self.icecaps{1}.results.TransientSolution);
 			for i=1:length(self.mergedcaps)/2,
 				md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2};
-				icecaps=self.icecaps(self.range{2*(i-1)+2});
+				%icecaps=self.icecaps(self.range{2*(i-1)+2});
 				for j=1:length(self.icecaps{1}.results.TransientSolution),
 					for k=1:length(champs),
@@ -143,43 +144,98 @@
 			end
 		end % }}}
+		function n=ncaps(self) % {{{
+			n=length(self.icecaps);
+		end % }}}
+		function list=continents(self) % {{{
+			list={};
+			for  i=1:length(self.basins),
+				list{end+1}=self.basins{i}.continent;
+			end
+			list=unique(list);
+		end % }}}
+		function list=basinsfromcontinent(self,continent) % {{{
+			list={};
+			for  i=1:length(self.icecaps),
+				if strcmpi(self.basins{i}.continent,continent),
+					list{end+1}=self.basins{i}.name;
+				end
+			end
+			list=unique(list);
+		end % }}}
 		function addbasin(self,bas) % {{{
-		if ~strcmpi(class(bas),'basin')
-			error('addbasin method only takes a ''basin'' class object as input');
-		end;
-		self.basins{end+1}=bas;
+			if ~strcmpi(class(bas),'basin')
+				error('addbasin method only takes a ''basin'' class object as input');
+			end;
+			self.basins{end+1}=bas;
+		end % }}}
+		function intersections2d(self,varargin) % {{{
+
+			options=pairoptions(varargin{:}); 
+			force=getfieldvalue(options,'force',0);
+			
+			%initialize, to avoid issues of having more transitions than meshes. 
+			self.transitions={};
+			self.eltransitions={};
+
+			%for elements:
+			xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
+			ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
+			
+			for i=1:length(self.icecaps),
+				mdi=self.icecaps{i};
+		
+				%for elements: 
+				xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
+				yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
+		
+				disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+			
+				self.transitions{end+1}=meshintersect2d(self.earth.mesh.x,self.earth.mesh.y,mdi.mesh.x,mdi.mesh.y,'force',force);
+
+				self.eltransitions{end+1}=meshintersect2d(xe,ye,xei,yei,'force',force);
+			end
 		end % }}}
 		function intersections(self,varargin) % {{{
 
-		options=pairoptions(varargin{:}); 
-		force=getfieldvalue(options,'force',0);
-
-		%for elements:
-		xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
-		ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
-		ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
-
-		for i=1:length(self.icecaps),
-			mdi=self.icecaps{i};
-			mdi=TwoDToThreeD(mdi);
-	
-			%for elements: 
-			xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
-			yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
-			zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
-	
-			disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+			options=pairoptions(varargin{:}); 
+			force=getfieldvalue(options,'force',0);
+			
+			%initialize, to avoid issues of having more transitions than meshes. 
+			self.transitions={};
+			self.eltransitions={};
+
+			%for elements:
+			xe=self.earth.mesh.x(self.earth.mesh.elements)*[1;1;1]/3;
+			ye=self.earth.mesh.y(self.earth.mesh.elements)*[1;1;1]/3;
+			ze=self.earth.mesh.z(self.earth.mesh.elements)*[1;1;1]/3;
+			
+			for i=1:length(self.icecaps),
+				mdi=self.icecaps{i};
+				mdi=TwoDToThreeD(mdi,self.planet);
 		
-			self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
-
-			self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
-		end
-
+				%for elements: 
+				xei=mdi.mesh.x(mdi.mesh.elements)*[1;1;1]/3;
+				yei=mdi.mesh.y(mdi.mesh.elements)*[1;1;1]/3;
+				zei=mdi.mesh.z(mdi.mesh.elements)*[1;1;1]/3;
+		
+				disp(sprintf('Computing vertex intersections for basin %s',self.basins{i}.name));
+			
+				self.transitions{end+1}=meshintersect3d(self.earth.mesh.x,self.earth.mesh.y,self.earth.mesh.z,mdi.mesh.x,mdi.mesh.y,mdi.mesh.z,'force',force);
+
+				self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
+			end
 		end % }}}
 		function checkintersections(self) % {{{
-		flags=zeros(self.earth.mesh.numberofvertices,1);
-		for i=1:length(self.basins),
-			flags(self.transitions{i})=i;
-		end
-		plotmodel(self.earth,'data',flags,'coastline','on');
+			flags=zeros(self.earth.mesh.numberofvertices,1);
+			for i=1:length(self.basins),
+				flags(self.transitions{i})=i;
+			end
+			plotmodel(self.earth,'data',flags,'coastline','on');
+
+		end % }}}
+		function checkbasinconsistency(self) % {{{
+			for i=1:length(self.basins),
+				self.basins{i}.checkconsistency();
+			end
 
 		end % }}}
@@ -192,14 +248,16 @@
 			if iscell(continent),
 				if length(continent)==1,
-					 if strcmpi(continent{1},'all'),
-						 %need to transform this into a list of continents: 
-						 continent={}; 
-						 for i=1:length(self.basins),
-							 continent{end+1}=self.basins{i}.continent;
-						 end
-						 continent=unique(continent);
-					 end
+					if strcmpi(continent{1},'all'),
+						%need to transform this into a list of continents: 
+						continent={}; 
+						for i=1:length(self.basins),
+							continent{end+1}=self.basins{i}.continent;
+						end
+						continent=unique(continent);
+					else
+						%nothing to do: assume we have a list of continents
+					end
 				else
-					%nothing to do, we have a list of continents
+					%nothing to do: assume we have a list of continents
 				end
 			else
@@ -208,5 +266,5 @@
 					continent={}; 
 					for i=1:length(self.basins),
-						 continent{end+1}=self.basins{i}.continent;
+						continent{end+1}=self.basins{i}.continent;
 					end
 					continent=unique(continent);
@@ -219,25 +277,25 @@
 			if iscell(bas),
 				if length(bas)==1,
-					 if strcmpi(bas{1},'all'),
-						 %need to transform this into a list of basins: 
-						 baslist=[];
-						 for i=1:length(self.basins),
-							 if self.basins{i}.iscontinentany(continent{:}),
-								 baslist(end+1)=i;
-							 end
-						 end
-						 baslist=unique(baslist);
-					 else
-						 bas=bas{1};
-						 baslist=[];
-						 for i=1:length(self.basins),
-							 if self.basins{i}.iscontinentany(continent{:}),
-								 if self.basins{i}.isnameany(bas),
-									 baslist(end+1)=i;
-								 end
-							 end
-						 end
-
-					 end
+					if strcmpi(bas{1},'all'),
+						%need to transform this into a list of basins: 
+						baslist=[];
+						for i=1:length(self.basins),
+							if self.basins{i}.iscontinentany(continent{:}),
+								baslist(end+1)=i;
+							end
+						end
+						baslist=unique(baslist);
+					else
+					bas=bas{1};
+					baslist=[];
+					for i=1:length(self.basins),
+						if self.basins{i}.iscontinentany(continent{:}),
+							if self.basins{i}.isnameany(bas),
+								baslist(end+1)=i;
+							end
+						end
+					end
+
+					end
 				else
 					%we have a list of basin names: 
@@ -280,13 +338,13 @@
 		end % }}}
 		function addicecap(self,md) % {{{
-		if ~strcmpi(class(md),'model')
-			error('addicecap method only takes a ''model'' class object as input');
-		end
-		self.icecaps{end+1}=md;
+			if ~strcmpi(class(md),'model')
+				error('addicecap method only takes a ''model'' class object as input');
+			end
+			self.icecaps{end+1}=md;
 		end % }}}
 		function basinsplot3d(self,varargin) % {{{
-		for i=1:length(self.basins),
-			self.basins{i}.plot3d(varargin{:});
-		end
+			for i=1:length(self.basins),
+				self.basins{i}.plot3d(varargin{:});
+			end
 		end % }}}
 		function caticecaps(self,varargin) % {{{
@@ -300,5 +358,5 @@
 			models=self.icecaps;
 			for i=1:length(models),
-				models{i}=TwoDToThreeD(models{i});
+				models{i}=TwoDToThreeD(models{i},self.planet);
 			end
 			
@@ -326,8 +384,40 @@
 			self.earth=md;
 
+			%Create mesh radius:
+			self.earth.mesh.r=planetradius('earth')*ones(md.mesh.numberofvertices,1);
+
+		end % }}}
+		function caticecaps2d(self,varargin) % {{{
+
+			%recover options: 
+			options=pairoptions(varargin{:}); 
+			tolerance=getfieldvalue(options,'tolerance',1e-5);
+			loneedgesdetect=getfieldvalue(options,'loneedgesdetect',0);
+			models=self.icecaps;
+
+			%Plug all models together:
+			md=models{1}; 
+			for i=2:length(models),
+				md=modelmerge2d(md,models{i},'tolerance',tolerance);
+			end
+
+			%Look for lone edges if asked for it: {{{
+			if loneedgesdetect,
+				edges=loneedges(md);
+				hold on;
+				for i=1:length(edges),
+					ind1=edges(i,1);
+					ind2=edges(i,2);
+					plot([md.mesh.x(ind1),md.mesh.x(ind2)],[md.mesh.y(ind1),md.mesh.y(ind2)],'g*-');
+				end
+			end %}}}
+
+			%Plug into earth: 
+			self.earth=md;
+
 		end % }}}
 		function viscousiterations(self) % {{{
 			for  i=1:length(self.icecaps),
-				ic=self.icecaps{i}; 
+				ic=self.icecaps{i};
 				mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
 				for j=2:length(ic.results.TransientSolution)-1,
@@ -339,5 +429,5 @@
 		function maxtimestep(self) % {{{
 			for  i=1:length(self.icecaps),
-				ic=self.icecaps{i}; 
+				ic=self.icecaps{i};
 				mvi=length(ic.results.TransientSolution);
 				timei=ic.results.TransientSolution(end).time;
@@ -350,5 +440,5 @@
 		function transfer(self,string) % {{{
 			%Recover field size in one icecap: 
-			eval(['n=length(self.icecaps{1}.' string ');']);
+			eval(['n=size(self.icecaps{1}.' string ',1);']);
 			if n==self.icecaps{1}.mesh.numberofvertices,
 				eval(['self.earth.' string '=zeros(self.earth.mesh.numberofvertices,1);']);
@@ -356,5 +446,34 @@
 					eval(['self.earth.' string '(self.transitions{' num2str(i) '})=self.icecaps{' num2str(i) '}.' string ';']);
 				end
-			elseif n==self.icecaps{1} .mesh.numberofelements,
+			elseif n==(self.icecaps{1}.mesh.numberofvertices+1),
+				%dealing with a transient dataset. 
+				%check that all timetags are similar between all icecaps:  %{{{
+				for i=1:length(self.icecaps),
+					eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+					for j=(i+1):length(self.icecaps),
+						eval(['capfieldj= self.icecaps{' num2str(j) '}.' string ';']);
+						if ~isequal(capfieldi(end,:),capfieldj(end,:)),
+							error(['Time stamps for ' string ' field are different between icecaps ' num2str(i) ' and ' num2str(j)]);
+						end
+					end
+				end
+				eval(['capfield1= self.icecaps{1}.' string ';']); 
+				times=capfield1(end,:);
+				nsteps=length(times);
+				%}}}
+				%initialize:  %{{{
+				eval(['field=zeros(self.earth.mesh.numberofvertices+1,' num2str(nsteps) ');']);
+				eval(['field(end,:)=times;']); %transfer the times only, not the values
+				%}}}
+				%transfer all time fields: {{{
+				for i=1:length(self.icecaps),
+					eval(['capfieldi= self.icecaps{' num2str(i) '}.' string ';']);
+					for j=1:nsteps,
+						eval(['field(self.transitions{' num2str(i) '},' num2str(j) ')=capfieldi(1:end-1,' num2str(j) ');']); %transfer only the values, not the time.
+					end
+				end
+				eval(['self.earth.' string '=field;']); %do not forget to plug the field variable into its final location
+				%}}}
+			elseif n==self.icecaps{1}.mesh.numberofelements,
 				eval(['self.earth.' string '=zeros(self.earth.mesh.numberofelements,1);']);
 				for i=1:length(self.icecaps),
@@ -389,4 +508,12 @@
 			self.earth=ic;
 		end % }}}
+		function self=initializemodels(self) % {{{
+
+			for i=1:length(self.basins),
+				md=model();
+				md.miscellaneous.name=self.basins{i}.name;
+				self.addicecap(md);
+			end
+		end % }}}
 	end
 end
