	(* color algebra: we assume all color structures have been
	   reduced down to SUNT chains, with the following conventions:
	   - SUNT[a, i, j] = T^a_{ij}
	   - SUNT[a, b, c, ..., i, j] = (T^a T^b T^c ...)_{ij}
	   - SUNT[a, ..., 0, 0] = Tr T^a T^b T^c ...
	   - SUNT[i, j] = delta_{ij}	*)

sunT[i_Symbol, i_] := SUNN

sunT[_, i_Symbol, i_] = 0

sunT[a_Integer, a_, i_Symbol, i_] = 1/2

sunT[_Integer, _Integer, i_Symbol, i_] = 0


	(* apply the Fierz identity for SU(N):
	   T^a_{ij} T^a_{kl} = 1/2 (delta_{il} delta_{kj} -
	                             1/N delta_{ij} delta_{kl})
	   in various guises *)

sunT[t1___, a_Symbol, t2___, a_, t3___, i_, j_] :=
  (sunT[t1, t3, i, j] sunTrace[t2] -
    sunT[t1, t2, t3, i, j]/SUNN)/2

sunT[t1___, a_Symbol, t2___, i_, j_] sunT[t3___, a_, t4___, k_, l_] ^:=
  (sunT[t1, t4, i, l] sunT[t3, t2, k, j] -
    sunT[t1, t2, i, j] sunT[t3, t4, k, l]/SUNN)/2

sunT[a___, i_, j_Symbol] sunT[b___, j_, k_] ^:= sunT[a, b, i, k]


	(* Mathematica does not automatically match a^2 = a*a *)

sunT/: sunT[a___, i_, j_Symbol]^2 :=
  Level[{{a}, Reverse[{a}], {i, i}}, {2}, sunT]

sunT/: sunT[a___, i_Symbol, j_]^2 :=
  Level[{Reverse[{a}], {a, j, j}}, {2}, sunT]


sunTrace[] := SUNN

sunTrace[_] = 0

sunTrace[a__] := sunT[a, #, #]&[ Unique["col"] ]


	(* we need sunText and sunTrace to generate dummy indices
	   for the summation implied in the trace *)

sunText[i_Symbol, i_] := Sqrt[SUNN]

sunText[a_Symbol, a_, 0, 0] := Sqrt[SUNN^2 - 1]/2

sunText[a___, 0, 0] := sunTrace[a]

sunText[other__] := sunT[other]


conj[t_Times] := conj/@ t

conj[t_SUNT] := RotateLeft[Reverse[t], 2]


ColorSimplify[amp_, ampc_] := ColorSimplify[amp conj[ampc]]

ColorSimplify[amp_] :=
  Simplify[ Expand[amp /. SUNT -> sunText] /. sunT -> SUNT ]


(* try it: *)

(*SUNN = 3*)

col1 = ColorSimplify[ SUNT[g, c1, c2] SUNT[g, c3, c4] ]

t1 = SUNT[Glu1, Glu2, Col2, Col1]

t2 = SUNT[Glu2, Glu1, Col2, Col1]

col2a = ColorSimplify[t1, t1]

col2b = ColorSimplify[t1, t2]

col2c = ColorSimplify[t2, t1]

col2d = ColorSimplify[t2, t2]

