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

%load data
% xGPS = positions of GPS observations (km) -- distance from arbitrary origin
% dxGPS = vector (column vector) of velocities in x direction (horizontal)
% dzGPS = vector (column vector) of velocities in z direction (vertical)
% sigxGPS = vector (colmun vector) of uncertainties (1 sigma)
% sigzGPS = vector (colmun vector) of uncertainties (1 sigma)
% xLEV = positions of leveling observations (km) -- distance from arbitrary origin
% dzLEV = vector (column vector) of velocities in z direction (vertical)
% sigzLEV = vector (colmun vector) of uncertainties (1 sigma)
% xrefLEV = reference point for leveling data

%replace file synthetic.mat with a mat file containing the above data 
load synthetic.mat %this file is generated with make_synthetic.m

%specify parameters
% convergence = convergence rate
% gamma = smoothing parameter (increasing gamma decreases roughness) 
convergence=80;
gamma=.3; 

%specify fault geometry
%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=150; %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];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

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