nhe=20;
nve=20;

%enter fault geometry
%faults=[length, width, *depth, dip, strike(degrees), *north offset, *east offset]
%*depth to top edge, north and east offsets refer to location of center of top edge
faults=[100 100 0 -150 0 0 0];
    
pm=[];

for k=1:size(faults,1)
%specify components of slip to be calculate ([strike-slip,dip-slip,opening]) -- e.g. [0 1 0] means dip slip only
dis_geom1  = [faults(k,:), [1 1 0]];
dis_geom = movefault(dis_geom1);  % move the fault so that the coordinates of the midpoint refer to the
											 % fault bottom as in Okada

	
%% Create slip patches
pm = [pm; patchfault(dis_geom(1,1:7),nhe,nve)];

end


%create site coordinates
x=linspace(-100,100,8);
y=linspace(-100,100,8);
[X,Y]=meshgrid(x,y);
xysites=[X(:) Y(:)];

coslip=zeros(2*nhe*nve,1);
for k=7:13
coslip(nhe*nve+nve*k+1:nhe*nve+nve*k+5)=5;
end

%create rate-state frictional parameters
sigAmB=.5*10^6*ones(size(pm,1),1);

%pre-earthquake sliprate
prerate=zeros(2*nhe*nve,1);
prerate(end/2+1:end)=0.01;  %m/yr

%construct stress and displacement matrices
setup_hotfriction

%construct synthetic data
stopt=1;
nafterpatches=size(pm,1);
numsites=size(xysites,1);

%positions before the earthquake
east_pos=.1*rand(numsites,1);
north_pos=.1*rand(numsites,1);
up_pos=.1*rand(numsites,1);

%interseismic velocities (m/yr)
inter_east=.00025*xysites(:,1);
inter_north=zeros(size(inter_east));
inter_up=zeros(size(inter_east));

%coseismic displacements
codisps=dispG*coslip;

%afterslip
costress=stressG*coslip;
[allslips,alltimes]=get_slip_hotfriction(prerate,costress,nafterpatches,stressG,stopt,sigAmB);
alltimes(end)=stopt;

cumslip=cumsum(allslips,2);
for k=1:size(cumslip,2)
    disps(:,k)=dispG*cumslip(:,k); 
end

%viscous flow
H=30;tR=1.5;
t_visco=0:stopt/10:stopt;
[Ux,Uy,Uz]=get_viscous_multitimes(coslip,H,stopt,tR,xysites',pm);


%construct time series as sum afterslip + pre-eq. position + interseismic
for k=1:size(disps,1)/3
    sigeast=.003;
    signorth=0.0015;
    sigup=0.008;
    disps(3*k-2,:)=disps(3*k-2,:)+sigeast*randn(size(disps(3*k-2,:)))+east_pos(k)+codisps(3*k-2)+inter_east(k)*alltimes+interp1(t_visco,Ux(:,k)',alltimes);
    disps(3*k-1,:)=disps(3*k-1,:)+signorth*randn(size(disps(3*k-1,:)))+north_pos(k)+codisps(3*k-1)+inter_north(k)*alltimes+interp1(t_visco,Uy(:,k)',alltimes);
    disps(3*k,:)=disps(3*k,:)+sigup*randn(size(disps(3*k,:)))+up_pos(k)+codisps(3*k)+inter_up(k)*alltimes+interp1(t_visco,Uz(:,k)',alltimes);
    
    InvertSeries{k,1}=[alltimes(150:end)' disps(3*k-2,150:end)' sigeast*ones(size(alltimes(150:end)'))];
    InvertSeries{k,2}=[alltimes(150:end)' disps(3*k-1,150:end)' signorth*ones(size(alltimes(150:end)'))];
    InvertSeries{k,3}=[alltimes(150:end)' disps(3*k,150:end)' sigup*ones(size(alltimes(150:end)'))];
end

%some sites will have no pre-earthquake positions
for k=1:6:numsites
    east_pos(k)=nan;
    north_pos(k)=nan;
    up_pos(k)=nan;
end


%site names
names=['01';'02';'03';'04';'05';'06';'07';'08';'09';...
    '10';'11';'12';'13';'14';'15';'16';'17';'18';'19';...
    '20';'21';'22';'23';'24';'25';'26';'27';'28';'29';...
    '30';'31';'32';'33';'34';'35';'36';'37';'38';'39';...
    '40';'41';'42';'43';'44';'45';'46';'47';'48';'49';...
    '50';'51';'52';'53';'54';'55';'56';'57';'58';'59';...
    '60';'61';'62';'63';'64'];

    
%time of earthquake
teq=0;


