function T1 = op_by_table(T1, T2, op, algo)
% OP_BY_TABLE Compute T1(i) := T1(i) op T2(i) elementwise.
% T1 = op_by_table(T1, T2, op, algo)
% 
% op can be one of '*', '/' or '+'.
% If T2's domain is smaller than T1's, it will be extended first.

if ~isequal(T1.domain, T2.domain)
  B = extend_domain(T2.T, T2.domain, T2.sizes, T1.domain, T1.sizes);
else
  %B = T2.T; % if T2.sizes = [1 1 2 1], size(T2.T) = [2 1] which is incompatible
  B = extend_domain(T2.T, T2.domain, T2.sizes, T1.domain, T1.sizes);
end
if op == '*'
  T1.T = T1.T .* B;
elseif op == '/'
  %tiny = exp(-700);
  %T1.T = T1.T ./ (B + tiny);
  % Replace 0s by 1s before dividing.
  B = B + (B==0);
  T1.T = T1.T ./ B;
elseif op == '+'
  T1.T = T1.T + B;
else
  error(['unrecognized op ' op]);
end


