function [dist,x,y] = bivar_norm(m,c,bin,plt);

[ndim,nd] = size(m);
maxx = zeros(1,nd);
maxy = zeros(1,nd);
minx = zeros(1,nd);
miny = zeros(1,nd);
densfact = 3.5; % determines how far out in the density to measure (in SD units)
for i = 1:nd
    maxx(i) = m(1,i) + densfact*sqrt(c(1,1,i));
    maxy(i) = m(2,i) + densfact*sqrt(c(2,2,i));
    minx(i) = m(1,i) - densfact*sqrt(c(1,1,i));
    miny(i) = m(2,i) - densfact*sqrt(c(2,2,i));
end

xlims = [min(minx) max(maxx)];
ylims = [min(miny) max(maxy)];

x = xlims(1):bin:xlims(2);
y = ylims(1):bin:ylims(2);
lenx = length(x);
leny = length(y);

for id = 1:nd
    rdcov = sqrt(det(c(:,:,id))); % for normalizing
    incov = inv(c(:,:,id));
    moo = m(:,id);
    for iy = leny:-1:1
        for ix = 1:lenx
            xy = [x(ix); y(iy)];
            zz = xy-moo;
            stat_dist = zz'*incov*zz;
%             if stat_dist > densfact^2
%                 dist(iy,ix,id) = 0;
%             else
                dist(iy,ix,id) = (1/(2*pi*rdcov)) * exp(-0.5*stat_dist);
%             end
        end
    end
end