M = randn(3);
M=M+M'
[X, Y, Z] = quadraytrace(M, [-10,10],[-10,10],[-10,10], 10000);
scatter3(X, Y, Z, rand(numel(X),3)*4)
daspect([1 1 1])
