	function [U1,U2, x] = visc_dip_cycle(H,xtip,ztip,L,dip,MT,T,t,tR,N,mu);

% INPUT:
%	H	= Layer thickness
%	xs	= x position of point source
%	zs	= Depth of point source, < H.
%	M	= moment tensor.
%	t   = time since earthquake
%	tR  = relaxation time 
%	N   = Number of points in FFT
%
% OUTPUT:
%	u1	= velocity
%	u2	= velocity
%	x	= x-coordinate corresponding to u
%
% Propagator matrix method for viscoelastic 2-D dip-slip faults
% Kaj Johnson and Paul Segall, sept. 2001

xmax = 20*H;       %range of x
delta = 2*xmax/N;  %sampling interval in x
kmax = 0.5/delta;  %maximum wavenumber
x = linspace(-xmax, xmax,N);
ks = linspace(-kmax,kmax,N)*(2*pi);  
% factor of 2*pi accounts for different definition of transform

% some material constants
nu = 0.25;
lam=-2*nu.*mu./(2*nu-1);
K=lam+2*mu/3;   %bulk modulus
g=lam+2*mu;
B = 1/tR;

u1 = zeros(size(ks));
u2 = zeros(size(ks));

rg=3*10^-3;
G=[1 0 0 0;0 1 0 0;0 0 1 0;0 -rg 0 1];

for j=1:N

	k = ks(j);
    
	%set up A matrix in equation dv=A*v+f
         A=[0  k 1/mu 0; ...
        -k*lam/g 0 0 1/g; ...
        4*k^2*mu*(lam+mu)/g 0 0 k*lam/g; ...
        0 0 -k 0];
  
  
  	%propagator from top of halfspace to surface
	z=0; z0=H;
	Ph = G*expm(A*(z-z0));
   
   %propagator from source to surface
   Pzs=G*Pintegrated(xtip,ztip,L,dip,lam,mu,ks(j));
   
   
	%body forces
	F = [i*MT(1,2)/mu; MT(2,2)/g; i*i*k*MT(1,1) + k*lam*MT(2,2)/g; 0];

	%compute roots
	[s1, s2, s3, a] = roots_right(k,mu,K,lam,B,Ph);
	%[s1, s2, s3, a] = roots_kaj(k,mu,K,B,Ph);
   
   
	%% RETURN VISCOELASTIC PART ONLY
	%% FIRST root:
	V = makeV_right(k,K,lam,s1,B,H,mu);
   M = Ph(3:4,:)*V;
	Adj(1,1) = M(2,2); Adj(2,2) = M(1,1); 
		Adj(1,2) = -M(1,2); Adj(2,1) = -M(2,1);
	Num1 = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den1 = a*s1*(s1-s2)*(s1-s3);
   Tterm=exp(-s1*T)/(exp(-s1*T)-1);
  % term1 = Num1*(exp(s1*t)-1)*Tterm/Den1;
   term1 = Num1*s1*exp(s1*t)*Tterm/Den1;
  
	%% SECOND root:
	V = makeV_right(k,K,lam,s2,B,H,mu);
	M = Ph(3:4,:)*V;
	Adj(1,1) = M(2,2); Adj(2,2) = M(1,1); 
		Adj(1,2) = -M(1,2); Adj(2,1) = -M(2,1);
	Num2 = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den2 = a*s2*(s2-s1)*(s2-s3);
	Tterm=exp(-s2*T)/(exp(-s2*T)-1);
  % term2 = Num2*(exp(s2*t)-1)*Tterm/Den2;
    term2 = Num2*s2*exp(s2*t)*Tterm/Den2;
   
   
	%% THIRD root:
   V = makeV_right(k,K,lam,s3,B,H,mu);
	M = Ph(3:4,:)*V;
	Adj(1,1) = M(2,2); Adj(2,2) = M(1,1); 
		Adj(1,2) = -M(1,2); Adj(2,1) = -M(2,1);
	Num3 = -Ph(1:2,:)*V*Adj*Pzs(3:4,:)*F;
	Den3 =  a*s3*(s3-s1)*(s3-s2);
	Tterm=exp(-s3*T)/(exp(-s3*T)-1);
 %  term3 = Num3*(exp(s3*t)-1)*Tterm/Den3;
  term3 = Num3*s3*exp(s3*t)*Tterm/Den3;


Fhat = term1 + term2 + term3;

u1(j) = -i*Fhat(1);
u2(j) = -Fhat(2);
   
end %j (loop over wave number)


%filter out high frequency noise at large k
temp=real(u1(N/2+1:N));
temp1=mean(abs(temp(end-10:end)));

if temp1>1e-3

j=0;
q=1;
while q>10^-5
   j=j+1;
   q=mean(abs(temp(5*(j-1)+1:5*j)));
end
u1(N/2+1+j*5:end)=0;
u1(1:N/2-j*5)=0;
u2(N/2+1+j*5:end)=0;
u2(1:N/2-j*5)=0;
end

      %inverse FT into physical space

    %divide by sampling interval to turn analytical transform
    %shift negative frequencies to N/2+1 to N-1

    %inverse FFT
	u1s=ifft(fftshift(u1/delta));
	u2s=ifft(fftshift(u2/delta));

    %shift the solution back
	U1=real(fftshift(u1s));%
	U2=real(fftshift(u2s));


  


