function v=antiplane_viscflt_main(x,H,D,eta,t,tR,T,s_rate,npsc,nti)

% v=antiplaneviscflt(x,H,D,eta,t,tR,T,npsc,nti)
%
% Calculates velocities at points in x for a fault with cyclic earthquakes down
% to depth D and creeping as a linear viscous fault zone from depth D to H.
% Fault is infintely long strike-slip fault (antiplane) in an elastic layer
% over Maxwell viscoelastic half-space. Implementation of solution outlined
% in Johnson and Segall, 2004, Viscoelastic earthquake cycle models with deep
% stress-driven creep along the San Andreas fault system, JGR, 109,
% 10.1029/2004JB003096
%
% inputs:
%		x = (vector) locations of observations (km)
%		H = (scalar) thickness of elastic layer (km)	
%	 	D = (scalar) depth to locked fault/creeping fault intersection (km)
%  s_rate = (scalar) long term slip rate (coseismic slip/recurrence time) (m/yr)
%	  eta = (scalar) viscosity/width of fault zone (normalized by half-space
%	        viscosity) (1/m)
%	    t = (scalar) time since last earthquake (years)
%      tR = (scalar) relaxationtime of half-space (defined as 2*eta/mu, eta viscosity,
%           mu shear modulus) (years)
%       T = (scalar) recurrence interval (years)
%    npsc = (scalar) number of patches on creeping fault (experiment to determine
%    how many needed, 10 generally works well)
%     nti = (scalar) number of time increments (experiment to determine how many
%           needed, 100 generally works well but more may be needed for
%           low fault zone viscosities and observation times early in the 
%           cycle. Fewer might be used for observation times late in the cycle)
% output:
%		v = (vector) velocities at points in x (m/yr)
%
%  Calls functions: makeTs.m, savpres_depth.m, savpres_patch.m
% Kaj Johnson, last modified May 30, 2003


T=T/tR;  %normalize times by relaxation time of half-space
t=t/tR;

ztipc=D;

mu=1;

delt=T/nti; %length of time increments
tspace=linspace(delt,T,nti)-.55*delt;  %incremental times

Lc=H-ztipc;	%length of coseismic fault
dzc=Lc/npsc;  %length of fault patches	
d2=ztipc:dzc:H-dzc;  %fault patch ends
d1=ztipc+dzc:dzc:H;
center=ztipc+dzc/2:dzc:H-dzc/2;  %center of fault patches

%get stresses acting on lower fault due to slip on upper fault
[dummy, SIG] = savpres_depth(0,-center,0,ztipc,tspace,H,1,T);
[v, dummy] = savpres_patch(x,0,ztipc,tspace,H,1,T);
%get stresses acting on lower fault due to slip on lower fault patches
for k=1:npsc
   [dummy, SIGL(:,k,:)] = savpres_depth(0,-center,d2(k),d1(k),tspace,H,1,T);
   [vL(:,k,:), dummy] = savpres_patch(x,d2(k),d1(k),tspace,H,1,T);
end

%convert from stress rate to stress
SIG=SIG*delt;
SIGL=SIGL*delt;

%calculate "driving" stress (i.e. stresses due to repeated coseismic slip)
	%elastic stresses on patches due to coseismic slip
	cosigma=-mu/(2*pi)*(1./(-center+ztipc)-1./(-center-ztipc)-1./(-center+0)+1./(-center-0));
   SIG(:,1)=SIG(:,1)+T*cosigma'; %coseismic elastic stress + viscoelastic stresses at first time interval
	for k=2:size(SIG,2)
   	SIG(:,k)=SIG(:,k)+SIG(:,k-1);  %accumulate stresses over the time intervals
	end
	sigma=SIG(:);  %vector of "driving stresses"


%elastic component of stresses on patches due to slip on patches
for j=1:length(d1)
s13(j,:)=-T*mu/(2*pi)*(1./(-center(j)+d1)-1./(-center(j)-d1)-1./(-center(j)+d2)+1./(-center(j)-d2));
end

%setup G matrix in sliprate=inv(ETA-delt*G)*sigma;
order=[];
for k=1:nti
   order=[order;(k-1):-1:1 nti:-1:k];
end
for loop=1:nti
   temp=[];
   for k=1:nti
      temp=[temp SIGL(:,:,order(loop,k))];
   end
   G(npsc*(loop-1)+1:npsc*loop,:)=temp;  %viscoelastic part only
end
%add in elastic part
S13=zeros(nti*npsc,nti*npsc);
 for loop=1:nti;
    S13(npsc*(loop-1)+1:npsc*loop,npsc*(loop-1)+1:npsc*loop)=s13;
end
G=G+S13;
%accumulate stresses over time intervals
for k=1:nti-1
    G(npsc*k+1:npsc*(k+1),:)=G(npsc*k+1:npsc*(k+1),:)+G(npsc*(k-1)+1:npsc*(k),:);
end

%invert system to find slip rates at all times on all patches
ETA=eye(size(G))*eta;
sliprate=inv(ETA-delt*G)*sigma;
sliprate=reshape(sliprate,npsc,nti);
slip=sliprate*delt;  %convert from slip rate to slip

%convert from slip per cycle to slip per tR
slip=slip*T;

%%%calculate velocity
order=[];
for k=1:nti
   order=[order;(k):-1:1 nti:-1:k+1];
end

S=zeros(length(x),npsc,nti);
vel=zeros(length(x),nti);

%elastic part of velocity (slip is over time interval of length delt, so divide displacement by delt)
for j=1:npsc
  ve(:,j)=1/pi*(atan(x'/d2(j))-atan(x'/d1(j)))/delt;
end
   
SLIP=slip';
for loop=1:nti
   for k=1:nti
      S(:,:,k)=repmat(SLIP(order(loop,k),:),length(x),1);
   end
   temp=sum(vL.*S,3);  %viscoelastic contribution to velocity
   temp=sum(temp,2)/T;  %divide by T to convert to slip rate
   vel(:,loop)=temp + v(:,loop)+ ve*SLIP(loop,:)'; %slip on patches + slip on coseismic fault + elastic part
end

%interpolate to find velocity at observation times
temp=interp1(tspace,vel',t);
v=temp';

v=v*s_rate;  %scale by coseismic slip magnitude

