/* Contract indices in expression */

/* NOTE: buildq works like splicing a list with quasiquote in lisp; escaped
items are defined at start of quote */
contr(expr, var) ::= (
  if listp(var) and length(var) > 1 then buildq(
    [expr, head : first(var), tail : rest(var)],
    sum(contr(expr, tail), head, 1, 4))
  elseif listp(var) and length(var) = 1 then buildq(
    [expr, head : first(var)], sum(expr, head, 1, 4))
  elseif listp(var) then (
    expr)
  else buildq([expr, var],
    sum(expr, var, 1, 4)));

/* Differentiate with respect to nth coordinate */
cdiff(func, coord) := diff(func, ct_coords[coord])$

/* Wave operator */
box_op(func) := (
  expand(contr(ug[μ,ν] * cdiff(cdiff(func, ν), μ), [μ,ν])
    - contr(ug[μ,ν] * mcs[μ,ν,ρ] * cdiff(func, ρ), [μ,ν,ρ])))$

mappend(f, list) := apply(append, map(f, list));

ensure_list(x) := if listp(x) then x else [x];

/*
WARNING: Preppended "-" means covariant index *not* contravariant like in
         itensor package

TODO: Handle sums
NOTE: Does not work on complex expressions, just indexed arrays
*/
cov_diff(sim,i) ::= (
  if atom(sim) then (
    buildq([sim, i], cdiff(sim, i)))
  elseif op(sim) = "+" then (
    buildq([terms : args(sim), i],
      apply("+", map(lambda([x], cov_diff(x, i)), terms))))
  else block(
    [ten : op(sim), inds : args(sim), covp, terms],
    if not listp(inds) then (
      error(printf(false, "Expected indexed object, got ~a", inds))),
    covp : map(lambda([x], not atom(x) and op(x)="-"), inds),
    inds : map(lambda([x], if atom(x) then x else first(args(x))), inds),
    terms : map(lambda([sgn,index],
        block([k : gensym()],
          /* Covariant */
          if sgn then buildq(
            [ten, inds : subst(index=k, inds), i, index, k],
            - sum(mcs[i, index, k] * ten[splice(inds)],k,1,4))
          /* Contravariant */
          else buildq(
            [ten, inds : subst(index=k, inds), i, index, k],
            sum(mcs[i, k, index] * ten[splice(inds)],k,1,4)))),
      covp, inds),
    buildq(
      [i, ten, inds, terms : apply("+", terms)],
      cdiff(ten[splice(inds)], i) + terms)));