function y =  rbfDLaplace (rr)   % Double Laplace of rbf, see rbflap2.mws
global RBFscale;
global RBFpar;
global RBFtype;   % 1=MQ, 2=G, 3=IMQ, 4=W2, 5=TPS, 6=W4
r=rr/RBFscale;
if RBFtype==1
    y=(-(-8+8*r.*r+r.^4).*(1+r.*r).^(-7/2))/RBFscale^4;     % Multiquadric
    return;
end
if RBFtype==2
    y=(16*(r.^4-4*r.*r+2).*exp(-r.*r))/RBFscale^4;       % Gaussian
    return;
end
if RBFtype==3
    y=3*((8-24*r.*r+3*r.^4).*(1+r.*r).^(-9/2))/RBFscale^4;    % inverse Multiquadric
    return;
end
if RBFtype==4
  y=(60*(-16*r+15*r.^2+3).*max(0,(1-r))).*(1/r)./(RBFscale^4);   % Wendland 2
   return;
end
if RBFtype==5
    y=0;                % TPS
    return;
end
if RBFtype==6
    y=(6720*(4*r-1).*(3*r-2).*max(0,(1-r)).^2)/RBFscale^4;  % Wendland C4
    return;
end
if RBFtype==7
    y=16.*RBFpar.^2.*(2*r.*r-RBFpar.*RBFpar)./((r+RBFpar).^4.*RBFscale^4);      % regularized log = AFS of Laplace in R^2
    return;
end

error('RBF type not implemented')