This function generates a pair of correlated gaussian variates, with
mean zero, correlation coefficient rho and standard deviations
sigma_x and sigma_y in the x and y directions.
The probability distribution for bivariate gaussian random variates is,
p(x,y) dx dy = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp (-(x^2 + y^2 - 2 \rho x y)/2\sigma_x^2\sigma_y^2 (1-\rho^2)) dx dy
for x,y in the range -\infty to +\infty. The
correlation coefficient rho should lie between 1 and
-1.