function [x,y,z,d] = dRandPoints( n, amplitude, nx, nd )

% dRandPoints - gives (Gaussian) random points and distances
%	[x,y,z,d] = dRandPoints( n, amplitude, nx, nd )
%	n  - size of the system
%	amplitude - deviation of the produced coordinates
%	nx - deviation from initial coordinates
%	nd - deviation from accurate distances
%	x,y,z - points coordinates with additional noise nx
%	d - matrix of distances between points with noise nd
%	random initial points are normaly distributed with std=100

if nargin < 4, nd = 0; end;
if nargin < 3, nx = 0; end;
if nargin < 2, amplitude = 100; end;

x = amplitude * randn(n,1);
y = amplitude * randn(n,1);
z = amplitude * randn(n,1);

deltax = x*ones(1,n) - ones(n,1)*x';
deltay = y*ones(1,n) - ones(n,1)*y';
deltaz = z*ones(1,n) - ones(n,1)*z';
d = sqrt(deltax.^2 + deltay.^2 + deltaz.^2);

x = x+nx*randn(size(x));
y = y+nx*randn(size(x));
z = z+nx*randn(size(x));
d = d+nd*randn(size(d));
d = d - diag(diag(d));	%setting diagonal elements to 0

return;