function f = multigaussian(x,mu,S)
  % Returns a multivariate Gaussian density
  % mu must be a 1xd row vector, S a dxd square covariance matrix
  % x a nxd matrix, with one point per row
  % returns a nx1 column vector of densities, calculated at these points
  [n,d] = size(x);
  xmu = op(x,'-',mu) ;
  f = (1/2*pi)^(d/2)*exp(-sum((xmu*inv(S)).*xmu,2)/2)/sqrt(det(S));