function W = dist_double_power_wts(L,k1,k2,b);

%DIST_POWER_WTS computes power weights of the form
% 
%    Wij = [1 - (dij/b)^k1]^k2 , dij <= b
%        = 0  ,  dij > b
% 
% for any integer k1,k2 > 1, where b > Dmin = bandwidth. W is then 
% normalized by its maximum eigenvalue.

%Written by: TONY E. SMITH, 9/29/07

%INPUTS:     L = (Nx2) coordinate matrix (Xi,Yi),i=1:N
%        k1,k2 = integers  > 1
%            b = bandwidth  
%                (NOTE: will get error message if b < Dmin)
%
%
%OUTPUTS: W = (NxN) weight matrix (in sparse form).                 

N = length(L(:,1)) ;

D = distance_mat(L);  %Make distance matrix

d = max(max(D));

d_pos = min(min(D + d*eye(N)));  %calculates min positive distance  

if b < d_pos        
    
    error(['bandwidth must be larger than smallest positive distance = ',num2str(d_pos)]);
    
end

U = ones(N); %Unit column vector

W = (U - (D/b).^k1).^k2;

Z = (D <= b);

W = W.*Z; %makes all wij's with dij > b equal to zero.

max_eig = max(eig(W));

W = W/max_eig;