function [Lap, Lap_inv] = modelwt4(pm1,pm2,pm3,pm4,nve,nve4, nhe, delx, dely, surf)
%%
%%	 [Lap, Lap_inv] = modelwt(nve, nhe, delx, dely, surf);
%%
%%  Model weighting matrix to compute Laplacian and its inverse
%%  for smoothing in slip inversions.
%%
%%  Input:
%%		nve 	= number of vertical elements
%%		nhe 	= number of horizontal elements
%%		delx 	= length of elements in along strike dimension
%%		dely 	= length of elements in dip dimension
%%		surf 	= 1 if fault breaks free surface ~= 1 otherwise
%%  Output:
%%		Lap 	= finite difference Laplacian in two dimensions
%%		Lap_inv = inverse of Lap

	nhe1=size(pm1,1)/nve;
   nhe2=size(pm2,1)/nve;
   nhe3=size(pm3,1)/nve;
   nhe4=size(pm4,1)/nve4;

	ngrid = nhe * nve + nhe4*nve4;

	Lap = zeros(ngrid);  xpartd = zeros(ngrid); ypartd = zeros(ngrid);
	
%       x-derivative for central part of grid
%       (exclude left & right edges)
  for i  = nve+1:(nhe*nve-nve)
		xpartd(i,i-nve) = + 1.0;
		xpartd(i,i)     = - 2.0;
		xpartd(i,i+nve) = + 1.0;
 	end
    
    for i  = 1+nve:size(pm4,1)-nve
  xpartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+i,size(pm1,1)+size(pm2,1)+size(pm3,1)+i-nve) = + 1.0;
	xpartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+i,size(pm1,1)+size(pm2,1)+size(pm3,1)+i)     = - 2.0;
	xpartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+i,size(pm1,1)+size(pm2,1)+size(pm3,1)+i+nve) = + 1.0;
   end

%       y-derivative for central part of grid
%       (exclude top & bottom edges)
	temp = zeros(nve);
	for i = 2:nve-1
		temp(i,i-1:i+1) = [1 -2 1];
	end
	for j = 1:nhe
		k = (j-1)*nve; 
		ypartd( k+1:k+nve, k+1:k+nve ) = temp;
	end
   
    for k = 0:nhe4-1
      for i  = 2:nve-1
         cnt=nve*k+i;
   	     	ypartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt,size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt-1) = + 1.0;
		ypartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt,size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt)     = - 2.0;
		ypartd(size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt,size(pm1,1)+size(pm2,1)+size(pm3,1)+cnt+1) = + 1.0;
   	end
   end

   
%% 	top edge y-derivative
%%	 for slip breaking the surface
        for i  = 1:nhe
		k = (i-1)*nve+1;
		if surf == 1 
			ypartd( k, k:k+1) = [-1 1];
		else
			ypartd( k, k:k+1) = [-2 1];
		end
 	end


% top edge y-derivative for plane 4
tmp=size(pm1,1)+size(pm2,1)+size(pm3,1);
for i = 1:nhe4
	
	k = i*nve;
	j=(i-1)*nve+1;
	
	ypartd(tmp+j,tmp+j+1) = + 1.0;
	ypartd(tmp+j,tmp+j)     = - 2.0;
	ypartd(tmp+j,k) = + 1.0;
end


%%	left edge x-derivative on plane 1
        for i  = 1:nve
		xpartd(i,i) 	    = - 2.0;
		xpartd(i,i+nve)     = + 1.0;
 	end
	
%%	right edge x-derivative on plane 3
        for i  = (nhe-1)*nve+1:nhe*nve
		xpartd(i,i) 	    = - 2.0;
		xpartd(i,i-nve)     = + 1.0;
 	end

   
% right edge x-derivative for plane 4
   for i  = ((nhe+nhe4)-1)*nve+1:(nhe+nhe4)*nve
		xpartd(i,i) 	    = - 2.0;
		xpartd(i,i-nve)     = + 1.0;
 	end
    
    % left edge x-derivative for plane 4
    tmp=size(pm1,1)+size(pm2,1)+size(pm3,1);
    for i  = 1:nve
		xpartd(tmp+i,tmp+i) 	    = - 2.0;
		xpartd(tmp+i,tmp+i+nve)     = + 1.0;
 	end
    
%%	bottom edge y-derivative
%for planes 2, 3, and 4
for i  = nhe1+1:nhe+nhe4
    k = i*nve; 
	  ypartd( k, k-1:k) = [1 -2];
end
% for plane 1

for i = 1:nhe1
	k = i*nve;
	
	ypartd(k,size(pm1,1)+size(pm2,1)+size(pm3,1)+(i-1)*nve+1) = + 1.0;
	ypartd(k,k)     = - 2.0;
	ypartd(k,k-1) = + 1.0;
end

	Lap = xpartd/delx^2 + ypartd/dely^2;
	Lap_inv = inv(Lap);
