%load the data
%d is vector of displacements (e,n,u)
%sig is vector of standard deviations
load d.txt; dtemp=[d(:,2) d(:,1) d(:,3)]; d=dtemp'; d=d(:);
load sig.txt; sig=sig'; sig=sig(:);
dsig=diag(sig);

%load site names
fid = fopen('sites.txt','r');
ff = fscanf(fid, '%s', [4 inf]);
sites = ff'; fclose(fid);
Nsites = size(sites,1);

%load site coorinates
load coords.txt;
llh = zeros(size(coords'));  llh(1,:) = coords(:,1)'; llh(2,:) = coords(:,2)';
%transform from lat,long to x,y    
or = coords(61,:)';		%% origin of fault model is AF26
origin = [or(1), or(2)];
xy = llh2localxy(llh, origin);

%setup fault geometry
nhe = 20; %number of fault elements along horizontal dimension
nve = 8; %number of fault elements along vertical dimension

%vector of fault parameters
%faults=[length width depth(top edge) dip strike east(center of top edge) north(center of top edge)
faults=[70.0000 60.0000 0 -180+23 3.3000 4.5+.897 -28-0.1945];
dis_geom1  = [faults, [1 1 0]];
dis_geom = movefaultnew(dis_geom1);

% setup fault patches
pm1 =patchfault(dis_geom(1,1:7),nhe,nve);
pm1ss = [pm1 ones(size(pm1(:,1))) zeros(size(pm1(:,1))) zeros(size(pm1(:,1)))];
pm1ds = [pm1 zeros(size(pm1(:,1))) ones(size(pm1(:,1))) zeros(size(pm1(:,1)))];

% Compute Kernel for distributed slip calculation
nu = 0.25; %poisson's ratio for elastic half-space
for i=1:size(pm1ss,1)
	tmp1 = disloc(pm1ss(i,:)', xy', nu);
	tmp2 = disloc(pm1ds(i,:)', xy', nu);
	G11(:,i)=tmp1(:);
   G12(:,i)=tmp2(:);
end
G = [G11, G12]; %Kernel

%calculate discrete Laplacian operator
[Lap, Lap_inv] = modelwt(nve, nhe, dis_geom(1,1)/nhe,dis_geom(1,2)/nve, 1);
Lapp = [Lap, zeros(size(Lap)); zeros(size(Lap)), Lap];


%do non-negative least squares inversion
gamma=1;  %SPECIFY SMOOTHING PARAMETER
data=[d;zeros(size(Lapp,1),1)];
A=[G;Lapp/(gamma^2)];
%fast non-negative least squares
s=fnnls(A'*A,A'*data); %s is vector of slips (strike-slip,dip-slip)

%calculate least squares displacement
dhat = G*s; %vector of computed displacements
r = (d - dhat); %residual
rnorm = r'*r; %norm of residual

%plot estimated slip distribution
plotpatchslip(pm1, s(1:size(Lap,1))',nve)
title('strike slip')
plotpatchslip(pm1, s(size(Lap,1)+1:end)',nve)
title('dip slip')
%net slip
ns=sqrt(s(1:size(Lap,1)).^2+s(size(Lap,1)+1:end).^2)
plotpatchslip(pm1, ns',nve)
title('net slip')
hold on
%plot slip vectors
slipvects2(pm1, s(1:size(Lap,1)), s(size(Lap,1)+1:end),nve,faults)
title('slip vectors')

%plot fault plane and displacement vectors
plotgeom