xplot=linspace(0.001,max(xGPS)+10);
%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,[xplot;0*xplot;0*xplot],1,.25);
    Gixp(:,k)=Ui(1,:)';
    Gizp(:,k)=Ui(3,:)';
    [UiLEVref,D,S,flag]=disloc3dR12(m,[xrefLEV;0;0],1,.25);
    GizLEVp(:,k)=(Ui(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,[xplot;0*xplot;0*xplot],1,.25);
        Gxp{j}(:,k)=sign(dip(j))*U(1,:)';
        Gzp{j}(:,k)=sign(dip(j))*U(3,:)';
        [ULEVref,D,S,flag]=disloc3dR12(m,[xrefLEV;0;0],1,.25);
        GzLEVp{j}(:,k)=(U(3,:)-ULEVref(3,:))';
end
end


G=[-Gixp;-Gizp;-GizLEVp];

backslipixp=zeros(size(Gixp,1),length(N));
backslipizp=zeros(size(Gixp,1),length(N));
backslipizLEVp=zeros(size(GizLEVp,1),length(N));
for k=1:length(N)
   G=[G [-Gxp{k};-Gzp{k};-GzLEVp{k}]];
    for j=1:k
        backslipixp(:,k)=backslipixp(:,k)-sum(Gixp(:,index(j,:)'),2);
        backslipizp(:,k)=backslipizp(:,k)-sum(Gizp(:,index(j,:)'),2);
        backslipizLEVp(:,k)=backslipizLEVp(:,k)-sum(GizLEVp(:,index(j,:)'),2);
 end
    forwardslipxp(:,k)=sum(Gxp{k},2);
    forwardslipzp(:,k)=sum(Gzp{k},2);
    forwardslipzLEVp(:,k)=sum(GzLEVp{k},2);
end

ss_vx=backslipixp+forwardslipxp; %steady state velocity in x direction
ss_vz=backslipizp+forwardslipzp; %steady state velocity in z direction
ss_vzLEV=backslipizLEVp+forwardslipzLEVp; 

G=[G [ss_vx;ss_vz;ss_vzLEV]];

dhat=G*backslip;
dhatxGPS=dhat(1:length(xplot))-convergence;
dhatzGPS=dhat(length(xplot)+1:2*length(xplot));
dhatzLEV=dhat(2*length(xplot)+1:end);

%plot fit to data
figure;
subplot(211)
errorbar(xGPS,dxGPS,2*sigxGPS,'.')
hold on
plot(xplot,dhatxGPS,'r')
xlabel('distance from origin (km)')
ylabel('horizontal velocity')
title('fit to GPS data')
subplot(212)
errorbar(xGPS,dzGPS,2*sigzGPS,'.')
hold on
plot(xplot,dhatzGPS,'r')
xlabel('distance from origin (km)')
ylabel('vertical velocity')

figure;
subplot(211)
errorbar(xLEV,dzLEV,2*sigzLEV,'.')
hold on
plot(xplot,dhatzLEV,'r')
xlabel('distance from origin (km)')
ylabel('vertical velocity')
title('fit to leveling data')

%plot slip rates

steady=convergence*ones(Ni,1);
for k=1:length(N)
    steady(index(k,:))=convergence-sum(backslip(end-(length(N)-k):end));
end

figure
plot(steady,-ztipi)
hold on
plot(steady-backslip(1:Ni),-ztipi,'r')
plot(backslip(1:Ni),-ztipi,'g')
    xlabel('slip rate')
    ylabel('depth')
title('interface')
legend('long-term','interseismic','back-slip')

figure
subplot(221)
    plot(backslip(end-(length(N)-1))*ones(size(ztip{1})),-ztip{1})
    hold on
    plot(backslip(end-(length(N)-1))*ones(length(ztip{1}),1)-backslip(Ni+1:Ni+N(1)),-ztip{1},'r')
    plot(backslip(Ni+1:Ni+N(1)),-ztip{1},'g')
    xlabel('slip rate')
    ylabel('depth')
    title('fault 1')

for k=2:length(N)
    if mod(k,4)==1
        figure
    end
   if mod(k,4)==0
        subplot(224)
   else
    subplot(2,2,k)
   end
    plot(backslip(end-(length(N)-k))*ones(size(ztip{k})),-ztip{k})
    hold on
    plot(backslip(end-(length(N)-k))*ones(length(ztip{k}),1)-backslip(Ni+sum(N(1:k-1))+1:Ni+sum(N(1:k))),-ztip{k},'r')
    plot(backslip(Ni+sum(N(1:k-1))+1:Ni+sum(N(1:k))),-ztip{k},'g')
    xlabel('slip rate')
    ylabel('depth')
    title(['fault ' num2str(k)])
end

legend('long-term','interseismic','back-slip')