	function [v, gamma] = savpres_depth(x,x2,d1,d2,t,H,tR,T);

%	 [v, gamma] = savpres(x,t,H,tR,T);
%
% Postseismic displacements due to slip on fault of
% depth H, in an elastic layer of thickness H, over
% a Maxwell viscoelastic half-space. Savage (1990)
%
%	tR = Maxwell relaxation time in same units as time
%	T = recurrence interval

v = zeros(length(x2),length(t));
if nargout == 2
	gamma = zeros(length(x2),length(t));	
end

% main loop
gs = makeTs(t/tR, T/tR);

for n=1:15
   if sum(isnan(gs(:,n)))==0
      g=gs(:,n)';
		g = T.*g.*exp(-t)./factorial(n-1);
      
      D=d2;
      Fn2=-(-atan(x./(-D - 2*H*n + x2)) + atan(x./(D - 2*H*n + x2)) - atan(x./(-D + 2*H*n + x2)) + ...
         atan(x./(D + 2*H*n + x2)))/(2*pi);
      D=d1;
      Fn1=-(-atan(x./(-D - 2*H*n + x2)) + atan(x./(D - 2*H*n + x2)) - atan(x./(-D + 2*H*n + x2)) + ...
         atan(x./(D + 2*H*n + x2)))/(2*pi);
      Fn=Fn2-Fn1;
      
      
      v = v + Fn'*g;
      
  		if nargout == 2
           D=d2;
			  Gn2 = -((D + 2*H*n - x2)./(x.^2 + (D + 2*H*n - x2).^2) + (D - 2*H*n + x2)./(x.^2 + (D - 2*H*n + x2).^2) + ... 
  			(D - 2*H*n - x2)./(x.^2 + (-D + 2*H*n + x2).^2) + (D + 2*H*n + x2)./(x.^2 + (D + 2*H*n + x2).^2))/(2*pi);
		   D=d1;
			  Gn1 = -((D + 2*H*n - x2)./(x.^2 + (D + 2*H*n - x2).^2) + (D - 2*H*n + x2)./(x.^2 + (D - 2*H*n + x2).^2) + ... 
  			(D - 2*H*n - x2)./(x.^2 + (-D + 2*H*n + x2).^2) + (D + 2*H*n + x2)./(x.^2 + (D + 2*H*n + x2).^2))/(2*pi);
	       Gn=Gn2-Gn1;
	      gamma = gamma + Gn'*g;	
		end
      
   end
end




