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

%fault 1 is upper ramp
%fault 4 is lower decollement
%faults 2 and 3 make northern bend
dip1=26; %dip of fault 1
dip4=0; %dip of fault 4
D=7.7; %depth to top of fault 4

%define some parameters
W=D/sin(dip1*pi/180);
tmp1=W*cos(dip1*pi/180)*cos(3.3*pi/180);
tmp2=W*cos(dip1*pi/180)*sin(3.3*pi/180);


faults(1,:)=[70 W 0 -180+dip1 3.3 3.5175 -38.15];
faults(2,:)=[8.58 W 0 -180+23 35 9.397 1.8055];
faults(3,:)=[25 W 0 -180+23 88 24.77 7];
faults(4,:)=[70 30 D -180+dip4 3.3 3.52+tmp1 -38.15-tmp2];


%shift the E and N to local Taiwan coordinate system
faults(:,6)=faults(:,6)+.897;
faults(:,7)=faults(:,7)-0.1945;
	 
%change faults into Okada format
dis_geom1  = [faults, [1 1 0; 1 1 0; 1 1 0; 1 1 0]];
dis_geom = movefaultnew(dis_geom1);
   
     
 %% Setup Fault Geometry and Smoothing operators
        nhe1 = 10; nve1 = 5;
        pm1 =patchfault(dis_geom(1,1:7),nhe1,nve1);
        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)))];
        
        nhe2 = 3; nve2 = 5;
        pm2 =patchfault(dis_geom(2,1:7),nhe2,nve2);
        pm2ss = [pm2 ones(size(pm2(:,1))) zeros(size(pm2(:,1))) zeros(size(pm2(:,1)))];
        pm2ds = [pm2 zeros(size(pm2(:,1))) ones(size(pm2(:,1))) zeros(size(pm2(:,1)))];

        
        nhe3 = 6; nve3 = 5;
        pm3 =patchfault(dis_geom(3,1:7),nhe3,nve3);
        pm3ss = [pm3 ones(size(pm3(:,1))) zeros(size(pm3(:,1))) zeros(size(pm3(:,1)))];
        pm3ds = [pm3 zeros(size(pm3(:,1))) ones(size(pm3(:,1))) zeros(size(pm3(:,1)))];

        nhe4 = 10; nve4 = 5;
        pm4 =patchfault(dis_geom(4,1:7),nhe4,nve4);
        pm4ss = [pm4 ones(size(pm4(:,1))) zeros(size(pm4(:,1))) zeros(size(pm4(:,1)))];
        pm4ds = [pm4 zeros(size(pm4(:,1))) ones(size(pm4(:,1))) zeros(size(pm4(:,1)))];

        pm = [pm1; pm2; pm3; pm4];


%% Compute Kernel for distributed slip calculation
        nu = 0.25;
 
       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
       for i=1:size(pm2ss,1)
          tmp1 = disloc(pm2ss(i,:)', xy, nu);
          tmp2 = disloc(pm2ds(i,:)', xy, nu);
          G21(:,i)=tmp1(:);
          G22(:,i)=tmp2(:);
       end
       for i=1:size(pm3ss,1)
          tmp1 = disloc(pm3ss(i,:)', xy, nu);
          tmp2 = disloc(pm3ds(i,:)', xy, nu);
          G31(:,i)=tmp1(:);
          G32(:,i)=tmp2(:);
       end
       for i=1:size(pm4ss,1)
          tmp1 = disloc(pm4ss(i,:)', xy, nu);
          tmp2 = disloc(pm4ds(i,:)', xy, nu);
          G41(:,i)=tmp1(:);
          G42(:,i)=tmp2(:);
       end


        %%%NEED TO REWRITE THE LAPLACIAN; This is an approximation assuming
                %%equal element sizes for all elements
		%also does not measure nearest neigbors but uses indicies only
       [Lap, Lap_inv] = modelwt4(pm1,pm2,pm3, pm4, nve1,nve4, nhe1+nhe2+nhe3, dis_geom(1,1)/nhe3, dis_geom(1,2)/nve3, 1);
     % [Lap, Lap_inv] = modelwt_shallow(pm1,pm2,pm3, pm4, nve1,nve4, nhe1+nhe2+nhe3, dis_geom(3,1)/nhe3, dis_geom(3,2)/nve3, 1);
  %    [Lap, Lap_inv, zeroslip] = attachhzpg(pm1,pm2,pm3, pm4, nve1,nve4, nhe1+nhe2+nhe3, dis_geom(3,1)/nhe3, dis_geom(3,2)/nve3, 1);
  % 	[Lap, Lap_inv, zeroslip] = detachhzpg(pm1,pm2,pm3, pm4, nve1,nve4, nhe1+nhe2+nhe3, dis_geom(3,1)/nhe3,dis_geom(3,2)/nve3, 1);                                             
                                             
	G = [G11,G21, G31, G41, G12, G22, G32, G42];
	Lapp = [Lap, zeros(size(Lap)); zeros(size(Lap)), Lap];





%make covariance matrix
%not this is not the covariance matrix!! dcov contains the uncertainties (sigma) not sigma^2
dcov=diag(dcov(1:Nsites*3));
	
%nnlsq4
nnlsq4_cycle
TIME=cputime-TIME