	function [v, gamma] = savpres_patch(x,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(x),length(t));
if nargout == 2
	gamma = zeros(length(x),length(t));	
end
x2=0;
% main loop
gs = makeTs(t/tR, T/tR);


	for n=1:15
      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




%% dislocations below
g = ones(size(g));
for k=1:500;
   Fn=1/pi*(atan(x/(2*(n+k)*H-d2))-atan(x/(2*(n+k)*H-d1))+atan(x/(2*(n+k)*H+d1))-atan(x/(2*(n+k)*H+d2)));
%   Gn=-(1/((-d1 + 2*H*(k + n))*(1 + x^2/(-d1 + 2*H*(k + n))^2))) +...
%      1/((d1 + 2*H*(k + n))* ...
%    (1 + x^2/(d1 + 2*H*(k + n))^2)) + ... 
%  1/((-d2 + 2*H*(k + n))* ...
%    (1 + x^2/(-d2 + 2*H*(k + n))^2)) - ... 
%  1/((d2 + 2*H*(k + n))*(1 + x^2/(d2 + 2*H*(k + n))^2));
%  gamma = gamma + Gn'*g;	
	v = v + Fn'*g;
 end
 
 

