%This program generates a synthetic data set and saves the synthetic data
%in synthetic.mat. This mat file is loaded into setup_inversion.m
%This program also plots the fault geometry
%You can change the input parameters.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT SECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%locations of surface observations
xGPS=linspace(10,200,200);
sigGPS=2;  %uncertainties in mm/yr
xLEV=linspace(50,100,50);
sigLEV=2;  %uncertainties in mm/yr
xrefLEV=50; %reference coordinate for leveling data

%fault geometry
dipi=15; %dip of interface
%positive dip indicates dipping to the east
%negative dip indicates dipping to the west
dipi=15; %dip of interface
%list dips for each fault in a row vector -- must list in order from east
%to west
dip=[40 -40 45 -60];

%positions of faults
xi=0;    %position of trace of interface
xbi=180; %position of bottom of interface
%list positions for each fault in a row vector -- must list in order from east
%to west
xflts=[20 50 70 120];

%number of fault patches
Ni=50;  %number patches on interface
%list number of pathes for each fault in a row vector-- must list in order from east
%to west
N=[10 10 10 10];

%specify slip rates
convergence=80;  %convergence rate in mm/yr
%list slip for each fault in a row vector-- must list in order from east
s=[15 10 20 10];

%number of locked patches down from the surface for each fault
%list number of locked patches for each fault in a row vector-- must list in order from east
nlp=[3 1 4 5];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END OF INPUT SECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Li=(xbi-xi)/cos(dipi*pi/180); %length of interface
pLi=Li/Ni;  %patch length for interface

%length of faults using law of sines
for k=1:length(N)
    if dip(k)>0
        L(k)=(xflts(k)-xi)*sin(dipi*pi/180)/sin((dip(k)-dipi)*pi/180);  %length of fault 1 using law of sines
    else
        L(k)=(xflts(k)-xi)*sin(dipi*pi/180)/sin((180-abs(dip(k))-dipi)*pi/180);  
    end
end
%patch length for faults
pL=L./N;

%set up patches
delxi=pLi*cos(dipi*pi/180);
delzi=pLi*sin(dipi*pi/180);
xtipi=xi:delxi:xi+(Ni-1)*delxi;
ztipi=0:delzi:(Ni-1)*delzi;
for k=1:length(N)
    delx(k)=sign(dip(k))*pL(k)*cos(dip(k)*pi/180);
    delz(k)=abs(pL(k)*sin(dip(k)*pi/180));
    xtip{k}=xflts(k):delx(k):xflts(k)+(N(k)-1)*delx(k);
    ztip{k}=0:delz(k):(N(k)-1)*delz(k);
end

%cut interface into segments
index(1,:)=xtipi<(xtip{1}(end)+delx(1)); %patches on segment 1 of interface
for k=2:length(N)
    index(k,:)=xtipi<(xtip{k}(end)+delx(k)) & xtipi>(xtip{k-1}(end)+delx(k-1)); %patches on segment k of interface
end

%displacements due to slip on interface
for k=1:Ni
    m=[50000 pLi ztipi(k) dipi+180 0 xtipi(k) 0 0 1 0]';
    [Ui,D,S,flag]=disloc3dR12(m,[xGPS;0*xGPS;0*xGPS],1,.25);
    Gix(:,k)=Ui(1,:)';
    Giz(:,k)=Ui(3,:)';
    [UiLEV,D,S,flag]=disloc3dR12(m,[xLEV;0*xLEV;0*xLEV],1,.25);
    [UiLEVref,D,S,flag]=disloc3dR12(m,[xrefLEV;0;0],1,.25);
    GizLEV(:,k)=(UiLEV(3,:)-UiLEVref(3,:))';
end  

%displacements due to slip on faults
for j=1:length(N)
    for k=1:N(j)
        m=[50000 pL(j) ztip{j}(k) dip(j)+90+90*sign(dip(j)) 0 xtip{j}(k) 0 0 1 0]';
        [U,D,S,flag]=disloc3dR12(m,[xGPS;0*xGPS;0*xGPS],1,.25);
        Gx{j}(:,k)=sign(dip(j))*U(1,:)';
        Gz{j}(:,k)=sign(dip(j))*U(3,:)';
        [ULEV,D,S,flag]=disloc3dR12(m,[xLEV;0*xLEV;0*xLEV],1,.25);
        [ULEVref,D,S,flag]=disloc3dR12(m,[xrefLEV;0;0],1,.25);
        GzLEV{j}(:,k)=sign(dip(j))*(ULEV(3,:)-ULEVref(3,:))';
end
end

%matrix of velocities
G=[-Gix;-Giz;-GizLEV];

backslipix=zeros(size(Gix,1),length(N));
backslipiz=zeros(size(Gix,1),length(N));
for k=1:length(N)
    G=[G [-Gx{k};-Gz{k};-GzLEV{k}]];
    backslipix(:,k)=-sum(Gix(:,index(k,:)'),2)*sum(s(k:length(N)));
    backslipiz(:,k)=-sum(Giz(:,index(k,:)'),2)*sum(s(k:length(N)));
    backslipizLEV(:,k)=-sum(GizLEV(:,index(k,:)'),2)*sum(s(k:length(N)));
    forwardslipx(:,k)=sum(Gx{k},2)*s(k);
    forwardslipz(:,k)=sum(Gz{k},2)*s(k);
    forwardslipzLEV(:,k)=sum(GzLEV{k},2)*s(k);
end

ss_vx=-convergence+sum(backslipix+forwardslipx,2); %steady state velocity in x direction
ss_vz=sum(backslipiz+forwardslipz,2); %steady state velocity in z direction
ss_vzLEV=sum(backslipizLEV+forwardslipzLEV,2); 

%backslip
bslip=zeros(size(G,2),1);
bslip(1:20)=convergence-sum(s);
bslip(Ni+1:Ni+nlp(1))=s(1);
for k=2:length(N)
    bslip(Ni+sum(N(1:k-1))+1:Ni+sum(N(1:k-1))+nlp(k))=s(k);
end
bs=G*bslip;
bs_vx=bs(1:length(xGPS));
bs_vz=bs(length(xGPS)+1:2*length(xGPS));
bs_vzLEV=bs(2*length(xGPS)+1:end);

dx=ss_vx+bs_vx;
dz=ss_vz+bs_vz;
dzLEV=ss_vzLEV+bs_vzLEV;

%add normally distributed error
dxGPS=dx+sigGPS*(randn(size(dx))-.5);
dzGPS=dz+sigGPS*(randn(size(dx))-.5);
dzLEV=dzLEV+sigLEV*(randn(size(dzLEV))-.5);

sigxGPS=sigGPS*ones(size(dxGPS));
sigzGPS=sigGPS*ones(size(dzGPS));
sigzLEV=sigLEV*ones(size(dzLEV));

%plot fault geometry
figure
hold on
plot(xtipi,-ztipi,'r')
for k=1:length(N)
plot(xtip{k},-ztip{k})
end
axis equal

%plot data
figure;
subplot(211)
errorbar(xGPS,dxGPS,2*sigxGPS,'.')
xlabel('distance from origin (km)')
ylabel('horizontal velocity')
subplot(212)
errorbar(xGPS,dzGPS,2*sigzGPS,'.')
xlabel('distance from origin (km)')
ylabel('vertical velocity')

figure;
subplot(211)
errorbar(xLEV,dzLEV,2*sigzLEV,'.')
xlabel('distance from origin (km)')
ylabel('vertical velocity')
title('leveling data')

save synthetic xGPS xLEV dxGPS dzGPS dzLEV sigxGPS sigzGPS sigzLEV 
