% boxmuller.m
% demonstrates Box-Muller algorithm to create 
% normally distributed random numbers

n = 1e6; z1 = rand(1,n); z2 = rand(1,n);

r1 = sqrt(-2*log(z1)).*cos(2*pi*z2); 
r2 = sqrt(-2*log(z1)).*sin(2*pi*z2);

dx = 0.1; x = -10:dx:10; % create x range

ps1 = hist(r1,x);        % create histogram, centers given by x
ps1 = ps1/sum(ps1)*1/dx; % normalize

ps2 = hist(r2,x);        % create histogram, centers given by x
ps2 = ps2/sum(ps2)*1/dx; % normalize

plot(x,ps1,x,ps2);       % plot