function [g1, h1, K1] = divide_by_canonical(g1, h1, K1, g2, h2, K2, domain1, domain2, sizes1, sizes2)
% DIVIDE_BY_CANONICAL Compute pot1 /= pot2 for Gaussian potentials in canonical form.
% [g1, h1, K1] = divide_by_canonical(g1, h1, K1, g2, h2, K2, domain1, domain2, sizes1, sizes2)

g1 = g1 - g2;
if sum(sizes2) ~= 0 % ~isempty(domain2)
  mask = find_equiv_posns(domain2, domain1); % dom2 is smaller
  u = block(mask, sizes1);
  h1(u) = h1(u) - h2;
  K1(u, u) = K1(u, u) - K2;
end
