function Lap = Laplacian(Ni,N)
%%
%%  Model weighting matrix to compute Laplacian for smoothing in slip inversions.
%%
%%  Input:
%%		Ni      = number of elements on interface
%%		N       = vector containing number of elements on each thrust
%%  Output:
%%		Lap 	= finite difference Laplacian in one dimension

numthrusts=length(N);
nel = Ni+sum(N);
Lap = zeros(nel);
    
%       derivative for interface
%       (exclude first & last patch)
for i  = 2:Ni-1
    Lap(i,i-1) = + 1.0;
    Lap(i,i)     = - 2.0;
    Lap(i,i+1) = + 1.0;
end
% 	top edge derivative for fault breaking surface
    Lap( 1, 1:2) = [-1 1];
%	bottom edge derivative - smooth to zero
    Lap( Ni, Ni-1:Ni) = [1 -2];

NN=cumsum([0 N]);
for k=1:numthrusts
    for i  = Ni+NN(k)+2:Ni+NN(k+1)-1
        Lap(i,i-1) = + 1.0;
        Lap(i,i)     = - 2.0;
        Lap(i,i+1) = + 1.0;
    end
    % 	top edge derivative for fault breaking surface
        Lap( Ni+NN(k)+1, Ni+NN(k)+1:Ni+NN(k)+2) = [-1 1];
    %	bottom edge derivative - smooth to zero
        %Lap( Ni+NN(k+1), Ni+NN(k+1)-1:Ni+NN(k+1)) = [1 -2];
%        Lap( Ni+NN(k+1), Ni+NN(k+1)-1:Ni+NN(k+1)) = [-1 1];
end
    

    
            