Index: /issm/trunk/src/m/classes/version/7.7/model.m
===================================================================
--- /issm/trunk/src/m/classes/version/7.7/model.m	(revision 8727)
+++ /issm/trunk/src/m/classes/version/7.7/model.m	(revision 8728)
@@ -26,5 +26,4 @@
 		 bamg=struct();
 		 dim=0;
-		 planet=0;
 		 numberofelements=0;
 		 numberofnodes=0;
Index: /issm/trunk/src/m/classes/version/7.7/planet.m
===================================================================
--- /issm/trunk/src/m/classes/version/7.7/planet.m	(revision 8728)
+++ /issm/trunk/src/m/classes/version/7.7/planet.m	(revision 8728)
@@ -0,0 +1,48 @@
+%PLANET class definition
+%
+%   Usage:
+%      md = planet(varargin)
+
+classdef planet < model
+    properties (SetAccess=public) %Planet fields
+		 % {{{1
+		 %Planet specific fields
+		 r=NaN;
+		 theta=NaN;
+		 phi=NaN;
+		 %}}}
+	 end
+	 methods
+		function md=planetmesh(md,varargin) %{{{1
+		%PLANETMESH: build 2d shell mesh
+		%
+		% Usage: md=planetmesh(md,'method','mixed','radius',6378000,'angleresol',1);
+		%        md=planetmesh(md,'method','tria','shape','iso','radius',6378000,'refinement',5);
+		%
+
+		%recover options
+		options=pairoptions(varargin{:});
+
+		method=getfieldvalue(options,'method','mixed');
+		
+		if strcmpi(method,'mixed'),
+			%recover radius and angleresol: 
+			radius=getfieldvalue(options,'radius',6378000); %earth radius as default
+			angleresol=getfieldvalue(options,'angleresol',10); %10 degree resolution
+			
+			%call mixed mesh 
+			md=planetmixedmesh(md,radius,angleresol);
+		else
+			%recover radius, shape and level of refinmenet
+			radius=getfieldvalue(options,'radius',6378000); %earth radius as default
+			refinement=getfieldvalue(options,'refinement',5); %refine 5 times
+			shape=getfieldvalue(options,'shape','ico'); 
+			
+			%call triangular mesh
+			md=planettrimesh(md,shape,radius,refinement);
+		end
+
+		end
+		%}}}
+	 end
+ end
Index: /issm/trunk/src/m/model/plot/applyoptions.m
===================================================================
--- /issm/trunk/src/m/model/plot/applyoptions.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/applyoptions.m	(revision 8728)
@@ -222,5 +222,5 @@
 	if exist(options,'colorbartitle'),
 		set(get(c,'title'),'FontSize',getfieldvalue(options,'colorbarfontsize',fontsize),'String',getfieldvalue(options,'colorbartitle'),...
-			'Color',getfieldvalue(options,'FontColor','k'),'Interpreter',getfieldvalue(options,'Interpreter','none'));
+			'Color',getfieldvalue(options,'FontColor','k'),'Interpreter',getfieldvalue(options,'Interpreter','none'),'Interpreter','Tex');
 	end
 	if exist(options,'colorbarYLabel'),
Index: /issm/trunk/src/m/model/plot/plot_mesh.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_mesh.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/plot_mesh.m	(revision 8728)
@@ -22,5 +22,5 @@
 	patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor',edgecolor,'linewidth',linewidth);
 else
-	if ~md.planet,
+	if ~isplanet,
 		A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); E=elements(:,5); F=elements(:,6);
 		patch( 'Faces', [A B C],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor',edgecolor,'linewidth',linewidth);
@@ -30,7 +30,7 @@
 		patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor',edgecolor,'linewidth',linewidth);
 	else
-		A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); 
+		A=elements(:,1); B=elements(:,2); C=elements(:,3); 
+		if (size(elements,2)==4), D=elements(:,4); else D=C; end
 		patch( 'Faces', [A B C D],  'Vertices', [x y z],'FaceVertexCData', [1 1 1],'FaceColor','none','EdgeColor',edgecolor,'linewidth',linewidth);
-
 	end
 end
Index: /issm/trunk/src/m/model/plot/plot_qmuhistnorm.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_qmuhistnorm.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/plot_qmuhistnorm.m	(revision 8728)
@@ -38,4 +38,6 @@
 if exist(options,'cdfleg'), cdfleg=getfieldvalue(options,'cdfleg'); qmuoptions=[qmuoptions ',''cdfleg'',''' cdfleg '''']; end
 if exist(options,'nrmplt'), nrmplt=getfieldvalue(options,'nrmplt'); qmuoptions=[qmuoptions ',''nrmplt'',''' nrmplt '''']; end
+if exist(options,'EdgeColor'), EdgeColor=getfieldvalue(options,'EdgeColor'); qmuoptions=[qmuoptions ',''EdgeColor'',''' EdgeColor '''']; end
+if exist(options,'FaceColor'), FaceColor=getfieldvalue(options,'FaceColor'); qmuoptions=[qmuoptions ',''FaceColor'',''' FaceColor '''']; end
 
 %use qmu plot
Index: /issm/trunk/src/m/model/plot/plot_section.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_section.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/plot_section.m	(revision 8728)
@@ -25,6 +25,6 @@
 %Get number of curves and generate random colors
 numcurves=size(data,2);
-colorm=getfieldvalue(options,'colormap','jet');
-color=eval([ colorm '(numcurves);']);
+colorm=getfieldvalue(options,'colormap','lines');
+color=eval([ colorm '(numcurves);'])
 options=removefield(options,'colormap',0); %back to default colormap
 
@@ -35,79 +35,39 @@
 end
 
-%Loop over number of curves
-for i=1:numcurves,
-
-	[datai datatype]=processdata(md3d,data(:,i),options);
-
-	%resolution
-	if exist(options,'resolution'),
-		resolution=getfieldvalue(options,'resolution');
-	else %Default resolution
-		resolution=[1000 10*md.numlayers];
-		disp(['plot_section warning: no resolution specified, use default resolution: [horizontal_resolution vertical_resolution]=[' num2str(resolution)  ']']);
-	end
-
-	%Compute section value
-	[elements,x,y,z,s,data_s]=SectionValues(md,datai,getfieldvalue(options,'sectionvalue'),resolution);
-
-	%units
-	if exist(options,'unit'),
-		unit=getfieldvalue(options,'unit');
-		x=x*unit;
-		y=y*unit;
-		z=z*unit;
-		s=s*unit;
-	end
-
-	%2D
-	if is2d
-
-		%Show Section if requested by user
-		if exist(options,'showsection')
-
-			%compute number of labels
-			numlabels=min(getfieldvalue(options,'showsection'),length(s));
-			shift=fix(length(s)/numlabels);
-
-			%plot labels on current graph
-			hold on
-			text(s(1),data_s(1),'1','backgroundcolor',[0.8 0.9 0.8])
-			for i=2:numlabels-1
-				text(s(1+(i-1)*shift),data_s(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
-			end
-			text(s(end),data_s(end),'end','backgroundcolor',[0.8 0.9 0.8])
-
-			%plot section only with labels
-			subplot(nlines,ncols,index2)
-			plot_unit(x_m,y_m,z_m,elements_m,data(:,i),is2d,isplanet,datatype,options)
-			hold on
-			text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
-			for i=2:numlabels-1
-				text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
-			end
-			text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
-			plot(x,y,'-r')
-			axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
-			view(2)
+%read contours: 
+profiles=expread(getfieldvalue(options,'sectionvalue'),1);
+numprofiles=length(profiles);
+
+%Loop over number of profiles: 
+for profile_i=1:numprofiles,
+	profile=profiles(profile_i);
+
+	%Loop over number of curves
+	for i=1:numcurves,
+
+		[datai datatype]=processdata(md3d,data(:,i),options);
+
+		%resolution
+		if exist(options,'resolution'),
+			resolution=getfieldvalue(options,'resolution');
+		else %Default resolution
+			resolution=[1000 10*md.numlayers];
+			disp(['plot_section warning: no resolution specified, use default resolution: [horizontal_resolution vertical_resolution]=[' num2str(resolution)  ']']);
 		end
 
-		%plot section value
-		hold on;
-		subplot(nlines,ncols,index1)
-		%subplot(1,3,[2 3])
-		if i==1,
-			plot(s,data_s,'-','color','k','LineWidth',getfieldvalue(options,'linewidth',1))
-		elseif i==2,
-			plot(s,data_s,':','color',color(i,:),'LineWidth',getfieldvalue(options,'linewidth',1))
-		else
-			plot(s,data_s,'--','color',color(i,:),'LineWidth',getfieldvalue(options,'linewidth',1))
+		%Compute section value
+		[elements,x,y,z,s,data_s]=SectionValues(md,datai,profile,resolution);
+
+		%units
+		if exist(options,'unit'),
+			unit=getfieldvalue(options,'unit');
+			x=x*unit;
+			y=y*unit;
+			z=z*unit;
+			s=s*unit;
 		end
 
-
-		%3D
-	else
-		%plot section value
-		%if user requested view2: 2d plot with curvilinear coordinate
-		if (getfieldvalue(options,'view',3)==2 )
+		%2D
+		if is2d
 
 			%Show Section if requested by user
@@ -120,13 +80,13 @@
 				%plot labels on current graph
 				hold on
-				text(s(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
+				text(s(1),data_s(1),'1','backgroundcolor',[0.8 0.9 0.8])
 				for i=2:numlabels-1
-					text(s(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
-				end
-				text(s(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
+					text(s(1+(i-1)*shift),data_s(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
+				end
+				text(s(end),data_s(end),'end','backgroundcolor',[0.8 0.9 0.8])
 
 				%plot section only with labels
 				subplot(nlines,ncols,index2)
-				plot_unit(x_m,y_m,z_m,elements_m,data(:,i),is2d,datatype,options)
+				plot_unit(x_m,y_m,z_m,elements_m,data(:,i),is2d,isplanet,datatype,options)
 				hold on
 				text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
@@ -140,44 +100,87 @@
 			end
 
+			%plot section value
+			hold on;
 			subplot(nlines,ncols,index1)
-			A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4); 
-			patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
-
+			%subplot(1,3,[2 3])
+			plot(s,data_s,'color',color(i,:),'LineWidth',getfieldvalue(options,'linewidth',1))
+
+
+			%3D
 		else
-
-			%Show Section if requested by user
-			if exist(options,'showsection')
-
-				%compute number of labels
-				numlabels=min(getfieldvalue(options,'showsection'),length(s));
-				shift=fix(length(x)/numlabels);
-
-				%plot labels on current graph
-				hold on
-				text(x(1),y(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
-				for i=2:numlabels-1
-					text(x(1+(i-1)*shift),y(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
-				end
-				text(x(end),y(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
-
-				%plot section only with labels
-				subplot(nlines,ncols,index2)
-				plot_unit(x_m,y_m,z_m,elements_m,data,is2d,datatype,options)
-				hold on
-				text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
-				for i=2:numlabels-1
-					text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
-				end
-				text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
-				plot(x,y,'-r')
-				axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
-				view(2)
+			%plot section value
+			%if user requested view2: 2d plot with curvilinear coordinate
+			if (getfieldvalue(options,'view',3)==2 )
+
+				%Show Section if requested by user
+				if exist(options,'showsection')
+
+					%compute number of labels
+					numlabels=min(getfieldvalue(options,'showsection'),length(s));
+					shift=fix(length(s)/numlabels);
+
+					%plot labels on current graph
+					hold on
+					text(s(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
+					for i=2:numlabels-1
+						text(s(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
+					end
+					text(s(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
+
+					%plot section only with labels
+					subplot(nlines,ncols,index2)
+					plot_unit(x_m,y_m,z_m,elements_m,data(:,i),is2d,datatype,options)
+					hold on
+					text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
+					for i=2:numlabels-1
+						text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
+					end
+					text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
+					plot(x,y,'-r')
+					axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
+					view(2)
+				end
+
+				subplot(nlines,ncols,index1)
+				A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4); 
+				patch( 'Faces', [A B C D], 'Vertices', [s z zeros(length(s),1)],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
+
+			else
+
+				%Show Section if requested by user
+				if exist(options,'showsection')
+
+					%compute number of labels
+					numlabels=min(getfieldvalue(options,'showsection'),length(s));
+					shift=fix(length(x)/numlabels);
+
+					%plot labels on current graph
+					hold on
+					text(x(1),y(1),z(1),'1','backgroundcolor',[0.8 0.9 0.8])
+					for i=2:numlabels-1
+						text(x(1+(i-1)*shift),y(1+(i-1)*shift),z(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
+					end
+					text(x(end),y(end),z(end),'end','backgroundcolor',[0.8 0.9 0.8])
+
+					%plot section only with labels
+					subplot(nlines,ncols,index2)
+					plot_unit(x_m,y_m,z_m,elements_m,data,is2d,datatype,options)
+					hold on
+					text(x(1),y(1),'1','backgroundcolor',[0.8 0.9 0.8])
+					for i=2:numlabels-1
+						text(x(1+(i-1)*shift),y(1+(i-1)*shift),num2str(i),'backgroundcolor',[0.8 0.9 0.8])
+					end
+					text(x(end),y(end),'end','backgroundcolor',[0.8 0.9 0.8])
+					plot(x,y,'-r')
+					axis([min(md.x)-1 max(md.x)+1 min(md.y)-1 max(md.y)+1])
+					view(2)
+				end
+
+				subplot(nlines,ncols,index1)
+				A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4); 
+				patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
+				view(3)
+
 			end
-
-			subplot(nlines,ncols,index1)
-			A=elements(:,1); B=elements(:,2); C=elements(:,3);  D=elements(:,4); 
-			patch( 'Faces', [A B C D], 'Vertices', [x y z],'FaceVertexCData',data_s,'FaceColor','interp','EdgeColor','none');
-			view(3)
-
 		end
 	end
Index: /issm/trunk/src/m/model/plot/plot_unit.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_unit.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/plot_unit.m	(revision 8728)
@@ -42,5 +42,6 @@
 				patch( 'Faces', [C A D F],'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
 			else
-				A=elements(:,1); B=elements(:,2); C=elements(:,3); D=elements(:,4); 
+				A=elements(:,1); B=elements(:,2); C=elements(:,3); 
+				if size(elements,2)==4, D=elements(:,4); else D=C; end
 				patch( 'Faces', [A B C D],  'Vertices', [x y z],'FaceVertexCData', data(:),'FaceColor','interp','EdgeColor',edgecolor);
 			end
Index: /issm/trunk/src/m/model/plot/processmesh.m
===================================================================
--- /issm/trunk/src/m/model/plot/processmesh.m	(revision 8727)
+++ /issm/trunk/src/m/model/plot/processmesh.m	(revision 8728)
@@ -92,5 +92,5 @@
 end
 
-if md.planet,
+if isa(md,'planet'),
 	isplanet=1;
 else
Index: /issm/trunk/src/m/planet/mesh_refine_tri4.m
===================================================================
--- /issm/trunk/src/m/planet/mesh_refine_tri4.m	(revision 8728)
+++ /issm/trunk/src/m/planet/mesh_refine_tri4.m	(revision 8728)
@@ -0,0 +1,146 @@
+function [ FV ] = mesh_refine_tri4(FV)
+
+% mesh_refine_tri4 - creates 4 triangle from each triangle of a mesh
+%
+% [ FV ] = mesh_refine_tri4( FV )
+%
+% FV.vertices   - mesh vertices (Nx3 matrix)
+% FV.faces      - faces with indices into 3 rows
+%                 of FV.vertices (Mx3 matrix)
+% 
+% For each face, 3 new vertices are created at the 
+% triangle edge midpoints.  Each face is divided into 4
+% faces and returned in FV.
+%
+%        B
+%       /\
+%      /  \
+%    a/____\b       Construct new triangles
+%    /\    /\       [A,a,c]
+%   /  \  /  \      [a,B,b]
+%  /____\/____\     [c,b,C]
+% A	     c	   C    [a,b,c]
+% 
+% It is assumed that the vertices are listed in clockwise order in
+% FV.faces (A,B,C above), as viewed from the outside in a RHS coordinate
+% system.
+% 
+% See also: mesh_refine, sphere_tri, sphere_project
+% 
+
+
+% ---this method is not implemented, but the idea here remains...
+% This can be done until some minimal distance (D) of the mean 
+% distance between vertices of all triangles is achieved.  If
+% no D argument is given, the function refines the mesh once.
+% Alternatively, it could be done until some minimum mean 
+% area of faces is achieved.  As is, it just refines once.
+
+
+% $Revision: 1.1 $ $Date: 2004/11/12 01:32:35 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% History:  05/2002, Darren.Weber_at_radiology.ucsf.edu, created
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+tic;
+fprintf('...refining mesh (tri4)...')
+
+% NOTE
+% The centroid is located one third of the way from each vertex to 
+% the midpoint of the opposite side. Each median divides the triangle 
+% into two equal areas; all the medians together divide it into six 
+% equal parts, and the lines from the median point to the vertices 
+% divide the whole into three equivalent triangles.
+
+% Each input triangle with vertices labelled [A,B,C] as shown
+% below will be turned into four new triangles:
+%
+% Make new midpoints
+% a = (A+B)/2
+% b = (B+C)/2
+% c = (C+A)/2
+%
+%        B
+%       /\
+%      /  \
+%    a/____\b       Construct new triangles
+%    /\    /\       [A,a,c]
+%   /  \  /  \      [a,B,b]
+%  /____\/____\     [c,b,C]
+% A	     c	   C    [a,b,c]
+%
+
+% Initialise a new vertices and faces matrix
+Nvert = size(FV.vertices,1);
+Nface = size(FV.faces,1);
+V2 = zeros(Nface*3,3);
+F2 = zeros(Nface*4,3);
+
+for f = 1:Nface,
+    
+    % Get the triangle vertex indices
+    NA = FV.faces(f,1);
+    NB = FV.faces(f,2);
+    NC = FV.faces(f,3);
+    
+    % Get the triangle vertex coordinates
+    A = FV.vertices(NA,:);
+    B = FV.vertices(NB,:);
+    C = FV.vertices(NC,:);
+    
+    % Now find the midpoints between vertices
+    a = (A + B) ./ 2;
+    b = (B + C) ./ 2;
+    c = (C + A) ./ 2;
+    
+    % Find the length of each median
+    %A2blen = sqrt ( sum( (A - b).^2, 2 ) );
+    %B2clen = sqrt ( sum( (B - c).^2, 2 ) );
+    %C2alen = sqrt ( sum( (C - a).^2, 2 ) );
+    
+    % Store the midpoint vertices, while
+    % checking if midpoint vertex already exists
+    [FV, Na] = mesh_find_vertex(FV,a);
+    [FV, Nb] = mesh_find_vertex(FV,b);
+    [FV, Nc] = mesh_find_vertex(FV,c);
+    
+    % Create new faces with orig vertices plus midpoints
+    F2(f*4-3,:) = [ NA, Na, Nc ];
+    F2(f*4-2,:) = [ Na, NB, Nb ];
+    F2(f*4-1,:) = [ Nc, Nb, NC ];
+    F2(f*4-0,:) = [ Na, Nb, Nc ];
+    
+end
+
+% Replace the faces matrix
+FV.faces = F2;
+
+t=toc; fprintf('done (%5.2f sec)\n',t);
+
+return
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [FV, N] = mesh_find_vertex(FV,vertex)
+
+    Vn = size(FV.vertices,1);
+    Va = repmat(vertex,Vn,1);
+    Vexist = find( FV.vertices(:,1) == Va(:,1) & ...
+                   FV.vertices(:,2) == Va(:,2) & ...
+                   FV.vertices(:,3) == Va(:,3) );
+    if Vexist,
+        if size(Vexist) == [1,1],
+            N = Vexist;
+        else,
+            msg = sprintf('replicated vertices');
+            error(msg);
+        end
+    else
+        FV.vertices(end+1,:) = vertex;
+        N = size(FV.vertices,1);
+    end
+
+return
Index: /issm/trunk/src/m/planet/planetmixedmesh.m
===================================================================
--- /issm/trunk/src/m/planet/planetmixedmesh.m	(revision 8728)
+++ /issm/trunk/src/m/planet/planetmixedmesh.m	(revision 8728)
@@ -0,0 +1,99 @@
+function md=planetmixedmesh(md,radius,angleresol)
+%PLANETMIXEDMESH: build 2d shell mesh
+%
+% Usage: md=planetmixedshell(md,radius,angleresol)
+%
+
+conv=pi/180;
+
+r=radius;
+theta=(0:angleresol:360)';
+theta=theta*conv;
+phi=(0:angleresol:180)';
+phi=phi(2:end-1)*conv;
+
+nphi=length(phi);
+nthetha=length(theta);
+
+x=zeros(nphi*nthetha,1);
+y=zeros(nphi*nthetha,1);
+z=zeros(nphi*nthetha,1);
+
+for i=1:nphi,
+	phii=phi(i);
+	x((i-1)*nthetha+1:i*nthetha)=r.*cos(theta).*sin(phii);
+	y((i-1)*nthetha+1:i*nthetha)=r.*sin(theta).*sin(phii);
+	z((i-1)*nthetha+1:i*nthetha)=r.*cos(phii);
+end
+
+quads=zeros((nthetha-1)*(nphi-1),4);
+
+%build elements
+for i=1:nphi-1,
+	for j=1:nthetha-1,
+		count=(i-1)*(nthetha-1)+j;
+	
+		A=(i-1)*nthetha+j;
+		B=(i-1)*nthetha+j+1;
+		C=(i)*nthetha+j+1;
+		D=(i)*nthetha+j;
+		
+		quads(count,:)=[A B C D];
+	end
+end
+
+%now, add polar elements:
+%first north: phi = 0;
+x=[r.*cos(0).*sin(0);x];
+y=[r.*sin(0).*sin(0);y];
+z=[r.*cos(0);z];
+
+%add 1 to quads: 
+quads=quads+1;
+
+%add tria elements:
+trias=zeros(nthetha-1,4);
+
+for i=1:nthetha-1,
+	A=1;
+	B=i+1;
+	C=i+2;
+	trias(i,:)=[A B C NaN];
+end
+
+quads=[trias;quads];
+
+
+% now add south pole: 
+phii=180*conv;
+x=[x;r.*cos(phii).*sin(phii)];
+y=[y;r.*sin(phii).*sin(phii)];
+z=[z;r.*cos(phii)];
+nods=length(x);
+
+%add tria elements:
+trias=zeros(nthetha-1,4);
+
+start=nods-nthetha;
+for i=1:nthetha-1,
+	A=start+i-1;
+	B=start+i;
+	C=nods;
+	trias(i,:)=[A B C NaN];
+end
+
+quads=[quads;trias];
+
+
+md.elements=quads;
+md.x=x;
+md.y=y;
+md.z=z;
+md.r=sqrt(x.^2+y.^2+z.^2);
+md.theta=acos(z./r);
+md.phi=atan2(y,x);
+
+md.numberofnodes=length(md.x);
+md.numberofelements=size(md.elements,1);
+
+md.dim=3;
Index: /issm/trunk/src/m/planet/planettrimesh.m
===================================================================
--- /issm/trunk/src/m/planet/planettrimesh.m	(revision 8728)
+++ /issm/trunk/src/m/planet/planettrimesh.m	(revision 8728)
@@ -0,0 +1,20 @@
+function md=planettrimesh(md,shape,radius,refinement)
+%PLANETTRIMESH: build 2d shell mesh
+%
+% Usage: md=planettrimesh(md,shape,radius,refinement)
+%
+
+results = sphere_tri(shape,refinement,radius);
+md.x=results.vertices(:,1);
+md.y=results.vertices(:,2);
+md.z=results.vertices(:,3);
+md.elements=results.faces;
+
+md.r=sqrt(md.x.^2+md.y.^2+md.z.^2);
+md.theta=acos(md.z./md.r);
+md.phi=atan2(md.y,md.x);
+
+md.numberofnodes=length(md.x);
+md.numberofelements=size(md.elements,1);
+
+md.dim=3;
Index: /issm/trunk/src/m/planet/runme.m
===================================================================
--- /issm/trunk/src/m/planet/runme.m	(revision 8728)
+++ /issm/trunk/src/m/planet/runme.m	(revision 8728)
@@ -0,0 +1,4 @@
+% 5 -> level of refinment
+% 1000 -> radius
+FV = sphere_tri('ico',5,1000); 
+patch('vertices',FV.vertices,'faces',FV.faces,'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]);
Index: /issm/trunk/src/m/planet/sphere_project.m
===================================================================
--- /issm/trunk/src/m/planet/sphere_project.m	(revision 8728)
+++ /issm/trunk/src/m/planet/sphere_project.m	(revision 8728)
@@ -0,0 +1,66 @@
+function V = sphere_project(v,r,c)
+
+% sphere_project - project point X,Y,Z to the surface of sphere radius r
+% 
+% V = sphere_project(v,r,c)
+% 
+% Cartesian inputs:
+% v is the vertex matrix, Nx3 (XYZ)
+% r is the sphere radius, 1x1 (default 1)
+% c is the sphere centroid, 1x3 (default 0,0,0)
+%
+% XYZ are converted to spherical coordinates and their radius is
+% adjusted according to r, from c toward XYZ (defined with theta,phi)
+% 
+% V is returned as Cartesian 3D coordinates
+% 
+
+% $Revision: 1.1 $ $Date: 2004/11/12 01:32:36 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% History:  06/2002, Darren.Weber_at_radiology.ucsf.edu, created
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if ~exist('v','var'),
+    msg = sprintf('SPHERE_PROJECT: No input vertices (X,Y,Z)\n');
+    error(msg);
+end
+
+X = v(:,1);
+Y = v(:,2);
+Z = v(:,3);
+
+if ~exist('c','var'),
+    xo = 0;
+    yo = 0;
+    zo = 0;
+else
+    xo = c(1);
+    yo = c(2);
+    zo = c(3);
+end
+
+if ~exist('r','var'), r = 1; end
+
+% alternate method is to use unit vector of V
+% [ n = 'magnitude(V)'; unitV = V ./ n; ]
+% to change the radius, multiply the unitV
+% by the radius required.  This avoids the
+% use of arctan functions, which have branches.
+
+
+% Convert Cartesian X,Y,Z to spherical (radians)
+theta = atan2( (Y-yo), (X-xo) );
+phi   = atan2( sqrt( (X-xo).^2 + (Y-yo).^2 ), (Z-zo) );
+% do not calc: r = sqrt( (X-xo).^2 + (Y-yo).^2 + (Z-zo).^2);
+
+%   Recalculate X,Y,Z for constant r, given theta & phi.
+R = ones(size(phi)) * r;
+x = R .* sin(phi) .* cos(theta);
+y = R .* sin(phi) .* sin(theta);
+z = R .* cos(phi);
+
+V = [x y z];
+
+return
Index: /issm/trunk/src/m/planet/sphere_tri.m
===================================================================
--- /issm/trunk/src/m/planet/sphere_tri.m	(revision 8728)
+++ /issm/trunk/src/m/planet/sphere_tri.m	(revision 8728)
@@ -0,0 +1,204 @@
+function [FV] = sphere_tri(shape,maxlevel,r,winding)
+
+% sphere_tri - generate a triangle mesh approximating a sphere
+% 
+% Usage: FV = sphere_tri(shape,Nrecurse,r,winding)
+% 
+%   shape is a string, either of the following:
+%   'ico'   starts with icosahedron (most even, default)
+%   'oct'   starts with octahedron
+%   'tetra' starts with tetrahedron (least even)
+%
+%   Nrecurse is int >= 0, setting the recursions (default 0)
+%
+%   r is the radius of the sphere (default 1)
+%
+%   winding is 0 for clockwise, 1 for counterclockwise (default 0).  The
+%   matlab patch command gives outward surface normals for clockwise
+%   order of vertices in the faces (viewed from outside the surface).
+%
+%   FV has fields FV.vertices and FV.faces.  The vertices 
+%   are listed in clockwise order in FV.faces, as viewed 
+%   from the outside in a RHS coordinate system.
+% 
+% The function uses recursive subdivision.  The first
+% approximation is an platonic solid, either an  icosahedron,
+% octahedron or a tetrahedron.  Each level of refinement 
+% subdivides each triangle face by a factor of 4 (see also 
+% mesh_refine).  At each refinement, the vertices are 
+% projected to the sphere surface (see sphere_project).
+% 
+% A recursion level of 3 or 4 is a good sphere surface, if
+% gouraud shading is used for rendering.
+% 
+% The returned struct can be used in the patch command, eg:
+% 
+% % create and plot, vertices: [2562x3] and faces: [5120x3]
+% FV = sphere_tri('ico',4,1);
+% lighting phong; shading interp; figure;
+% patch('vertices',FV.vertices,'faces',FV.faces,...
+%       'facecolor',[1 0 0],'edgecolor',[.2 .2 .6]);
+% axis off; camlight infinite; camproj('perspective');
+% 
+% See also: mesh_refine, sphere_project
+%
+
+
+
+% $Revision: 1.2 $ $Date: 2005/07/20 23:07:03 $
+
+% Licence:  GNU GPL, no implied or express warranties
+% Jon Leech (leech @ cs.unc.edu) 3/24/89
+% icosahedral code added by Jim Buddenhagen (jb1556@daditz.sbc.com) 5/93
+% 06/2002, adapted from c to matlab by Darren.Weber_at_radiology.ucsf.edu
+% 05/2004, reorder of the faces for the 'ico' surface so they are indeed
+% clockwise!  Now the surface normals are directed outward.  Also reset the
+% default recursions to zero, so we can get out just the platonic solids.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+eegversion = '$Revision: 1.2 $';
+fprintf('SPHERE_TRI [v %s]\n',eegversion(11:15)); tic
+
+if ~exist('shape','var') || isempty(shape),
+    shape = 'ico';
+end
+fprintf('...creating sphere tesselation based on %s\n',shape);
+
+% default maximum subdivision level
+if ~exist('maxlevel','var') || isempty(maxlevel) || maxlevel < 0,
+    maxlevel = 0;
+end
+
+% default radius
+if ~exist('r','var') || isempty(r),
+    r = 1;
+end
+
+if ~exist('winding','var') || isempty(winding),
+    winding = 0;
+end
+
+
+% -----------------
+% define the starting shapes
+
+shape = lower(shape);
+
+switch shape,
+case 'tetra',
+    
+    % Vertices of a tetrahedron
+    sqrt_3 = 0.5773502692;
+    
+    tetra.v = [  sqrt_3,  sqrt_3,  sqrt_3 ;   % +X, +Y, +Z  - PPP
+                -sqrt_3, -sqrt_3,  sqrt_3 ;   % -X, -Y, +Z  - MMP
+                -sqrt_3,  sqrt_3, -sqrt_3 ;   % -X, +Y, -Z  - MPM
+                 sqrt_3, -sqrt_3, -sqrt_3 ];  % +X, -Y, -Z  - PMM
+	
+    % Structure describing a tetrahedron
+    tetra.f = [ 1, 2, 3;
+                1, 4, 2;
+                3, 2, 4;
+                4, 1, 3 ];
+    
+    FV.vertices = tetra.v;
+    FV.faces    = tetra.f;
+    
+case 'oct',
+    
+    % Six equidistant points lying on the unit sphere
+    oct.v = [  1,  0,  0 ;  %  X
+              -1,  0,  0 ; 	% -X
+               0,  1,  0 ;  %  Y
+               0, -1,  0 ; 	% -Y
+               0,  0,  1 ; 	%  Z
+               0,  0, -1 ];	% -Z
+	
+    % Join vertices to create a unit octahedron
+    oct.f = [ 1 5 3 ;    %  X  Z  Y  -  First the top half
+              3 5 2 ;    %  Y  Z -X
+              2 5 4 ;    % -X  Z -Y
+              4 5 1 ;    % -Y  Z  X
+              1 3 6 ;    %  X  Y -Z  -  Now the bottom half
+              3 2 6 ;    %  Y  Z -Z
+              2 4 6 ;    % -X  Z -Z
+              4 1 6 ];   % -Y  Z -Z
+    
+    FV.vertices = oct.v;
+    FV.faces    = oct.f;
+    
+case 'ico',
+    
+    % Twelve vertices of icosahedron on unit sphere
+    tau = 0.8506508084; % t=(1+sqrt(5))/2, tau=t/sqrt(1+t^2)
+    one = 0.5257311121; % one=1/sqrt(1+t^2) , unit sphere
+    
+    ico.v( 1,:) = [  tau,  one,    0 ]; % ZA
+    ico.v( 2,:) = [ -tau,  one,    0 ]; % ZB
+    ico.v( 3,:) = [ -tau, -one,    0 ]; % ZC
+    ico.v( 4,:) = [  tau, -one,    0 ]; % ZD
+    ico.v( 5,:) = [  one,   0 ,  tau ]; % YA
+    ico.v( 6,:) = [  one,   0 , -tau ]; % YB
+    ico.v( 7,:) = [ -one,   0 , -tau ]; % YC
+    ico.v( 8,:) = [ -one,   0 ,  tau ]; % YD
+    ico.v( 9,:) = [   0 ,  tau,  one ]; % XA
+    ico.v(10,:) = [   0 , -tau,  one ]; % XB
+    ico.v(11,:) = [   0 , -tau, -one ]; % XC
+    ico.v(12,:) = [   0 ,  tau, -one ]; % XD
+    
+    % Structure for unit icosahedron
+    ico.f = [  5,  8,  9 ;
+               5, 10,  8 ;
+               6, 12,  7 ;
+               6,  7, 11 ;
+               1,  4,  5 ;
+               1,  6,  4 ;
+               3,  2,  8 ;
+               3,  7,  2 ;
+               9, 12,  1 ;
+               9,  2, 12 ;
+              10,  4, 11 ;
+              10, 11,  3 ;
+               9,  1,  5 ;
+              12,  6,  1 ;
+               5,  4, 10 ;
+               6, 11,  4 ;
+               8,  2,  9 ;
+               7, 12,  2 ;
+               8, 10,  3 ;
+               7,  3, 11 ];
+	
+    FV.vertices = ico.v;
+    FV.faces    = ico.f;
+end
+
+
+% -----------------
+% refine the starting shapes with subdivisions
+if maxlevel,
+    
+    % Subdivide each starting triangle (maxlevel) times
+    for level = 1:maxlevel,
+        
+        % Subdivide each triangle and normalize the new points thus
+        % generated to lie on the surface of a sphere radius r.
+        FV = mesh_refine_tri4(FV);
+        FV.vertices = sphere_project(FV.vertices,r);
+        
+        % An alternative might be to define a min distance
+        % between vertices and recurse or use fminsearch
+        
+    end
+end
+
+if winding,
+    fprintf('...returning counterclockwise vertex order (viewed from outside)\n');
+    FV.faces = FV.faces(:,[1 3 2]);
+else
+    fprintf('...returning clockwise vertex order (viewed from outside)\n');
+end
+
+t=toc; fprintf('...done (%6.2f sec)\n\n',t);
+
+return