function [U1,U2,x]=ElasVisco2D(H,ztip,xtip,dip,L,slip,t,nu,elastic,gravity)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[U1,U2,x]=ElasVisco2D(H,ztip,xtip,dip,L,slip,t)
%
%Calculates displacements at time t due to uniform slip on a dislocation in
%an elastic layer overlying a Maxwell viscoelastic half-space
%
%input parameters:
%H         thickness of elastic layer (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)
%slip      amount of slip (m)
%t    	  time normalized by tR (tR = relaxation time)
%            NOTE: relaxation time is defined as tR = 2*eta/mu where eta is viscosity
%                  and mu is shear modulus
%elastic = 0 means visco part of solution only (only postseismic deformation)
%elastic = 1 means visco+elastic (coseismic+postseismic)
%
%gravity = 0 means no gravity
%gravity = 1 means inculde gravity
%
%outputs:
%U1        horizontal displacements (m)
%U2        vertical displacements (m)
%x         horizontal position (km)
%
%Note: 1. Fault must be in elastic layer. Solution is not valid if any part of the fault
%extends into the viscoelastic half-space. 2. The solution is only valid within a distance
%of about 4*H from the fault tip.
%  
%Calls functions EdgeDisp.m and visc_dip.m
%
%Kaj Johnson -- last modified Dec. 11, 2002

Ho=H;
tR=1;   %times normalized by tR (relaxation time)

%normalize lengths by H
xtip=xtip/H;
ztip=ztip/H;
L=L/H;
H=1;
mu=1;
dip=dip*pi/180;
D=ztip+L*sin(dip);
N = 1000;		%number of points in FFT

%% 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 continuation of fault) intersects origin
x_intercept=xtip-ztip/tan(dip);
xtip=xtip-x_intercept;
[U1,U2, x] = visc_dip(H,xtip,ztip,L,dip,M,t,tR,N,mu,nu,gravity);

%elastic part of solution
if elastic == 1
m1=[xtip -ztip dip*180/pi slip];
m2=[xtip+L*cos(dip) -(ztip+L*sin(dip)) dip*180/pi -slip];
x2=zeros(size(x));
[u11, u21] =EdgeDisp(m1,x,x2,.25);
[u12, u22] =EdgeDisp(m2,x,x2,.25);
u1=u11+u12;u2=u21+u22;

U1=U1+u1;
U2=U2+u2;
end


x=x+x_intercept; %shift x-coordinate back
x=x*Ho;   %re-dimensionalize H

