Index: /issm/trunk-jpl/src/m/contrib/adhikari/greens.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/greens.m	(revision 28265)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/greens.m	(revision 28265)
@@ -0,0 +1,59 @@
+function [vert, horz, accl] = greens(dist,love_h,love_l,love_k); 
+% compute greens functions for vertical & horizontal crustal motion, acceleration and tilt 
+% see Ferrell, 1972, equations 40, 45, 46-49
+% 
+% usage: 
+%		[vert, horz, accl] = greens(dist,love_h,love_l,love_k); 
+%		
+%		dist = grid points along the distance away from the disc center [m] 
+%		vert = Green's function for vertical crustal motion 
+%		horz = Green's function for horizontal crustal motion
+%		accl = Green's function for acceleration 
+%		love_h = load Love numbers h for vertical crustal motion. 
+%		love_l = load Love numbers l for horizontal crustal motion. 
+%		love_k = load Love numbers k for gravitational potential 
+%		
+
+	% compute P(x), dP(x)/dx, d2P(x)/dx2
+	%---------------------------------------------------------------------
+	theta=km2deg(dist/1000)';
+	ang = theta/180*pi; 
+	alpha=cos(ang);
+	m=length(alpha);
+	n=length(love_h)-1; % number of degrees. 
+	p_value = p_polynomial_value(m,n,alpha);
+	p_prime = p_polynomial_prime(m,n,alpha);
+	%---------------------------------------------------------------------
+	
+	love_h_inf = love_h(end); 
+	love_l_inf = love_l(end)*n; 
+	love_k_inf = love_k(end)*n; 
+	
+	for jj=1:n+1
+		n_d = jj-1; % n degree 
+		if n_d==0
+			coeff_horz(jj) = love_l(jj) - love_l_inf/(n_d+1e-12); 
+		else
+			coeff_horz(jj) = love_l(jj) - love_l_inf/n_d; 
+		end
+		coeff_acc(jj) = n_d + 2*love_h(jj) - (n_d + 1)*love_k(jj); 
+		coeff_love_h(jj) = love_h(jj) - love_h_inf; 
+		coeff_love_k(jj) = n_d*love_k(jj) - love_k_inf; 
+	end
+
+	% vertical. 
+	vert = 0.5*love_h_inf./sin(ang/2) + sum(bsxfun(@times,p_value,coeff_love_h),2); 
+	%vert = sum(bsxfun(@times,p_value,love_h'),2); 
+
+	% horizontal. 
+	horz = -cos(ang/2).*(1 + 2*sin(ang/2)) ./ (2*sin(ang/2) .* (1 + sin(ang/2))) * love_l_inf +...
+		+ sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,coeff_horz)),2); 
+	%horz = sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,love_l')),2); 
+
+	% acceleration. 
+	accl = -0.25./sin(ang/2) +...
+		+ love_h_inf./sin(ang/2) + 2*sum(bsxfun(@times,p_value,coeff_love_h),2) +...
+		- love_k_inf*0.5./sin(ang/2) - sum(bsxfun(@times,p_value,coeff_love_k),2) +...
+		- sum(bsxfun(@times,p_value,love_k'),2); % last term is stable by itself!  
+	%accl = sum(bsxfun(@times,p_value,coeff_acc),2); 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/wahr.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/wahr.m	(revision 28264)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/wahr.m	(revision 28265)
@@ -1,4 +1,4 @@
-function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); 
-% a function to compute vertical and horizontal crustal motion... 
+function [vert, horz, accl] = wahr(disc_rad,xi,love_h,love_l,love_k); 
+% a function to compute 3D crustal motion and change in gravitational acceleration... 
 % ...for a disc load, inspired by Figure 1 of Wahr et al., 2013. 
 % 
@@ -8,8 +8,10 @@
 %		vert = vertical crustal motion [m] 
 %		horz = horizontal crustal motion [m] 
+%		accl = change in gravitational acceleration [m/s2] 
 %		disc_rad = disc radius [m] => set to 20 km to replicate the Wahr experiment. 
 %		xi = grid points along the distance away from the disc center [m] 
-%		love_h = load Love numbers h for the vertical crustal motion. 
-%		love_l = load Love numbers l for the horizontal crustal motion. 
+%		love_h = load Love numbers h for vertical crustal motion. 
+%		love_l = load Love numbers l for horizontal crustal motion. 
+%		love_k = load Love numbers l for gravitational potential. 
 %		
 
@@ -31,28 +33,37 @@
 	% 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)); 
+
+	%for jj=2:n+1
+	%	nnn = jj-1; 
+	%	tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); 
+	%end
+
+	for jj=1:n+1
+		n_d = jj-1; % n degree 
+		coeff(jj) = 1/(2*(jj-1)+1); 
+		coeff_accl(jj) = (n_d + 2*love_h(jj) - (n_d+1)*love_k(jj)) /  (2*n_d + 1); 
+		if jj==1
+			tau(jj) = 0.5*(1-cos(disc_rad)); % tau_0 
+		else
+			tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); 
+		end
 	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); 
 
-	disc=sum(bsxfun(@times,p_value,tau'),2); 
+	g1 = -sum(bsxfun(@times,p_value,tau.*coeff.*love_h'),2); 
+	g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,tau.*coeff.*love_l')),2); 
+	% gravitational potential 
+	accl = -sum(bsxfun(@times,p_value,tau.*coeff_accl),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; 
+	% constants 
+	const = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; % [m] 
+	const_accl = 4*pi*(6.67408*10^-11); % [m/s2]
 
 	% vertical and horizontal solutions
-	vert = g1*coeff; % m 
-	horz = g5*coeff; % m 
+	vert = g1*const; % m g1 and g1_check give the same results. 
+	horz = g5*const; % m
+	accl = accl*const_accl; % [m/s2]