function output = wsqs(args);  %, pts, vario, np);
global pts vario np; 
range = args(1);
sill = args(2);

% n = length(np);

% Gaussian variog.

theor = sill*(1.5*(pts/range) - 0.5*(pts/range).^3).*(pts < range) + ...
          sill*(pts >= range);

output = sum((theor - vario).^2.*np);