function pot1 = op_by_cgpot(pot1, pot2, op)
% OP_BY_CGPOT Compute pot1 op= pot2 for Conditional Gaussian potentials.
% pot1 = op_by_cgpot(pot1, pot2, op)
% 
% op is '*' or '/'


pot1 = cg_mom_to_can(pot1);
pot2 = cg_mom_to_can(pot2);

mask = find_equiv_posns(pot2.ddom, pot1.ddom);
if op == '*'
  for i=1:pot1.dsize
    if isempty(pot2.ddom)
      j = 1;
    else
      sub = ind2subv(pot1.dsizes, i);
      j = subv2ind(pot2.dsizes, sub(mask));
    end
    pot1i = pot1.pot{i};
    pot2j = pot2.pot{j};
    [pot1.pot{i}.g, pot1.pot{i}.h, pot1.pot{i}.K] =  multiply_by_canonical(...
        pot1i.g, pot1i.h, pot1i.K, pot2j.g, pot2j.h, pot2j.K, ...
        pot1.cdom, pot2.cdom, pot1.csizes, pot2.csizes);
  end
elseif op=='/'
  for i=1:pot1.dsize
    if isempty(pot2.ddom)
      j = 1;
    else
      sub = ind2subv(pot1.dsizes, i);
      j = subv2ind(pot2.dsizes, sub(mask));
    end
    pot1i = pot1.pot{i};
    pot2j = pot2.pot{j};
    [pot1.pot{i}.g, pot1.pot{i}.h, pot1.pot{i}.K] = divide_by_canonical(...
        pot1i.g, pot1i.h, pot1i.K, pot2j.g, pot2j.h, pot2j.K, ...
        pot1.cdom, pot2.cdom, pot1.csizes, pot2.csizes);
  end
else
  error(['unrecognized op ' op]);
end
