function P = gaussian(k)
%               P = gaussian(k)
% Calculates the probability that a measurement is within
% k standard deviations from the mean for a Guassian density.
% It assumes zero mean and identity covariance matrix.
% Input: k - integrate out to k standard deviations.
%        P - probability for measurement to fall within k
%            standard deviations from the (zero) mean.
% Integrate only over the first quadrant.
% Use the midpoint rule for integration.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 100; % Specify number of grid points.
d = k/N;
y = 0.5*d:d:k;
x = y;
for i = 1:N-1
    up = round(N*(sqrt(1-(i/N)*(i/N))));
    X = x(1:up);
    ex = exp(-0.5*X.*X);
    s(i) = sum(ex);
end
ey = exp(-0.5*y(1:N-1).*y(1:N-1))*2*k*k/(pi*N*N);
P = sum(ey.*s);