% Test for a gaussian random number generator

% mu,sig,storedval
mu = 0;
sig = 2;%sqrt(0.2);
N = 500;

res = zeros(2*N,1);

for i=0:(N-1)
    v1=0;v2=0;rsq=0;
    while (rsq >= 1.0 || rsq == 0)
        v1 = 2*rand() - 1;
        v2 = 2*rand() - 1;
        rsq = v1*v1 + v2*v2;
    end

    fac = sqrt(-2*log(rsq)/rsq);
    storedval = v1*fac;
    res(2*i+1) = mu+sig*v2*fac;
    res(2*i+2) = mu+sig*v1*fac;
end

figure();
subplot(1,2,1);
hist(res,100);
subplot(1,2,2);
[h,x] = hist(res,50);
h = h/(2*N*(x(2)-x(1)));
int = trapz(x,h);
hN = 1/(sig*sqrt(2*pi))*exp(-0.5*(x-mu).*(x-mu)/(sig*sig));
plot(x,h,'b',x,1/(sig*sqrt(2*pi))*exp(-0.5*(x-mu).*(x-mu)/(sig*sig)),'r--');
title(['Integral over R : ',num2str(int)]);
