function [u1,u2,x] = drive_visc_cycle(H,ztip,xtip,dip,L,SlipRate,mu,t,T)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[u1,u2,x] = drive_visc_cycle(H,ztip,xtip,dip,L,SlipRate,mu,t,T)
%
% Program to calulate surface velocities due to periodic slip on a 2D fault in
% an elastic layer overlying a viscoelastic half-space
%
% NOTES: 1. relaxation time is defined as tR = 2*eta/mu where eta is viscosity
%           and mu is shear modulus 
%        2. Fault must be in elastic layer. Solution is not valid if any part of the fault
%           extends into the viscoelastic half-space. 
%        3. The solution is only valid within a distance
%           of about 4*H from the fault tip.
%%%input parameters:
%H             depth (POSITIVE) to bottom of layer1 (km)
%ztip          depth (POSITIVE) of fault tip (km)
%xtip          horizontal position of fault tip (km)
%dip           dip of fault (degrees)
%L             length of fault (km)
%SlipRate      slip rate -- SlipRate*T = displacement
%mu				shear modulus
%t  	       time normalized by tR (tR = relaxation time)
%T             recurrance interval normalized by tR
%
%%%Outputs:
%u1,u2         horizontal and vertical VELOCITIES (per relaxation time) at 500 points
%              spaced evenly from -40H to 40H
%              NOTE: Solution is only correct within about +/-5*H km from the fault
%
% Kaj Johnson, June 2002
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


slip = SlipRate*T;

tR = 1;		% relaxation time
dip=dip*pi/180;
D=ztip+L*sin(dip);

N = 500;




%% set up moment tensor
M = zeros(2,2);
mus = mu;
M(1,1)=slip*cos(dip)*sin(dip)*2*mus;
M(1,2)=slip*sin(dip)*sin(dip)*mus-slip*cos(dip)*cos(dip)*mus;
M(2,2)=-slip*sin(dip)*cos(dip)*2*mus;

%shift coordinate system such that fault (or projection of fault) intersects origin
x_intercept=xtip-ztip/tan(dip);
xtip=xtip-x_intercept;

[u1,u2, x] = visc_dip_cycle(H,xtip,ztip,L,dip,M,T,t,tR,N,mu);

x=x+x_intercept; %shift x-coordinate back
