function f = bnormpdf(x,th,th2,p);
%
%
%
%
%

m11 = th(1);
m12 = th(2);
sig11 = th(3);
sig12 = th(4);
rho1 = th(5);

[n n2] = size(x);
x(:,1) = (x(:,1)-m11*ones(n,1))/sig11;
x(:,2) = (x(:,2)-m12*ones(n,1))/sig12;
quad = (x(:,1).*x(:,1) + x(:,2).*x(:,2) - 2*rho1*x(:,1)*x(:,2))/(1-rho1^2);
f = exp(-quad/2)/(2*pi*sig11*sig12);

if nargin >= 3

  m21 = th2(1);
  m22 = th2(2);
  sig21 = th2(3);
  sig22 = th2(4);
  rho2 = th2(5);

  [n n2] = size(x);
  x(:,1) = (x(:,1)-m21*ones(n,1))/sig21;
  x(:,2) = (x(:,2)-m22*ones(n,1))/sig22;
  quad = (x(:,1).*x(:,1) + x(:,2).*x(:,2) - 2*rho2*x(:,1)*x(:,2))/(1-rho2^2);
  f2 = exp(-quad/2)/(2*pi*sig21*sig22);

  f = (1-p)*f + p*f2;

end
