Index: /issm/trunk-jpl/examples/EsaWahr/RoundDomain.exp
===================================================================
--- /issm/trunk-jpl/examples/EsaWahr/RoundDomain.exp	(revision 22793)
+++ /issm/trunk-jpl/examples/EsaWahr/RoundDomain.exp	(revision 22793)
@@ -0,0 +1,49 @@
+## Name:
+## Icon:0
+# Points Count Value
+41 1
+# X pos Y pos
+20000 0
+19754 3128.7
+19021 6180.3
+17820 9079.8
+16180 11756
+14142 14142
+11756 16180
+9079.8 17820
+6180.3 19021
+3128.7 19754
+0 20000
+-3128.7 19754
+-6180.3 19021
+-9079.8 17820
+-11756 16180
+-14142 14142
+-16180 11756
+-17820 9079.8
+-19021 6180.3
+-19754 3128.7
+-20000 0
+-19754 -3128.7
+-19021 -6180.3
+-17820 -9079.8
+-16180 -11756
+-14142 -14142
+-11756 -16180
+-9079.8 -17820
+-6180.3 -19021
+-3128.7 -19754
+-0 -20000
+3128.7 -19754
+6180.3 -19021
+9079.8 -17820
+11756 -16180
+14142 -14142
+16180 -11756
+17820 -9079.8
+19021 -6180.3
+19754 -3128.7
+20000 0
+
+
+
Index: /issm/trunk-jpl/examples/EsaWahr/runme.m
===================================================================
--- /issm/trunk-jpl/examples/EsaWahr/runme.m	(revision 22792)
+++ /issm/trunk-jpl/examples/EsaWahr/runme.m	(revision 22793)
@@ -1,15 +1,38 @@
 
 clear all;
-steps=[4];
+steps=[1:6]; % [1:6]; 
 
-if any(steps==1)
-	disp('   Step 1: Mesh creation');
+if any(steps==0)
+	disp('   Step 0: Mesh creation');
 
-	%Generate initial uniform mesh 
+	% Generate initial uniform mesh 
 	md=roundmesh(model,100000,10000);  % Domain radius and element size (meters) 
 
 	save ./Models/EsaWahr.Mesh md;
 	
+	plotmodel(md,'data','mesh');
+	%export_fig('Fig0.pdf'); 
+
+end
+
+if any(steps==1)
+	disp('   Step 1: Anisotropic mesh creation');
+
+	% Generate initial uniform mesh 
+	md=roundmesh(model,100000,1000);  % Domain radius and element size (meters) 
+	%plotmodel(md,'data','mesh'); 
+
+	% Define a fake field to have anisotropic meshing 
+	disc_radius=20*1000; % km => m 
+	rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2);	% radial distance [m] 
+	field = abs(rad_dist-disc_radius); 
+	%plotmodel(md,'data',field); 
+
+	md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m 
+
+	save ./Models/EsaWahr.Mesh md;
+	
 	plotmodel (md,'data','mesh');
+	%export_fig('Fig1.pdf'); 
 
 end
@@ -20,6 +43,4 @@
 
 	% primer 
-	nv=md.mesh.numberofvertices; 
-	ne=md.mesh.numberofelements; 
 	rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; 
 
@@ -30,5 +51,5 @@
 
 	% Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly  
-	md.esa.deltathickness = zeros(ne,1); 
+	md.esa.deltathickness = zeros(md.mesh.numberofelements,1); 
 	disc_radius=20; % km 
 	rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000;	% radial distance in km 
@@ -38,4 +59,5 @@
 	
 	plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]');
+	%export_fig('Fig2.pdf'); 
 
 end
@@ -67,4 +89,5 @@
 	% Miscellaneous: 
 	md.miscellaneous.name='EsaWahr';
+	%% IGNORE BUT DO NOT DELETE %% 
 	
 	save ./Models/EsaWahr.Parameterization md; 
@@ -86,29 +109,19 @@
 	md=solve(md,'Esa');
 
-	save ./Models/EsaWahr.EsaSolutions md; 
+	save ./Models/EsaWahr.Solution md; 
 
 end
 
-return; 
+if any(steps==5)
+	disp('   Step 5: Plot solutions');
+	md = loadmodel('./Models/EsaWahr.Solution');
 
-% high-res mesh...
-% plot step. 
-% analytical step. 
-
-	
+	% m => mm 
 	vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
-	horz_n = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
-	horz_e = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+	horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 
+	horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 
 	horz = sqrt(horz_n.^2+horz_e.^2); 
 
-	% grid data from the center 
-	xi=[0:100:60000]; % 100 m resolution 
-	yi=zeros(1,length(xi)); 
-	vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
-	horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
-
-	return; 
-
-	% figure S1a {{{ 
+	% plot 
 	set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
 	figure('Position', [100, 100, 800, 600]);
@@ -116,5 +129,5 @@
 		'xTickLabel#all',[],'yTickLabel#all',[],...
 		'colormap#all','jet','colorbarcornerposition#all','south',...
-		'expdisp#all','../Exp/RoundDomain.exp',...
+		'expdisp#all','./RoundDomain.exp',...
 		'gridded#all',1,...
 		'axispos#1',[0.02 0.505 0.47 0.47],...
@@ -132,28 +145,55 @@
 		'axispos',[0.505 0.02 0.47 0.47],...
 		'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); 
-	% }}} 
-	%export_fig('./Fig/Fig_S1a.pdf'); 
 
-	% figure S1b {{{
-	set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
+	%export_fig('Fig5.pdf'); 
+
+end
+
+if any(steps==6)
+	disp('   Step 6: Compare results against Wahr semi-analytic solutions');
+	md = loadmodel('./Models/EsaWahr.Solution');
+
+	% numerical solutions 
+	% m => mm 
+	vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
+	horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 
+	horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 
+	horz = sqrt(horz_n.^2+horz_e.^2); 
+	% grid data from the center 
+	xi=[0:500:100000]; % 500 m resolution 
+	yi=zeros(1,length(xi)); 
+	vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
+	horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
+
+	% semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) 
+	disc_radius = 20*1000; % 20 km 
+	[vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l);
+
+	% plot 
+	set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); 
 	figure1=figure('Position', [100, 100, 700, 400]);
 	ylabel_1={'0',' ','1','','2','','3',''}; 
 	axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],...
-		'XTick',[0:10:60],'xlim',[0 60],...
+		'XTick',[0:10:100],'xlim',[0 100],...
 		'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); 
 		box(axes1,'on'); hold(axes1,'all'); grid on; 
 		xlabel(axes1,'Radial distance [km]'); 
 		ylabel(axes1,'Displacement [mm]');
-		plot([20 20],[0 3.5],'-.k','linewidth',2,'parent',axes1); 
-		h1=plot(xi/1000,vert_track*917/1000,'color',[0.8500 0.3250 0.0980],'linewidth',3,'parent',axes1); 
-		h2=plot(xi/1000,horz_track*917/1000,'b','linewidth',3,'parent',axes1); 
+		plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); 
+		% analytic soln 
+		h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); 
+		h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); 
+		% ISSM 
+		h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); 
+		h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); 
 		% first box 
 		ag1 = gca;
-		leg1a = legend(ag1,[h1,h2],'Vertical (upward)','Horizontal (away from disc)',1);
-		set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',24); 
+		leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)',1);
+		set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); 
 		% 
 	set(gcf,'color','w');
-	% }}} 
-	%export_fig('./Fig/Fig_S1b.pdf'); 
+	
+	%export_fig('Fig6.pdf'); 
 
+end
 
Index: /issm/trunk-jpl/examples/EsaWahr/wahr.m
===================================================================
--- /issm/trunk-jpl/examples/EsaWahr/wahr.m	(revision 22793)
+++ /issm/trunk-jpl/examples/EsaWahr/wahr.m	(revision 22793)
@@ -0,0 +1,46 @@
+% compute semi-analytic solutions for a disc load 
+function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); 
+
+	disc_rad = disc_rad/1000; % km 
+	% compute P(x), dP(x)/dx, d2P(x)/dx2
+	%---------------------------------------------------------------------
+	% compute p_value 
+	theta=km2deg(xi/1000)';
+	ang = theta/180*pi; 
+	alpha=cos(ang);
+	m=length(alpha);
+	n=length(love_h)-1; 
+	p_value = p_polynomial_value(m,n,alpha);
+	p_prime = p_polynomial_prime(m,n,alpha);
+	%---------------------------------------------------------------------
+	nn=[0:n];
+	nn_plus_1=nn+1; 
+
+	% disc radius in degree 
+	disc_rad = km2deg(disc_rad)/180*pi; 
+	tau=zeros(size(love_h)); 
+	tau(1) = 0.5*(1-cos(disc_rad)); % tau_0 
+	p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad));
+	p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad));
+	for jj=2:n+1
+		nnn = jj-1; 
+		tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); 
+	end
+
+	const=zeros(size(love_h)); 
+	for jj=1:n+1
+		const(jj) = 1/(2*(jj-1)+1); 
+	end
+
+	disc=sum(bsxfun(@times,p_value,tau'),2); 
+
+	g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2); 
+	g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2); 
+
+	% coeff 
+	coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; 
+
+	% vertical and horizontal solutions in mm 
+	vert = g1*coeff*1000; % mm 
+	horz = g5*coeff*1000; % mm 
+