(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Input::Initialization:: *)
xAct`SymSpin`$Version={"0.1.1",{2022,06,02}}


(* ::Input::Initialization:: *)
(* SymSpin: A symmetric spinor algebra *)

(* Copyright (C) 2021 - 2022 Steffen Aksteiner and Thomas B\[ADoubleDot]ckdahl *)

(* This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License,or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307, USA. 
*)


(* ::Input::Initialization:: *)
(* :Title: SymSpin *)

(* :Authors: Steffen Aksteiner and Thomas B\[ADoubleDot]ckdahl *)

(* :Summary: Setup of a symmetric spinor algebra *)

(* :Brief Discussion:
   - Introduces a product on symmetric spinor spaces. *)
  
(* :Context: xAct`SymSpin` *)

(* :Package Version: 0.1.1 *)

(* :Copyright: Steffen Aksteiner and Thomas B\[ADoubleDot]ckdahl (2021 - 2022) *)

(* :History: See SymSpin.History *)

(* :Keywords: *)

(* :Source: SymSpin.nb *)

(* :Warning: *)

(* :Mathematica Version: 9.0 and later *)

(* :Limitations: *)
	
(* :Acknowledgements: We thank Simon Jacobsson for testing and suggestions. *)


(* ::Input::Initialization:: *)
If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`SymSpin`"];


(* ::Input::Initialization:: *)
BeginPackage["xAct`SymSpin`",{"xAct`SpinFrames`","xAct`SymManipulator`","xAct`TexAct`","xAct`xCoba`","xAct`Spinors`","xAct`xTensor`","xAct`xPerm`","xAct`xCore`"}]


(* ::Input::Initialization:: *)
Print[xAct`xCore`Private`bars]
Print["Package xAct`SymSpin`  version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2021, Steffen Aksteiner and Thomas B\[ADoubleDot]ckdahl, under the General Public License."];


(* ::Input::Initialization:: *)
Off[General::shdw]
xAct`SymSpin`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS\.b4 WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]


(* ::Input::Initialization:: *)
If[xAct`xCore`Private`$LastPackage==="xAct`SymSpin`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]]


(* ::Input::Initialization:: *)
InitSymSpin::usage="InitSymSpin[sigma] initializes the SymSpin package by setting $DefaultSpinBundle to be the VBundleOfSolderingForm of sigma";
$DefaultSpinBundle::usage="$DefaultSpinBundle is a global variable initialized by InitSymSpin[sigma] and set to be the VBundleOfSolderingForm of sigma";
HeldxTensorQ::usage="Displays tensors using PrintAs even if no indices are given.";
SymSpinor::usage ="SymSpinor[numunprim,numprim,symb,vb] is a symmetric spinor with numunprim unprimed and numprim primed indices and symbol symb (optional, blank by default) on the vector bundle vb (optional, $DefaultSpinBundle by default). numunprim and numprim can be numeric or more generally, constants.";
DefSymmetricSpinor::usage ="DefSymmetricSpinor[symb,numunprim,numprim,vb,printas] defines a symmetric spinor symb with numunprim unprimed and numprim primed indices on the vector bundle vb and with displayform printas. numunprim and numprim can be numeric or more generally, constants.";
MultScal::usage ="MultScal[scal,TT] multiplies thr scsalar field scal to the tensor head TT.";
DummyTensor::usage="DummyTensor[vbs_,sym_] is a tensor head element of the vector bundle vbs (e.g. by SlotsOfTensor) and symmetry sym (e.g. by SymmetryGroupOfTensor)";
TensorPlus::usage="TensorPlus[T1, T2, ...] is a Plus on tensor heads Ti.";
IsolateTensorPlus::usage="IsolateTensorPlus[eq,i] isolates the ith term of a TensorPlus equation with default value i=1.";
IsolateTensorPlusPattern::usage="IsolateTensorPlusPattern[eq,p] isolates the term containing p from an equation with TensorPlus.";
SetLinearOperatorRules::usage="SetLinearOperatorRules[op] sets op to be linear w.r.t. TensorPlus and MultScal.";
LinearOperatorQ::usage="LinearOperatorQ[op] gives True if op was set by SetLinearOperatorRules[op] ad False otherwise.";
CommuteOp::usage="CommuteOp[op1, op2] commutes operators op1 and op2. Examples are fundamental spinor operators or CommuteOp[SymMult,MultScal]";
SymMult::usage="SymMult[T1,numunprim,numprim,vb][T2] contracts numunprim unprimed and numprim primed indices (on the vector bundle vb, which is optional and $DefaultSpinBundle by default) of symmetric spinors T1, T2 and symmetrizes.";
SymMultRule::usage="Expands SymMult[T1,numunprim,numprim,vb][T2][inds] into indexed spinor expressions.";
SymHToSymMultRule::usage="Converts SymH terms involving two symmetric spinors into SymMult terms.";
SymMultLeibnizRules::usage="SymMultLeibnizRules[covd] commutes SymMult out of Fundamental spinor operators of covd (optional, First[SpinCovDsOfSolderingForm[SolderingFormOfVBundle[$DefaultSpinBundle]]] by default).";
MultScalLeibnizRules::usage="MultScalLeibnizRules[covd] commutes MultScal out of Fundamental spinor operators of covd (optional, First[SpinCovDsOfSolderingForm[SolderingFormOfVBundle[$DefaultSpinBundle]]] by default).";
SortSymMult::usage="SortSymMult[orderingfunc] sorts T1, T2 in SymMult[T1,m,n,spin][T2] according to orderingfunc.";
SortSymMultReverse::usage="Reverse ordering to SortSymMult";
SymMultToMultScalRule::usage="Converts SymMult[f,0,0,spin_][T] into MultScal[f[],T] if f is a scalar function.";
MultScalToSymMultRule::usage="MultScalToSymMultRule[spin] converts MultScal[f[],T]\[RuleDelayed]SymMult[f,0,0,spin][T]";
IrrDecomposeSymMult::usage="IrrDecomposeSymMult[T1,T2,{p_,q_},spin_] expands (p,q) contractions of T1, T2 into a TensorPlus of SymMult terms.";
SymMultCoeffFunc::usage="SymMultCoeffFunc[i_,r_,k_,t_,mm_,MM_] is the coefficient function appearing in sorting nested SymMult objects.";
CommuteSymMultRuleOut::usage="CommuteSymMultRuleOut[T2] commutes SymMult[T2,..] outside in SymMult[T1,t,u,spin][SymMult[T2,m,n,spin][T3]]";
CommuteSymMultRuleIn::usage="CommuteSymMultRuleIn[T1] commutes SymMult[T1,..] inside in SymMult[T1,t,u,spin][SymMult[T2,m,n,spin][T3]]";
ToGHPComponents::usage="ToGHPComponents[expr,dyad] returns a list of dyad components of the index free expr";
TermToIndexFree::usage="TermToIndexFree[term] is a function called from ToIndexFree that turns a single term into an index free expression. If special code is needed for some specific situation please add the definitions here.";
ToIndexFree::usage="ToIndexFree[expr] converts an indexed expression expr to an unindexed version. It is assumed that expr does not have any dummy indices and that it is completely symmetric.";
ToIndexed::usage="ToIndexed[expr] assumes that expr is an index free expression or an equation with indexfree expressions. The output is the same with appropriate indices inserted.";
TensorMinus::usage ="TensorMinus is automatically converted into the corresponding TensorPlus expression.";
IndexFreeEqToZeroLHS::usage ="IndexFreeEqToZeroLHS[eq] transforms the index free tensor equation eq to a form with a zero tensor in the left hand side.";
ExtractCoeffsIndexFree::usage="ExtractCoeffsIndexFree[expr, field] contracts expr with a test field and extracts the coefficients of all irreducible components of field times the test field. This is still experimental so please check the validity afterwards. expr can be an indexfree spinor, an index free equation or a list. field can indexfree spinor or a list.";
IsolateFirstInLHS::usage="IsolateFirstInLHS[eq] isolates the first term in the left hand side of an indexfree equation.";
IsolateFirstInRHS::usage="IsolateFirstInRHS[eq] isolates the first term in the right hand side of an indexfree equation.";


(* ::Input::Initialization:: *)
Begin["`Private`"]


(* ::Input::Initialization:: *)
$ContextPath


(* ::Input::Initialization:: *)
InitSymSpin[sigma_]:=$DefaultSpinBundle=VBundleOfSolderingForm[sigma]


(* ::Input::Initialization:: *)
SetAttributes[HeldxTensorQ,HoldAllComplete];
HeldxTensorQ[expr_]:=xTensorQ[Unevaluated[expr]];
MakeBoxes[tensor_?HeldxTensorQ,StandardForm]:=xAct`xTensor`Private`interpretbox[tensor,xAct`xTensor`Private`xTensorBox[tensor,{}]]


(* ::Input::Initialization:: *)
SymmetricSpinorOfArbitraryValenceQ[]:=SymmetricSpinorOfArbitraryValenceQ[$DefaultSpinBundle];


(* ::Input::Initialization:: *)
PrintAs[DummySymbol]^:="";
Tex[DummySymbol]^:="";
MakeDaggerSymbol[DummySymbol]^:=DummySymbol;


(* ::Input::Initialization:: *)
SymSpinor[numunprim_,numprim_,vb_?VBundleQ]:=SymSpinor[numunprim,numprim,DummySymbol,vb]
SymSpinor[numunprim_,numprim_,symb_]:=SymSpinor[numunprim,numprim,symb,$DefaultSpinBundle]
SymSpinor[numunprim_,numprim_]:=SymSpinor[numunprim,numprim,DummySymbol,$DefaultSpinBundle]


(* ::Input::Initialization:: *)
xTensorQ[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=True;
SlotsOfTensor[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=SymmetricOfValence[numunprim,numprim,vb];
SymmetryGroupOfTensor[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=If[And[IntegerQ@numunprim,IntegerQ@numprim,numunprim>=0,numprim>=0],xAct`SymManipulator`Private`CompatibleSym@SlotsOfTensor[SymSpinor[numunprim,numprim,symb,vb]],Throw@Message[SymmetryGroupOfTensor::error,"Asked for the symmetry group of spinor with valence numbers("<>ToString[numunprim]<>","<>ToString[numprim]<>")"]];
SymmetricSpinorOfArbitraryValenceQ[vb_][SymSpinor[numunprim_,numprim_,symb_,vb_]]^=True;
PrintAs[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=PrintAs[symb];
Tex[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=Tex[symb];
Dagger[SymSpinor[numunprim_,numprim_,symb_,vb_]]^:=SymSpinor[numprim,numunprim,MakeDaggerSymbol[symb],vb]


(* ::Input::Initialization:: *)
DefSymmetricSpinor[symb_,numunprim_,numprim_,vb_,printas_]:=With[{symbdg=MakeDaggerSymbol@symb},
xTensorQ[symb]^:=True;
SlotsOfTensor[symb]^:=Simplify@SymmetricOfValence[numunprim,numprim,vb];
SymmetryGroupOfTensor[symb]^:=If[And[IntegerQ@numunprim,IntegerQ@numprim,numunprim>=0,numprim>=0],xAct`SymManipulator`Private`CompatibleSym@SlotsOfTensor[symb],Throw@Message[SymmetryGroupOfTensor::error,"Asked for the symmetry group of spinor with valence numbers("<>ToString[numunprim]<>","<>ToString[numprim]<>")"]];
SymmetricSpinorOfArbitraryValenceQ[vb][symb]^=True;
DefInfo[symb]^:={"spinor", ""};
HostsOf[symb]^:={BaseOfVBundle[vb],vb};
ServantsOf[symb]^:={symbdg};
DependenciesOfTensor[symb]^:={BaseOfVBundle[vb]};
PrintAs[symb]^:=printas;
Tex[symb]^:=Tex[printas];
Dagger[symb]^:=symbdg;
xTensorQ[symbdg]^:=True;
SlotsOfTensor[symbdg]^:=SymmetricOfValence[numprim,numunprim,vb];
SymmetryGroupOfTensor[symbdg]^:=If[And[IntegerQ@numunprim,IntegerQ@numprim,numunprim>=0,numprim>=0],xAct`SymManipulator`Private`CompatibleSym@SlotsOfTensor[symbdg],Throw@Message[SymmetryGroupOfTensor::error,"Asked for the symmetry group of spinor with valence numbers("<>ToString[numprim]<>","<>ToString[numunprim]<>")"]];
SymmetricSpinorOfArbitraryValenceQ[vb][symbdg]^=True;
DefInfo[symbdg]^:={"spinor", ""};
HostsOf[symbdg]^:={BaseOfVBundle[Dagger[vb]],Dagger[vb]};
ServantsOf[symbdg]^:={};
DependenciesOfTensor[symbdg]^:={BaseOfVBundle[Dagger[vb]]};
PrintAs[symbdg]^:=OverscriptBox[printas,"_"];
Tex[symbdg]^:=Tex[OverscriptBox[printas,"_"]];
Dagger[symbdg]^:=symb;]


(* ::Input::Initialization:: *)
xTensorQ[MultScal[scal_,LL_?xTensorQ]]^=True;
SlotsOfTensor[MultScal[scal_,LL_?xTensorQ]]^:=SlotsOfTensor[LL];
SymmetryGroupOfTensor[MultScal[scal_,LL_?xTensorQ]]^:=SymmetryGroupOfTensor[LL];


(* ::Input::Initialization:: *)
MultScal[c_,expr_TensorPlus]:=MultScal[c,#]&/@expr;
MultScal[scal_,ZeroTensor[vbs_]]:=ZeroTensor[vbs];
MultScal[c1_,MultScal[c2_,TT_?xTensorQ]]:=MultScal[c1*c2,TT];
MultScal[1,expr_]:=expr;
MultScal[scal_,LL_][inds___]:=scal*LL[inds];


(* ::Input::Initialization:: *)
NumOfUnprimedSlots[spin_][MultScal[Scal_,BB_?xTensorQ]]^:=NumOfUnprimedSlots[spin][BB];
NumOfPrimedSlots[spin_][MultScal[Scal_,BB_?xTensorQ]]^:=NumOfPrimedSlots[spin][BB];


(* ::Input::Initialization:: *)
MultScal[0,TT_?xTensorQ]:=ZeroTensorOfTensor@TT;


(* ::Input::Initialization:: *)
Dagger[MultScal[scal_,LL_?xTensorQ]]^:=MultScal[Dagger@scal,DaggerTensor@LL]


(* ::Input::Initialization:: *)
PrintAs[MultScal[scal_Plus,LL_?xTensorQ]]^:=RowBox[{"(","(",MakeBoxes[scal,StandardForm],")",PrintAs[LL],")"}];
PrintAs[MultScal[-1,LL_?xTensorQ]]^:=RowBox[{"(","-",PrintAs[LL],")"}];PrintAs[MultScal[scal_,LL_?xTensorQ]]^:=RowBox[{"(",MakeBoxes[scal,StandardForm],PrintAs[LL],")"}];


(* ::Input::Initialization:: *)
PrintAsNoPar[MultScal[scal_Plus,LL_?xTensorQ]]^:=RowBox[{"(",MakeBoxes[scal,StandardForm],")",PrintAs[LL]}];
PrintAsNoPar[MultScal[-1,LL_?xTensorQ]]^:=RowBox[{"-",PrintAs[LL]}];
PrintAsNoPar[MultScal[scal_,LL_?xTensorQ]]^:=RowBox[{MakeBoxes[scal,StandardForm],PrintAs[LL]}];
PrintAsNoPar[expr___]:=PrintAs[expr]


(* ::Input::Initialization:: *)
Tex[MultScal[scal_,LL_?xTensorQ]]^:=StringJoin[xAct`TexAct`Private`TexOpen["("],xAct`TexAct`Private`TexFactor[scal/.xAct`TexAct`ToOrderedPlus],xAct`SymManipulator`Private`TexFundOp[LL],xAct`TexAct`Private`TexClose[")"]];


(* ::Input::Initialization:: *)
xAct`SymManipulator`Private`TexFundOp[MultScal[scal_,LL_?xTensorQ]]^:=StringJoin[xAct`TexAct`Private`TexFactor[scal/.xAct`TexAct`ToOrderedPlus],xAct`SymManipulator`Private`TexFundOp[LL]];


(* ::Input::Initialization:: *)
$KopExclude={};


(* ::Input::Initialization:: *)
ImposeMultScal[expr_Plus,scal_]:=ImposeMultScal[#,scal]&/@expr;
ImposeMultScal[expr_Equal,scal_]:=ImposeMultScal[#,scal]&/@expr;
ImposeMultScal[0,scal_]:=0;
ImposeMultScal[expr_,0]:=0;
ImposeMultScal[c_?ConstantQ*expr_,scal_]:=c*ImposeMultScal[expr,scal];
ImposeMultScal[TT_?xTensorQ[tinds___],scal_]:=MultScal[scal,TT][tinds];
ImposeMultScal[expr_Times,scal_]:=With[{factorlist=xAct`xTensor`Private`ListOfFactors@expr},Times@@Select[factorlist,Function[x,Not@AllTrue[$KopExclude,FreeQ[x,#]&]]]*ImposeMultScal2[Times@@Select[factorlist,Function[x,AllTrue[$KopExclude,FreeQ[x,#]&]]],scal]]
ImposeMultScal2[TT_?xTensorQ[tinds___],scal_]:=MultScal[scal,TT][tinds];
ImposeMultScal2[Times[TT_?xTensorQ[tinds1___],z__],scal_]:=Times[MultScal[scal,TT][tinds1],z];
ImposeMultScal2[z_,scal_]:=Times[z,scal]


(* ::Input::Initialization:: *)
xTensorQ[TensorPlus[TT_?xTensorQ,more__]]^:=True;
SlotsOfTensor[TensorPlus[TT_?xTensorQ,more__]]^:=SlotsOfTensor@TT;


(* ::Input::Initialization:: *)
TensorPlus[TT_?xTensorQ]:=TT;


(* ::Input::Initialization:: *)
TensorPlus[tensors___][inds___]:=Plus@@(#[inds]&/@List[tensors])


(* ::Input::Initialization:: *)
TensorPlus[z___,ZeroTensor[vbs_],more___]:=TensorPlus[z,more];
TensorPlus[z___,0,more___]:=TensorPlus[z,more];
TensorPlus[x___,TensorPlus[expr__],z___]:=TensorPlus[x,expr,z];


(* ::Input::Initialization:: *)
TensorPlus[z___,a_,b___,a_,d___]:=TensorPlus[z,MultScal[2,a],b,d];
TensorPlus[z___,a_,f___,MultScal[c_,a_],d___]:=TensorPlus[z,MultScal[c+1,a],f,d];
TensorPlus[z___,MultScal[c_,a_],f___,a_,d___]:=TensorPlus[z,MultScal[c+1,a],f,d];
TensorPlus[z___,MultScal[b_,a_],f___,MultScal[c_,a_],d___]:=TensorPlus[z,MultScal[b+c,a],f,d];


(* ::Input::Initialization:: *)
TensorPlus[lhs1_==rhs1_,lhs2_==rhs2_]:=TensorPlus[lhs1,lhs2]==TensorPlus[rhs1,rhs2];
TensorPlus[lhs1_==rhs1_,expr_]:=TensorPlus[lhs1,expr]==TensorPlus[rhs1,expr];


(* ::Input::Initialization:: *)
TensorMinus[lhs1_==rhs1_,lhs2_==rhs2_]:=TensorPlus[lhs1,MultScal[-1,lhs2]]==TensorPlus[rhs1,MultScal[-1,rhs2]];
TensorMinus[lhs1_==rhs1_,expr_]:=TensorPlus[lhs1,MultScal[-1,expr]]==TensorPlus[rhs1,MultScal[-1,expr]];
TensorMinus[expr1_,expr2_]:=TensorPlus[expr1,MultScal[-1,expr2]];


(* ::Input::Initialization:: *)
NumOfUnprimedSlots[spin_][TensorPlus[TT_?xTensorQ,more__]]^:=NumOfUnprimedSlots[spin]@TT;
NumOfPrimedSlots[spin_][TensorPlus[TT_?xTensorQ,more__]]^:=NumOfPrimedSlots[spin]@TT;


(* ::Input::Initialization:: *)
CoeffFromMultScal[MultScal[c_?ConstantQ, T_]] := c;
CoeffFromMultScal[_] := 1;


(* ::Input::Initialization:: *)
TensorPlusToList[expr_TensorPlus]:=List@@expr;
TensorPlusToList[ZeroTensor[__]]:={};
TensorPlusToList[expr_]:={expr};


(* ::Input::Initialization:: *)
IsolateTensorPlus[eq_Equal,i_:1]:=Module[{tmpexpr=TensorPlusToList@TensorMinus[eq[[1]],eq[[2]]]},
MultScal[1/CoeffFromMultScal[tmpexpr[[i]]],tmpexpr[[i]]]==MultScal[-CoeffFromMultScal[tmpexpr[[i]]]^(-1),TensorPlus@@Append[Delete[tmpexpr,{i}],ZeroTensorOfTensor@tmpexpr[[i]]]]]


(* ::Input::Initialization:: *)
IsolateTensorPlusPattern[eq_Equal,pat_]:=Module[{tmpexpr=TensorPlusToList@TensorMinus[eq[[1]],eq[[2]]],pos,i},
pos=Position[tmpexpr,pat];
If[Length@pos==0,Throw@Message[IsolateTensorPlusPattern::error, "Could not find "<>ToString@pat]];
i=First@First@pos;
MultScal[1/CoeffFromMultScal[tmpexpr[[i]]],tmpexpr[[i]]]==MultScal[-1/CoeffFromMultScal[tmpexpr[[i]]],TensorPlus@@Append[Delete[tmpexpr,{i}],ZeroTensorOfTensor@tmpexpr[[i]]]]];


(* ::Input::Initialization:: *)
IsolateFirstInLHS[Equal[lhs_,rhs_]]:=IsolateTensorPlusPattern[lhs==rhs,RemoveMultscalCoeff@First@TensorPlusToList[lhs]];


(* ::Input::Initialization:: *)
IsolateFirstInRHS[Equal[lhs_,rhs_]]:=IsolateTensorPlusPattern[lhs==rhs,RemoveMultscalCoeff@First@TensorPlusToList[rhs]];


(* ::Input::Initialization:: *)
PrintAs[TensorPlus[expr1_,expr2___]]^:=RowBox[Join[{xAct`xTensor`Private`interpretbox[expr1,PrintAsNoPar[expr1]]},Join@@(If[xAct`SpinFrames`Private`BoxStartWithMinusQ[#[[2]]],{xAct`xTensor`Private`interpretbox@@#},{"+",xAct`xTensor`Private`interpretbox@@#}]&/@({#,PrintAsNoPar[#]}&/@List[expr2]))]]


(* ::Input::Initialization:: *)
Tex[TensorPlus[expr___]]^:=StringJoin[Riffle[xAct`SymManipulator`Private`TexFundOp/@List[expr]," + "]];


(* ::Input::Initialization:: *)
xAct`TexAct`Private`TexFactor[expr_TensorPlus]:=StringJoin[xAct`TexAct`Private`TexOpen["("],Tex[expr],xAct`TexAct`Private`TexClose[")"]];


(* ::Input::Initialization:: *)
xAct`SymManipulator`Private`TexFundOp[expr_TensorPlus]:=xAct`TexAct`Private`TexFactor[expr];


(* ::Input::Initialization:: *)
SymmetricSpinorOfValenceQ[TensorPlus[arg__],k_,l_,spin_]^:=And@@(SymmetricSpinorOfValenceQ[#,k,l,spin]&/@List[arg]);


(* ::Input::Initialization:: *)
xTensorQ[DummyTensor[vbs_,sym_]]^=True;
SlotsOfTensor[DummyTensor[vbs_,sym_]]^:=vbs;
SymmetryGroupOfTensor[DummyTensor[vbs_,sym_]]^:=sym;
PrintAs[DummyTensor[vbs_,sym_]]^="";
Tex[DummyTensor[vbs_,sym_]]^="";
TensorToDummyTensor[TT_?xTensorQ]:=TT->DummyTensor[SlotsOfTensor[TT],SymmetryGroupOfTensor[TT]]


(* ::Input::Initialization:: *)
Dagger[DummyTensor[vbs_,sym_]]^:=DummyTensor[Reverse[Dagger/@vbs],xAct`SymManipulator`Private`CompatibleSym@Reverse[Dagger/@vbs]]/;CompatibleSymQ[vbs,sym]


(* ::Input::Initialization:: *)
Dagger[DummyTensor[vbs_,sym_][inds___]]^:=PlaceIndicesInTensor[Evaluate@Dagger[DummyTensor[vbs,sym]],DaggerIndex/@{inds}]


(* ::Input::Initialization:: *)
SetLinearOperatorRules[op_]:=Module[{},
LinearOperatorQ[op]^=True;
op[expr_TensorPlus,z___]:=op[#,z]&/@expr;
op[ZeroTensor[vbs_],z___]:=ZeroTensorOfTensor@op[DummyTensor[vbs,CompatibleSymmetric[vbs]],z];
op[MultScal[c_?ConstantQ,TT_?xTensorQ],z___]:=MultScal[c,op[TT,z]];
]


(* ::Input::Initialization:: *)
DaggerTensor[expr_]:=Which[HermitianQ[expr],expr,AntihermitianQ[expr],MultScal[-1,expr],True,Dagger[expr]]


Dagger[TensorPlus[arg___]]^:=DaggerTensor/@TensorPlus[arg]


(* ::Input::Initialization:: *)
ToTensorPlusRules={Plus[T1_?TensorQ,T2_?TensorQ]:>TensorPlus[T1,T2],c_?ConstantQ*T_?TensorQ:>MultScal[c,T]};


(* ::Input::Initialization:: *)
TensorPlusSortFunc[MultScal[a_,b_]]:=TensorPlusSortFunc[b];
TensorPlusSortFunc[b_]:=b;


(* ::Input::Initialization:: *)
SortTensorPlusRule=expr_TensorPlus:>SortBy[expr,TensorPlusSortFunc];


(* ::Input::Initialization:: *)
SymmetricDummySpinor[k_,l_,spin_:$DefaultSpinBundle]:=With[{vbs=Join[ConstantArray[-spin,k],ConstantArray[-Dagger[spin],l]]},DummyTensor[vbs,xAct`SymManipulator`Private`CompatibleSym@vbs]]


ExprToIndexFreeTensor[0]:=0;
ExprToIndexFreeTensor[expr_Plus]:=TensorPlus@@(ExprToIndexFreeTensor/@expr);
ExprToIndexFreeTensor[c_?ConstantQ*expr_]:=MultScal[c,ExprToIndexFreeTensor[expr]];
ExprToIndexFreeTensor[TT_?xTensorQ[]*expr_]:=MultScal[TT[],ExprToIndexFreeTensor[expr]];
ExprToIndexFreeTensor[TT_?xTensorQ[inds___]]:=TT;


ExprToIndexFreeTensor2[lia_List,vbs_]:=ExprToIndexFreeTensor[#,vbs]&/@lia


(* ::Input::Initialization:: *)
IndexFreeEqToZeroLHS[Equal[lhs_?xTensorQ,rhs_?xTensorQ]]:=ZeroTensorOfTensor[lhs]==TensorMinus[rhs,lhs]


(* ::Input::Initialization:: *)
xTension["SymSpin`",DefFundSpinOperators,"End"]=ExtraDefFundSpinOpFunction;


(* ::Input::Initialization:: *)
ExtraDefFundSpinOpFunction[covd_,options___]:=With[{divcd=xAct`SymManipulator`Private`DivName@covd, curlcd=xAct`SymManipulator`Private`CurlName@covd, curldgcd=xAct`SymManipulator`Private`CurlDgName@covd, twistcd=xAct`SymManipulator`Private`TwistName@covd,
spin=VBundleOfSolderingForm[SolderingFormOfSpinCovD[covd]]},
SetLinearOperatorRules@divcd;
SetLinearOperatorRules@curlcd;
SetLinearOperatorRules@curldgcd;
SetLinearOperatorRules@twistcd;
CommuteOp[divcd,curlcd]^=divcd[curlcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[k/(k+1),curlcd[divcd[TT]]],
MultScal[(l-2),SymMult[Dagger[Psi[covd]],0,3,spin][TT]],
MultScal[k,SymMult[Phi[covd],1,2,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,2,spin];
CommuteOp[curlcd,divcd]^=curlcd[divcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[(k+1)/k,divcd[curlcd[TT]]],
MultScal[-(l-2)(k+1)/k,SymMult[Dagger[Psi[covd]],0,3,spin][TT]],
MultScal[-(k+1),SymMult[Phi[covd],1,2,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,2,spin];
CommuteOp[divcd,curldgcd]^=divcd[curldgcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[l/(l+1),curldgcd[divcd[TT]]],
MultScal[(k-2),SymMult[Psi[covd],3,0,spin][TT]],
MultScal[l,SymMult[Phi[covd],2,1,spin][TT]]
]]/;EnoughIndsQSpin[TT,2,0,spin];
CommuteOp[curldgcd,divcd]^=curldgcd[divcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[(l+1)/l,divcd[curldgcd[TT]]],
MultScal[-(k-2)(l+1)/l,SymMult[Psi[covd],3,0,spin][TT]],
MultScal[-(l+1),SymMult[Phi[covd],2,1,spin][TT]]
]]/;EnoughIndsQSpin[TT,2,1,spin];
CommuteOp[curlcd,twistcd]^=curlcd[twistcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[l/(l+1),twistcd[curlcd[TT]]],
MultScal[k,SymMult[Psi[covd],1,0,spin][TT]],
MultScal[l,SymMult[Phi[covd],0,1,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,0,spin];
CommuteOp[twistcd,curlcd]^=twistcd[curlcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[(l+1)/l,curlcd[twistcd[TT]]],
MultScal[-k(l+1)/l,SymMult[Psi[covd],1,0,spin][TT]],
MultScal[-(l+1),SymMult[Phi[covd],0,1,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,1,spin];
CommuteOp[curldgcd,twistcd]^=curldgcd[twistcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[k/(k+1),twistcd[curldgcd[TT]]],
MultScal[l,SymMult[Dagger[Psi[covd]],0,1,spin][TT]],
MultScal[k,SymMult[Phi[covd],1,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,0,spin];
CommuteOp[twistcd,curldgcd]^=twistcd[curldgcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[(k+1)/k,curldgcd[twistcd[TT]]],
MultScal[-l(k+1)/k,SymMult[Dagger[Psi[covd]],0,1,spin][TT]],
MultScal[-(k+1),SymMult[Phi[covd],1,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,0,spin];
CommuteOp[divcd,twistcd,curlcd,curldgcd]^=divcd[twistcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[-(1/(k+1)+1/(l+1)),curlcd[curldgcd[TT]]],
MultScal[l*(l+2)/(l+1)^2,twistcd[divcd[TT]]],
MultScal[(l+2)(k-1)/(l+1),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[(l+2)l/(l+1),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(l+2)(k+2)/(l+1),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[ l (l-1)/(l+1),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[ l k/(l+1),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[- l(l+2)/(l+1),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,0,spin];
CommuteOp[twistcd,divcd,curlcd,curldgcd]^=twistcd[divcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[((1+l)*(2+k+l))/((1+k)*l*(2+l)),curlcd[curldgcd[TT]]],
MultScal[(1+l)^2/(l*(2+l)),divcd[twistcd[TT]]],
MultScal[-(k-1) (1+l)/l,SymMult[Psi[covd],2,0,spin][TT]],
MultScal[-l(1+l)/l,SymMult[Phi[covd],1,1,spin][TT]],
MultScal[(k+2)(1+l)/l,SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[-l(l-1)(1+l)/(l*(2+l)),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[-k(1+l)/(2+l),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[l+1,SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,1,spin];
CommuteOp[divcd,twistcd,curldgcd,curlcd]^=divcd[twistcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[-(1/(k+1)+1/(l+1)),curldgcd[curlcd[TT]]],
MultScal[k*(k+2)/(k+1)^2,twistcd[divcd[TT]]],
MultScal[(k+2)(l-1)/(k+1),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[(k+2)k/(k+1),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(k+2)(l+2)/(k+1),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[k (k-1)/(k+1),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[l k/(k+1),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-k(k+2)/(k+1),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,1,spin];
CommuteOp[twistcd,divcd,curldgcd,curlcd]^=twistcd[divcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
MultScal[((1+k)*(2+k+l))/((1+l)*k*(2+k)),curldgcd[curlcd[TT]]],
MultScal[(1+k)^2/(k*(2+k)),divcd[twistcd[TT]]],
MultScal[-(l-1)(1+k)/k,SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[-(1+k),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[(l+2)(1+k)/k,SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[-(k-1)(1+k)/(2+k),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[-l (1+k)/(2+k),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[(1+k),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,1,spin];
CommuteOp[curlcd,curldgcd,curldgcd,curlcd]^=curlcd[curldgcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
curldgcd[curlcd[TT]],
MultScal[(1/(k+1)-1/(l+1)),twistcd[divcd[TT]]],
MultScal[(k-1),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[l,SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(k+2),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[-(l-1),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[ -k,SymMult[Phi[covd],1,1,spin][TT]],
MultScal[(l+2),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,1,spin];
CommuteOp[curldgcd,curlcd,curlcd,curldgcd]^=curldgcd[curlcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[
curlcd[curldgcd[TT]],
MultScal[(1/(l+1)-1/(k+1)),twistcd[divcd[TT]]],
MultScal[-(k-1),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[-l,SymMult[Phi[covd],1,1,spin][TT]],
MultScal[(k+2),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[(l-1),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[k,SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(l+2),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,1,spin];
CommuteOp[curldgcd,curlcd,divcd,twistcd]^=curldgcd[curlcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[-(((1+k) (1+l))/(2+k+l)),divcd[twistcd[TT]]],
MultScal[(k (2+k) (1+l))/((1+k) (2+k+l)),twistcd[divcd[TT]]],
MultScal[((2+k) (-1+l) (1+l))/(2+k+l),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[(k (2+k) (1+l))/(2+k+l),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(((2+k) (1+l) (2+l))/(2+k+l)),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[((-1+k) k (1+l))/(2+k+l),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[(k l (1+l))/(2+k+l),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-((k (2+k) (1+l))/(2+k+l)),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,0,1,spin];
CommuteOp[curlcd,curldgcd,divcd,twistcd]^=curlcd[curldgcd[TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[-(((1+k) (1+l))/(2+k+l)),divcd[twistcd[TT]]],
MultScal[((1+k) l (2+l))/((1+l) (2+k+l)),twistcd[divcd[TT]]],
MultScal[((-1+k) (1+k) (2+l))/(2+k+l),SymMult[Psi[covd],2,0,spin][TT]],
MultScal[((1+k) l (2+l))/(2+k+l),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(((1+k) (2+k) (2+l))/(2+k+l)),SymMult[Lambda[covd],0,0,spin][TT]],
MultScal[((1+k) (-1+l) l)/(2+k+l),SymMult[Dagger[Psi[covd]],0,2,spin][TT]],
MultScal[(k (1+k) l)/(2+k+l),SymMult[Phi[covd],1,1,spin][TT]],
MultScal[-(((1+k) l (2+l))/(2+k+l)),SymMult[Lambda[covd],0,0,spin][TT]]
]]/;EnoughIndsQSpin[TT,1,0,spin];]


(* ::Input::Initialization:: *)
TermToIndexFree[TT_?xTensorQ[inds___]]:=If[SymmetricSpinorOfArbitraryValenceQ[$DefaultSpinBundle]@TT,TT,Throw@Message[TermToIndexFree::error,"The tensor "<>ToString[TT,InputForm]<>" is not symmetric."]];
TermToIndexFree[c_?ConstantQ*z_]:=MultScal[c,TermToIndexFree[z]];
TermToIndexFree[c_?xTensorQ[]*z_]:=MultScal[c[],TermToIndexFree[z]];
TermToIndexFree[(c_?xTensorQ[])^(i_)]:=MultScal[c[]^(i-1),c];
TermToIndexFree[(c_?xTensorQ[])^(i_)*z_]:=MultScal[c[]^i,TermToIndexFree[z]];


(* ::Input::Initialization:: *)
TermToIndexFree[expr_]:=Throw@Message[TermToIndexFree::error,"Can not turn "<>ToString[expr,InputForm]<>" into an unindexed form"]


(* ::Input::Initialization:: *)
ToIndexFree[l_List]:=ToIndexFree/@l;
ToIndexFree[lhs_==0]:=With[{indexfreelhs=ToIndexFree[lhs]},indexfreelhs==ZeroTensorOfTensor@indexfreelhs];
ToIndexFree[0==rhs_]:=With[{indexfreerhs=ToIndexFree[rhs]},ZeroTensorOfTensor@indexfreerhs==indexfreerhs];
ToIndexFree[eq_Equal]:=ToIndexFree/@eq;
ToIndexFree[expr_]:=TensorPlus@@(TermToIndexFree/@xAct`xTensor`Private`ListOfTerms[expr]);


ToIndexed[eq_Equal]:=ToIndexed/@eq;
ToIndexed[TT_?xTensorQ]:=GiveIndicesToTensor[TT]//.SymMultRule


(* ::Input::Initialization:: *)
SymMult[AA_?xTensorQ,unprim_,prim_]:=SymMult[AA,unprim,prim,$DefaultSpinBundle]


(* ::Input::Initialization:: *)
xTensorQ[SymMult[AA_?xTensorQ,unprim_,prim_,vb_][BB_?xTensorQ]]^=True;


(* ::Input::Initialization:: *)
SlotsOfTensor[SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]^:=
With[{k=NumOfUnprimedSlots[spin][AA]+NumOfUnprimedSlots[spin][BB]-2unprim,l=NumOfPrimedSlots[spin][AA]+NumOfPrimedSlots[spin][BB]-2prim},If[And[xAct`SymManipulator`Private`$SlotsOfTensorWarning,Or[k<0,l<0]],Throw@Message[SymMult::error,StringJoin["Not enough indices:",ToString@InputForm@SymMult[ AA,unprim,prim,spin][BB]]],SymmetricOfValence[k,l,spin,{}]]];


(* ::Input::Initialization:: *)
SymmetryGroupOfTensor[SymMult[AA_?xTensorQ,unprim_,prim_,vb_][BB_?xTensorQ]]^:=xAct`SymManipulator`Private`CompatibleSym@SlotsOfTensor[SymMult[AA,unprim,prim,vb][BB]];


(* ::Input::Initialization:: *)
SymMultRule=SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ][inds___]:>Module[{unprimdummies=DummyIn/@Table[spin,unprim],primdummies=DummyIn/@Table[Dagger@spin,prim],unprimedinds=Select[{inds},SymbolJoin[spin,"`pmQ"]],primedinds=Select[{inds},SymbolJoin[Dagger@spin,"`pmQ"]],numunprimA=NumOfUnprimedVBundles[spin][SlotsOfTensor[AA]],numprimA=NumOfPrimedVBundles[spin][SlotsOfTensor[AA]],numunprimB=NumOfUnprimedVBundles[spin][SlotsOfTensor[BB]],numprimB=NumOfPrimedVBundles[spin][SlotsOfTensor[BB]]},
ImposeSym[xAct`SymManipulator`Private`PlaceIndicesInSpinor[AA,{Sequence@@unprimdummies,Sequence@@(Take[unprimedinds,numunprimA-unprim]),Sequence@@primdummies,Sequence@@(Take[primedinds,numprimA-prim])},{spin,Dagger@spin}]xAct`SymManipulator`Private`PlaceIndicesInSpinor[BB,{Sequence@@(ChangeIndex/@unprimdummies),Sequence@@(Take[unprimedinds,-(numunprimB-unprim)]),Sequence@@(ChangeIndex/@primdummies),Sequence@@(Take[primedinds,-(numprimB-prim)])},{spin,Dagger@spin}],IndexList[inds]]];


(* ::Input::Initialization:: *)
PrintAs[SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]^:=RowBox[{If[Not@AtomQ[AA],"(",""],PrintAs[AA],If[Not@AtomQ[AA],")",""],UnderoverscriptBox["\[CircleDot]",If[OptionValue[DefFundSpinOperators,ShowValenceInfo],StringJoin[ToString[DisplayForm@NumOfUnprimedSlots[spin]@BB,StandardForm],",",ToString[DisplayForm@NumOfPrimedSlots[spin]@BB,StandardForm]],""],StringJoin[ToString[DisplayForm@unprim,StandardForm],",",ToString[DisplayForm@prim,StandardForm]]],If[Not@AtomQ[BB],"(",""],PrintAs[BB],If[Not@AtomQ[BB],")",""]}];


(* ::Input::Initialization:: *)
NumOfUnprimedSlots[spin_][SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]^:=NumOfUnprimedSlots[spin][AA]+NumOfUnprimedSlots[spin][BB]-2unprim;
NumOfPrimedSlots[spin_][SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]^:=NumOfPrimedSlots[spin][AA]+NumOfPrimedSlots[spin][BB]-2prim;


(* ::Input::Initialization:: *)
SymMult[expr_TensorPlus,unprim_,prim_,spin_][BB_]:=(SymMult[#1,unprim,prim,spin][BB]&)/@expr;
SymMult[AA_,unprim_,prim_,spin_][expr_TensorPlus]:=(SymMult[AA,unprim,prim,spin][#1]&)/@expr;
SymMult[ZeroTensor[vbs_],unprim_,prim_,spin_][BB_]:=ZeroTensorOfTensor[SymMult[DummyTensor[vbs,CompatibleSymmetric[vbs]],unprim,prim,spin][BB]];
SymMult[AA_,unprim_,prim_,spin_][ZeroTensor[vbs_]]:=ZeroTensorOfTensor[SymMult[AA,unprim,prim,spin]@DummyTensor[vbs,CompatibleSymmetric[vbs]]];
SymMult[MultScal[c_?ConstantQ,AA_?xTensorQ],unprim_,prim_,spin_][BB_]:=MultScal[c,SymMult[AA,unprim,prim,spin][BB]];
SymMult[AA_,unprim_,prim_,spin_][MultScal[c_?ConstantQ,BB_?xTensorQ]]:=MultScal[c,SymMult[AA,unprim,prim,spin][BB]];


(* ::Input::Initialization:: *)
CommuteOp[SymMult,MultScal]={SymMult[MultScal[s_,AA_],unprim_,prim_,spin_][BB_]:>MultScal[s,SymMult[AA,unprim,prim,spin][BB]],SymMult[AA_,unprim_,prim_,spin_][MultScal[s_,BB_]]:>MultScal[s,SymMult[AA,unprim,prim,spin][BB]]};


(* ::Input::Initialization:: *)
Dagger[SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]^:=SymMult[If[HermitianQ[AA],AA,Dagger[AA]],prim,unprim,spin][If[HermitianQ[BB],BB,Dagger[BB]]]


(* ::Input::Initialization:: *)
Dagger[SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ][inds__]]^:=xAct`SymManipulator`Private`PlaceIndicesInSpinor[Dagger[SymMult[AA,unprim,prim,spin][BB]],DaggerIndex/@{inds},{spin,Dagger@spin}]


(* ::Input::Initialization:: *)
Tex[x:(SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ])]:=StringJoin[xAct`TexAct`Private`TexOpen["("],xAct`SymManipulator`Private`TexFundOp@x,xAct`TexAct`Private`TexClose[")"]];


(* ::Input::Initialization:: *)
xAct`SymManipulator`Private`TexFundOp[SymMult[AA_?xTensorQ,unprim_,prim_,spin_][BB_?xTensorQ]]:=StringJoin[xAct`TexAct`Private`TexFactor[AA],If[OptionValue[DefFundSpinOperators,ShowValenceInfo],StringJoin["\\underset{",Tex[NumOfUnprimedSlots[spin][BB]/.ToOrderedPlus],",",Tex[NumOfPrimedSlots[spin][BB]/.ToOrderedPlus],"}"],""],"{\\overset{",StringJoin[Tex[unprim/.ToOrderedPlus],",",Tex[prim/.ToOrderedPlus]],"}{\\odot }}",xAct`SymManipulator`Private`TexFundOp[BB]];


(* ::Input::Initialization:: *)
SpinorOfValenceQ[k_,l_,spin_][TT_]:=And[xTensorQ[TT],NumOfUnprimedSlots[spin][TT]==k,NumOfPrimedSlots[spin][TT]==l]


(* ::Input::Initialization:: *)
CommuteOp[SymMult[AA_,1,1,spin_],SymMult[AA_,0,1,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[AA,0,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
If[k>0,MultScal[k/(k+1),SymMult[AA,0,1,spin][SymMult[AA,1,1,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,1,1,spin][SymMult[AA,0,1,spin][TT]]]]/;EnoughIndsQSpin[TT,0,2,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,1,spin_],SymMult[AA_,1,1,spin_]]^:=(SymMult[AA,0,1,spin][SymMult[AA,1,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
MultScal[(k+1)/k,SymMult[AA,1,1,spin][SymMult[AA,0,1,spin][TT]]]]/;EnoughIndsQSpin[TT,1,2,spin])/;SpinorOfValenceQ[1,1,spin][AA];CommuteOp[SymMult[AA_,1,1,spin_],SymMult[AA_,1,0,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[AA,1,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
If[l>0,MultScal[l/(l+1),SymMult[AA,1,0,spin][SymMult[AA,1,1,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,1,1,spin][SymMult[AA,1,0,spin][TT]]]]/;EnoughIndsQSpin[TT,2,0,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,1,0,spin_],SymMult[AA_,1,1,spin_]]^:=(SymMult[AA,1,0,spin][SymMult[AA,1,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
MultScal[(l+1)/l,SymMult[AA,1,1,spin][SymMult[AA,1,0,spin][TT]]]]/;EnoughIndsQSpin[TT,2,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,1,spin_],SymMult[AA_,0,0,spin_]]^:=(SymMult[AA,0,1,spin][SymMult[AA,0,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
If[l>0,MultScal[l/(l+1),SymMult[AA,0,0,spin][SymMult[AA,0,1,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,0,1,spin][SymMult[AA,0,0,spin][TT]]]])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,1,0,spin_],SymMult[AA_,0,0,spin_]]^:=(SymMult[AA,1,0,spin][SymMult[AA,0,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
If[k>0,MultScal[k/(k+1),SymMult[AA,0,0,spin][SymMult[AA,1,0,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,1,0,spin][SymMult[AA,0,0,spin][TT]]]])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,0,spin_],SymMult[AA_,0,1,spin_]]^:=(SymMult[AA,0,0,spin][SymMult[AA,0,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
MultScal[(l+1)/l,SymMult[AA,0,1,spin][SymMult[AA,0,0,spin][TT]]]]/;EnoughIndsQSpin[TT,0,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,0,spin_],SymMult[AA_,1,0,spin_]]^:=(SymMult[AA,0,0,spin][SymMult[AA,1,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
MultScal[(k+1)/k,SymMult[AA,1,0,spin][SymMult[AA,0,0,spin][TT]]]]/;EnoughIndsQSpin[TT,1,0,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,1,spin_],SymMult[AA_,1,0,spin_],SymMult[AA_,1,0,spin_],SymMult[AA_,0,1],spin_]^:=(SymMult[AA,0,1,spin][SymMult[AA,1,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[SymMult[AA,1,0,spin][SymMult[AA,0,1,spin][TT]],MultScal[(1/(k+1)-1/(l+1)),SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT]]]]]/;EnoughIndsQSpin[TT,1,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,1,0,spin_],SymMult[AA_,0,1,spin_],SymMult[AA_,0,1,spin_],SymMult[AA_,1,0,spin_]]^:=(SymMult[AA,1,0,spin][SymMult[AA,0,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[SymMult[AA,0,1,spin][SymMult[AA,1,0,spin][TT]],MultScal[(1/(l+1)-1/(k+1)),SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT]]]]]/;EnoughIndsQSpin[TT,1,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,1,1,spin_],SymMult[AA_,0,0,spin_],SymMult[AA_,0,1,spin_],SymMult[AA_,1,0,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[-(1/(k+1)+1/(l+1)),SymMult[AA,0,1,spin][SymMult[AA,1,0,spin][TT]]],If[l>0,MultScal[l*(l+2)/(l+1)^2,SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT]]]]]/;EnoughIndsQSpin[TT,1,0,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,1,1,spin_],SymMult[AA_,0,0,spin_],SymMult[AA_,1,0,spin_],SymMult[AA_,0,1,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[-(1/(k+1)+1/(l+1)),SymMult[AA,1,0,spin][SymMult[AA,0,1,spin][TT]]],If[k>0,MultScal[k*(k+2)/(k+1)^2,SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT]]],ZeroTensorOfTensor@SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT]]]]]/;EnoughIndsQSpin[TT,0,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,0,spin_],SymMult[AA_,1,1,spin_],SymMult[AA_,0,1,spin_],SymMult[AA_,1,0,spin_]]^:=(SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[((1+l)*(2+k+l))/((1+k)*l*(2+l)),SymMult[AA,0,1,spin][SymMult[AA,1,0,spin][TT]]],MultScal[(1+l)^2/(l*(2+l)),SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT]]]]]/;EnoughIndsQSpin[TT,1,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];
CommuteOp[SymMult[AA_,0,0,spin_],SymMult[AA_,1,1,spin_],SymMult[AA_,1,0,spin_],SymMult[AA_,0,1,spin_]]^:=(SymMult[AA,0,0,spin][SymMult[AA,1,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},
TensorPlus[MultScal[((1+k)*(2+k+l))/((1+l)*k*(2+k)),SymMult[AA,1,0,spin][SymMult[AA,0,1,spin][TT]]],MultScal[(1+k)^2/(k*(2+k)),SymMult[AA,1,1,spin][SymMult[AA,0,0,spin][TT]]]]]/;EnoughIndsQSpin[TT,1,1,spin])/;SpinorOfValenceQ[1,1,spin][AA];


(* ::Input::Initialization:: *)
AllSymMultCommutatorRules[AA_,spin:$DefaultSpinBundle]:={CommuteOp[SymMult[AA,1,1,spin],SymMult[AA,0,1,spin]],CommuteOp[SymMult[AA,1,1,spin],SymMult[AA,1,0,spin]],CommuteOp[SymMult[AA,0,1,spin],SymMult[AA,0,0,spin]],CommuteOp[SymMult[AA,1,0,spin],SymMult[AA,0,0,spin]],CommuteOp[SymMult[AA,1,1,spin],SymMult[AA,0,0,spin],SymMult[AA,1,0,spin],SymMult[AA,0,1,spin]],CommuteOp[SymMult[AA,1,1,spin],SymMult[AA,0,0,spin],SymMult[AA,0,1,spin],SymMult[AA,1,0,spin]],CommuteOp[SymMult[AA,0,1,spin],SymMult[AA,1,0,spin],SymMult[AA,1,0,spin],SymMult[AA,0,1,spin]]}/;SpinorOfValenceQ[1,1,spin][AA];


(* ::Input::Initialization:: *)
CommuteOp[SymMult[AA_,i_,j_,spin_],SymMult[BB_,i_,j_,spin_]]^:=(SymMult[AA,i,j,spin][SymMult[BB,i,j,spin][TT_?xTensorQ]]:>SymMult[BB,i,j,spin][SymMult[AA,i,j,spin][TT]]/;AA=!=BB)/;SpinorOfValenceQ[1,1,spin][AA]&&SpinorOfValenceQ[1,1,spin][BB];


(* ::Input::Initialization:: *)
CommuteOp[SymMult[AA_,1,1,spin_],SymMult[BB_,0,1,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[BB,0,1,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},TensorPlus[MultScal[-1/(k+1),SymMult[BB,1,1,spin][SymMult[AA,0,1,spin][TT]]],MultScal[k (k+2)/(k+1)^2,SymMult[BB,0,1,spin][SymMult[AA,1,1,spin][TT]]]]]/;AA=!=BB)/;SpinorOfValenceQ[1,1,spin][AA]&&SpinorOfValenceQ[1,1,spin][BB];
CommuteOp[SymMult[AA_,1,1,spin_],SymMult[BB_,1,0,spin_]]^:=(SymMult[AA,1,1,spin][SymMult[BB,1,0,spin][TT_?xTensorQ]]:>With[{k=NumOfUnprimedSlots[spin]@TT,l=NumOfPrimedSlots[spin]@TT},TensorPlus[MultScal[-1/(l+1),SymMult[BB,1,1,spin][SymMult[AA,1,0,spin][TT]]],MultScal[l (l+2)/(l+1)^2,SymMult[BB,1,0,spin][SymMult[AA,1,1,spin][TT]]]]]/;AA=!=BB)/;SpinorOfValenceQ[1,1,spin][AA]&&SpinorOfValenceQ[1,1,spin][BB];


(* ::Input::Initialization:: *)
SymHToSymMultRule={x:SymH[{T1_,T2_ },group_,label_][inds___]:>Module[{dummybundles=VBundleOfIndex/@List@@FindDummyIndices@Evaluate@x,symmultexpr,symmulteq,TFirst,TSecond,spin=$DefaultSpinBundle},
{TFirst,TSecond}=SortBy[{T1,T2},Length@SlotsOfTensor@#&];
symmultexpr=GiveIndicesToTensor@SymMult[TFirst,NumOfUnprimedVBundles[spin]@dummybundles,NumOfPrimedVBundles[spin]@dummybundles][TSecond];
symmulteq=symmultexpr==ToCanonicalSym[symmultexpr/.SymMultRule];
ToCanonicalSym[x]/.EqToCompareRule@FlipEquation@symmulteq
],x:(T1_?(And[xTensorQ[#],Not@MetricQ[#]&&Not[#===delta]]&)[indsa__]T2_?(And[xTensorQ[#],Not@MetricQ[#]&&Not[#===delta]]&)[indsb__]):> Module[{dummybundles=VBundleOfIndex/@List@@FindDummyIndices@Evaluate@x,unprimed,primed,symmultexpr,symmulteq,TFirst,TSecond,spin=$DefaultSpinBundle},
unprimed=NumOfUnprimedVBundles[spin]@dummybundles;
primed=NumOfPrimedVBundles[spin]@dummybundles;
{TFirst,TSecond}=SortBy[{T1,T2},Length@SlotsOfTensor@#&];
If[unprimed+primed==0,x,
symmultexpr=GiveIndicesToTensor@SymMult[TFirst,unprimed,primed][TSecond];
symmulteq=symmultexpr==ToCanonical[symmultexpr/.SymMultRule];
If[symmulteq===True,x,x/.EqToCompareRule@FlipEquation@symmulteq]
]]};


(* ::Input::Initialization:: *)
SymMultLeibnizRules[]:=SymMultLeibnizRules[First[SpinCovDsOfSolderingForm[SolderingFormOfVBundle[$DefaultSpinBundle]]]]


(* ::Input::Initialization:: *)
SymMultLeibnizRules[covd_]:=
With[{
divcd=xAct`SymManipulator`Private`DivName@covd, 
curlcd=xAct`SymManipulator`Private`CurlName@covd, 
curldgcd=xAct`SymManipulator`Private`CurlDgName@covd, 
twistcd=xAct`SymManipulator`Private`TwistName@covd},
{twistcd[SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>
With[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi]},
TensorPlus[
MultScal[(-1)^(m+n),SymMult[\[CurlyPhi],m,n,spin]@twistcd[\[Phi]]],
MultScal[(-1)^(m+n)*n/(j+1),SymMult[\[CurlyPhi],m,n-1,spin]@curlcd[\[Phi]]],
MultScal[(-1)^(m+n)*m/(i+1),SymMult[\[CurlyPhi],m-1,n,spin]@curldgcd[\[Phi]]],
MultScal[(-1)^(m+n)*m*n/(i+1)/(j+1),SymMult[\[CurlyPhi],m-1,n-1,spin]@divcd[\[Phi]]],
SymMult[\[Phi],m,n,spin]@twistcd[\[CurlyPhi]],
MultScal[n/(l+1),SymMult[\[Phi],m,n-1,spin]@curlcd[\[CurlyPhi]]],
MultScal[m/(k+1),SymMult[\[Phi],m-1,n,spin]@curldgcd[\[CurlyPhi]]],
MultScal[m*n/(k+1)/(l+1),SymMult[\[Phi],m-1,n-1,spin]@divcd[\[CurlyPhi]]]
]],
curlcd[SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>
With[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi]},
TensorPlus[
MultScal[-(-1)^(m+n)(l-n)/(j+l-2n),SymMult[\[CurlyPhi],m,n+1,spin]@twistcd[\[Phi]]],
MultScal[(-1)^(m+n)*(j-n)(j+l-n+1)/(j+l-2n)/(j+1),SymMult[\[CurlyPhi],m,n,spin]@curlcd[\[Phi]]],
MultScal[-(-1)^(m+n)*m*(l-n)/(j+l-2n)/(i+1),SymMult[\[CurlyPhi],m-1,n+1,spin]@curldgcd[\[Phi]]],
MultScal[(-1)^(m+n)*m*(j-n)(j+l-n+1)/(j+l-2n)/(i+1)/(j+1),SymMult[\[CurlyPhi],m-1,n,spin]@divcd[\[Phi]]],
MultScal[-(j-n)/(j+l-2n),SymMult[\[Phi],m,n+1,spin]@twistcd[\[CurlyPhi]]],
MultScal[(l-n)(j+l-n+1)/(j+l-2n)/(l+1),SymMult[\[Phi],m,n,spin]@curlcd[\[CurlyPhi]]],
MultScal[-m*(j-n)/(j+l-2n)/(k+1),SymMult[\[Phi],m-1,n+1,spin]@curldgcd[\[CurlyPhi]]],
MultScal[m*(l-n)(j+l-n+1)/(j+l-2n)/(k+1)/(l+1),SymMult[\[Phi],m-1,n,spin]@divcd[\[CurlyPhi]]]
]],
curldgcd[SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>
With[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi]},
TensorPlus[
MultScal[-(-1)^(m+n)(k-m)/(i+k-2m),SymMult[\[CurlyPhi],m+1,n,spin]@twistcd[\[Phi]]],
MultScal[-(-1)^(m+n)*n*(k-m)/(i+k-2m)/(j+1),SymMult[\[CurlyPhi],m+1,n-1,spin]@curlcd[\[Phi]]],
MultScal[(-1)^(m+n)*(i-m)(i+k-m+1)/(i+k-2m)/(i+1),SymMult[\[CurlyPhi],m,n,spin]@curldgcd[\[Phi]]],
MultScal[(-1)^(m+n)*n*(i-m)(i+k-m+1)/(i+k-2m)/(i+1)/(j+1),SymMult[\[CurlyPhi],m,n-1,spin]@divcd[\[Phi]]],
MultScal[-(i-m)/(i+k-2m),SymMult[\[Phi],m+1,n,spin]@twistcd[\[CurlyPhi]]],
MultScal[-n*(i-m)/(i+k-2m)/(l+1),SymMult[\[Phi],m+1,n-1,spin]@curlcd[\[CurlyPhi]]],
MultScal[(k-m)(i+k-m+1)/(i+k-2m)/(k+1),SymMult[\[Phi],m,n,spin]@curldgcd[\[CurlyPhi]]],
MultScal[n*(k-m)(i+k-m+1)/(i+k-2m)/(k+1)/(l+1),SymMult[\[Phi],m,n-1,spin]@divcd[\[CurlyPhi]]]
]],
divcd[SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>
With[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi]},
TensorPlus[
MultScal[(-1)^(m+n)(k-m)(l-n)/(i+k-2m)/(j+l-2n),SymMult[\[CurlyPhi],m+1,n+1,spin]@twistcd[\[Phi]]],
MultScal[-(-1)^(m+n)(k-m)(j-n)(j+l-n+1)/(i+k-2m)/(j+l-2n)/(j+1),SymMult[\[CurlyPhi],m+1,n,spin]@curlcd[\[Phi]]],
MultScal[-(-1)^(m+n)(l-n)(i-m)(i+k-m+1)/(i+k-2m)/(j+l-2n)/(i+1),SymMult[\[CurlyPhi],m,n+1,spin]@curldgcd[\[Phi]]],
MultScal[(-1)^(m+n)(i-m)(j-n)(i+k-m+1)(j+l-n+1)/(i+k-2m)/(j+l-2n)/(i+1)/(j+1),SymMult[\[CurlyPhi],m,n,spin]@divcd[\[Phi]]],
MultScal[(i-m)(j-n)/(i+k-2m)/(j+l-2n),SymMult[\[Phi],m+1,n+1,spin]@twistcd[\[CurlyPhi]]],
MultScal[-(i-m)(l-n)(j+l-n+1)/(i+k-2m)/(j+l-2n)/(l+1),SymMult[\[Phi],m+1,n,spin]@curlcd[\[CurlyPhi]]],
MultScal[-(k-m)(j-n)(i+k-m+1)/(i+k-2m)/(j+l-2n)/(k+1),SymMult[\[Phi],m,n+1,spin]@curldgcd[\[CurlyPhi]]],
MultScal[(k-m)(l-n)(i+k-m+1)(j+l-n+1)/(i+k-2m)/(j+l-2n)/(k+1)/(l+1),SymMult[\[Phi],m,n,spin]@divcd[\[CurlyPhi]]]
]]}];


(* ::Input::Initialization:: *)
TwistMultScal[covd_,(T_?xTensorQ)[]]:=(xAct`SymManipulator`Private`TwistName@covd)[T]
TwistMultScal[covd_,Times[a_,rest__]]:=TensorPlus[MultScal[Times[rest],TwistMultScal[covd,a]],MultScal[a,TwistMultScal[covd,Times[rest]]]]
TwistMultScal[covd_,Power[a_,c_?ConstantQ]]:=MultScal[c*Power[a,c-1],TwistMultScal[covd,a]]
TwistMultScal[covd_,expr_Plus]:=TwistMultScal[covd,#]&/@TensorPlus@@expr
TwistMultScal[covd_,c_?ConstantQ]:=With[{spin=VBundlesOfCovD[covd][[2]],spindg=VBundlesOfCovD[covd][[3]]},ZeroTensor[{-spin,-spindg}]]


(* ::Input::Initialization:: *)
MultScalLeibnizRules[]:=MultScalLeibnizRules[First[SpinCovDsOfSolderingForm[SolderingFormOfVBundle[$DefaultSpinBundle]]]]


(* ::Input::Initialization:: *)
MultScalLeibnizRules[covd_]:=
With[{
divcd=xAct`SymManipulator`Private`DivName@covd, 
curlcd=xAct`SymManipulator`Private`CurlName@covd, 
curldgcd=xAct`SymManipulator`Private`CurlDgName@covd, 
twistcd=xAct`SymManipulator`Private`TwistName@covd,
spin=VBundleOfSolderingForm@SolderingFormOfSpinCovD[covd]},
{divcd[MultScal[expr_,T_]]:>TensorPlus[MultScal[expr,divcd[T]],SymMult[TwistMultScal[covd,expr],1,1,spin][T]],
curlcd[MultScal[expr_,T_]]:>TensorPlus[MultScal[expr,curlcd[T]],SymMult[TwistMultScal[covd,expr],0,1,spin][T]],
curldgcd[MultScal[expr_,T_]]:>TensorPlus[MultScal[expr,curldgcd[T]],SymMult[TwistMultScal[covd,expr],1,0,spin][T]],
twistcd[MultScal[expr_,T_]]:>TensorPlus[MultScal[expr,twistcd[T]],SymMult[TwistMultScal[covd,expr],0,0,spin][T]]}];


(* ::Input::Initialization:: *)
DefaultSortSymMultFunc[spin_]:={NumOfUnprimedSlots[spin][#]+NumOfPrimedSlots[spin][#],NumOfUnprimedSlots[spin][#],NumOfPrimedSlots[spin][#]}&


SortSymMult[] := SortSymMult[DefaultSortSymMultFunc[$DefaultSpinBundle]];
SortSymMultReverse[]:=SortSymMultReverse[DefaultSortSymMultFunc[$DefaultSpinBundle]];


(* ::Input::Initialization:: *)
SortSymMult[orderingfunc_]:={SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]:>MultScal[(-1)^(m+n),SymMult[\[CurlyPhi],m,n,spin][\[Phi]]]/;Not[OrderedQ[orderingfunc/@{\[Phi],\[CurlyPhi]}]],SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[Phi]_?xTensorQ]:>ZeroTensorOfTensor@SymMult[\[Phi],m,n,spin][\[Phi]]/;((-1)^(m+n)==-1)}


(* ::Input::Initialization:: *)
SortSymMultReverse[orderingfunc_]:={SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]:>MultScal[(-1)^(m+n),SymMult[\[CurlyPhi],m,n,spin][\[Phi]]]/;Not[OrderedQ[orderingfunc/@{\[CurlyPhi],\[Phi]}]],SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[Phi]_?xTensorQ]:>ZeroTensorOfTensor@SymMult[\[Phi],m,n,spin][\[Phi]]/;((-1)^(m+n)==-1)}


(* ::Input::Initialization:: *)
SymMultToMultScalRule=SymMult[\[Phi]_?xTensorQ,0,0,spin_][\[CurlyPhi]_?xTensorQ]:>MultScal[\[Phi][],\[CurlyPhi]]/;NumOfUnprimedSlots[spin][\[Phi]]+NumOfPrimedSlots[spin][\[Phi]]==0;


(* ::Input::Initialization:: *)
MultScalToSymMultRule[spin_]:={MultScal[(\[Phi]_?xTensorQ)[],\[CurlyPhi]_?xTensorQ]:>SymMult[\[Phi],0,0,spin][\[CurlyPhi]],MultScal[c_*(\[Phi]_?xTensorQ)[],\[CurlyPhi]_?xTensorQ]:>MultScal[c,SymMult[\[Phi],0,0,spin][\[CurlyPhi]]],
MultScal[Plus[(\[Phi]_?xTensorQ)[],d__],\[CurlyPhi]_?xTensorQ]:>TensorPlus[SymMult[\[Phi],0,0,spin][\[CurlyPhi]],MultScal[d,\[CurlyPhi]]],
MultScal[Plus[c_*(\[Phi]_?xTensorQ)[],d__],\[CurlyPhi]_?xTensorQ]:>TensorPlus[MultScal[c,SymMult[\[Phi],0,0,spin][\[CurlyPhi]]],MultScal[d,\[CurlyPhi]]]};


(* ::Input::Initialization:: *)
ProductSymmetric[Ainds_,Binds_,ADginds_,BDginds_]:=Module[{VBundles,SGSlist},
(* Compute the strong generating sets for the separate symmetries *)
VBundles=Select[First@First@Position[IndexList[Sequence@@Ainds,Sequence@@Binds,Sequence@@ADginds,Sequence@@BDginds],#,1]&/@#&/@{Ainds,Binds,ADginds,BDginds},(Length[#]>=2)&];
If[Length@VBundles==1,SGSlist=Symmetric/@VBundles,
SGSlist=Symmetric[#,Cycles]&/@VBundles;
(* Join the symmetries *)
While[Length[SGSlist]>=2,
SGSlist[[2]]=JoinSGS[SGSlist[[1]],SGSlist[[2]]];
SGSlist=Delete[SGSlist,1];
];];
(* Return the remaining SGS or the trivial SGS if the list is empty. *)
If[Length[SGSlist]==1,First@SGSlist,StrongGenSet[{}, GenSet[]]]
]


(* ::Input::Initialization:: *)
IrrDecomposeSymMult[\[Phi]_?xTensorQ,\[CurlyPhi]_?xTensorQ,{p_,q_},spin_]:=Module[{k,l,i,j,allinds,Ainds,Binds,Cinds,ADginds,BDginds,CDginds},
{i,j}=ValenceNumbersOfSpinor[\[Phi],spin];
{k,l}=ValenceNumbersOfSpinor[\[CurlyPhi],spin];
allinds=List@@GiveIndicesToTensor@SymMult[\[Phi],0,0,spin][\[CurlyPhi]];
Ainds=Take[allinds,i-p];
Binds=Take[Drop[allinds,i-p],k-p];
Cinds=Take[Drop[allinds,i+k-2p],p];
ADginds=Take[Drop[allinds,i+k],j-q];
BDginds=Take[Drop[allinds,i+k+j-q],l-q];
CDginds=Take[Drop[allinds,i+k+j+l-2q],q];
xAct`SymManipulator`Private`PlaceIndicesInSpinor[\[Phi],Join[Ainds,ChangeIndex/@Cinds,ADginds,ChangeIndex/@CDginds],{spin,Dagger@spin}]xAct`SymManipulator`Private`PlaceIndicesInSpinor[\[CurlyPhi],Join[Binds,Cinds,BDginds,CDginds],{spin,Dagger@spin}]==(-1)^(p+q)Sum[(-1)^(m+n)Binomial[i-p,m-p]Binomial[k-p,m-p]Binomial[j-q,n-q]Binomial[l-q,n-q]/Binomial[i+k-m-p+1,m-p]/Binomial[j+l-n-q+1,n-q]*ImposeSym[Product[First[MetricsOfVBundle[spin]][Ainds[[r]],Binds[[r]]],{r,1,m-p}]*Product[Dagger[First[MetricsOfVBundle[spin]]][ADginds[[r]],BDginds[[r]]],{r,1,n-q}]*xAct`SymManipulator`Private`PlaceIndicesInSpinor[SymMult[\[Phi],m,n,spin][\[CurlyPhi]],Flatten@{Drop[Ainds,m-p],Drop[Binds,m-p],Drop[ADginds,n-q],Drop[BDginds,n-q]},{spin,Dagger@spin}],IndexList@@Flatten@{Ainds,Binds,ADginds,BDginds},ProductSymmetric[Ainds,Binds,ADginds,BDginds]],{m,p,Min[i,k]},{n,q,Min[j,l]}]
]/;And[CompatibleSymQ[SlotsOfTensor@\[Phi],SymmetryGroupOfTensor@\[Phi]],CompatibleSymQ[SlotsOfTensor@\[CurlyPhi],SymmetryGroupOfTensor@\[CurlyPhi]]]


IrrDecomposeSymMult[\[Phi]_, \[CurlyPhi]_, {p_, q_}] := IrrDecomposeSymMult[\[Phi], \[CurlyPhi], {p, q}, $DefaultSpinBundle];


(* ::Input::Initialization:: *)
SymMultCoeffFunc[i_Integer,r_Integer,k_Integer,t_Integer,mm_Integer,MM_Integer]:=Sum[Sum[((-1)^(-p+t+q)*Binomial[k-mm,p]*Binomial[mm,MM-p-q]*Binomial[k-mm-p,q]*Binomial[-mm+r,-p+t]*Binomial[i-t,MM-p-q]*Binomial[-p+t,q])/(Binomial[1+i+k-MM-p,MM-p]*Binomial[MM-p,q]*Binomial[k-2*mm+r,t]),{q,Max[0,MM-mm-p,MM+t-i-p],Min[MM-p,k-mm-p,t-p]}],{p,Max[0,mm-r+t],Min[k-mm,MM,t]}];


(* ::Input::Initialization:: *)
SymMultCoeffFunc[i_,r_,k_,t_,mm_,MM_]:=FullSimplify@Sum[Sum[((-1)^(-p+t+q)*Binomial[k-mm,p]*Binomial[mm,MM-p-q]*Binomial[k-mm-p,q]*Binomial[-mm+r,-p+t]*Binomial[i-t,MM-p-q]*Binomial[-p+t,q])/(Binomial[1+i+k-MM-p,MM-p]*Binomial[MM-p,q]*Binomial[k-2*mm+r,t]),{q,0,MM-p}],{p,0,MM}];


(* ::Input::Initialization:: *)
CommuteSymMultRuleOut[\[Omega]_?xTensorQ,AllowSame_:False]:={SymMult[\[Phi]_?xTensorQ,t_,u_,spin_][SymMult[\[Omega],m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>Module[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi],r=NumOfUnprimedSlots[spin]@\[Omega],s=NumOfPrimedSlots[spin]@\[Omega]},
TensorPlus@@Flatten@Table[MultScal[SymMultCoeffFunc@@(Simplify/@{i,r,k,t,m,M})*SymMultCoeffFunc@@(Simplify/@{j,s,l,u,n,NN}),SymMult[\[Omega],t+m-M,u+n-NN,spin][SymMult[\[Phi],M,NN,spin][\[CurlyPhi]]]],{M,Max[0,t+m-r],Min[i,k]},{NN,Max[0,u+n-s],Min[j,l]}]]/;(Or[AllowSame,\[Phi]=!=\[Omega]]),SymMult[SymMult[\[Omega],m_,n_,spin_][\[CurlyPhi]_?xTensorQ],t_,u_,spin_][\[Phi]_?xTensorQ]:>Module[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi],r=NumOfUnprimedSlots[spin]@\[Omega],s=NumOfPrimedSlots[spin]@\[Omega]},
TensorPlus@@Flatten@Table[MultScal[(-1)^(t+u)*SymMultCoeffFunc@@(Simplify/@{i,r,k,t,m,M})*SymMultCoeffFunc@@(Simplify/@{j,s,l,u,n,NN}),SymMult[\[Omega],t+m-M,u+n-NN,spin][SymMult[\[Phi],M,NN,spin][\[CurlyPhi]]]],{M,Max[0,t+m-r],Min[i,k]},{NN,Max[0,u+n-s],Min[j,l]}]]/;(Or[AllowSame,\[Phi]=!=\[Omega]])}


(* ::Input::Initialization:: *)
CommuteSymMultRuleIn[\[Phi]_?xTensorQ]:=SymMult[\[Phi],t_,u_,spin_][SymMult[\[Omega]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ]]:>Module[{i=NumOfUnprimedSlots[spin]@\[Phi],j=NumOfPrimedSlots[spin]@\[Phi],k=NumOfUnprimedSlots[spin]@\[CurlyPhi],l=NumOfPrimedSlots[spin]@\[CurlyPhi],r=NumOfUnprimedSlots[spin]@\[Omega],s=NumOfPrimedSlots[spin]@\[Omega]},
TensorPlus@@Flatten@Table[MultScal[SymMultCoeffFunc@@(Simplify/@{i,r,k,t,m,M})*SymMultCoeffFunc@@(Simplify/@{j,s,l,u,n,NN}),SymMult[\[Omega],t+m-M,u+n-NN,spin][SymMult[\[Phi],M,NN,spin][\[CurlyPhi]]]],{M,0,Min[i,k]},{NN,0,Min[j,l]}]]


(* ::Input::Initialization:: *)
GHPComponentRules[SymMult[\[Phi]_?xTensorQ,m_,n_,spin_][\[CurlyPhi]_?xTensorQ],dyad_]:=GHPComponentRules[SymMult[\[Phi],m,n,spin][\[CurlyPhi]],dyad]=Module[{s,r,t,u,p,q,i=NumOfUnprimedSlotsDyad[dyad]@\[Phi],j=NumOfPrimedSlotsDyad[dyad]@\[Phi],k=NumOfUnprimedSlotsDyad[dyad]@\[CurlyPhi],l=NumOfPrimedSlotsDyad[dyad]@\[CurlyPhi]},
xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[\[CurlyPhi],dyad];
xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[\[Phi],dyad];
xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[SymMult[\[Phi],m,n,spin][\[CurlyPhi]],dyad];
If[$DyadCalcInfo,Print["Calculating the ",dyad,"components of ",SymMult[\[Phi],m,n,spin][\[CurlyPhi]]," in GHP form. (Improved code)"]];
Expand[Expand[Flatten[Table[xAct`SpinFrames`DyadComponentByNumberSymmetric[SymMult[\[Phi],m,n,spin][\[CurlyPhi]],{s,t},dyad]->Sum[Sum[(-1)^p*Binomial[i,r]Binomial[i-r,p]Binomial[r,m-p]Binomial[k,s+m-r]Binomial[k-s-m+r,m-p]Binomial[s+m-r,p]/(Binomial[i,m]Binomial[k,m]Binomial[m,p]Binomial[i+k-2m,s]),{p,0,m}]Sum[(-1)^q*Binomial[j,u]Binomial[j-u,q]Binomial[u,n-q]Binomial[l,t+n-u]Binomial[l-t-n+u,n-q]Binomial[t+n-u,q]/(Binomial[j,n]Binomial[l,n]Binomial[n,q]Binomial[j+l-2n,t]),{q,0,n}]xAct`SpinFrames`DyadComponentByNumberSymmetric[\[Phi],{r,u},dyad]xAct`SpinFrames`DyadComponentByNumberSymmetric[\[CurlyPhi],{s+m-r,t+n-u},dyad],{r,s+m-k,s+m},{u,t+n-l,t+n}],{s,0,i+k-2m},{t,0,j+l-2n}]]//.$GHPExtraRules]/.GHPComponentRules[\[Phi],dyad]/.GHPComponentRules[\[CurlyPhi],dyad]//.$GHPExtraRules]];


(* ::Input::Initialization:: *)
GHPComponentRules[lhs:(MultScal[s_,\[Phi]_?xTensorQ]),dyad_]:=GHPComponentRules[lhs,dyad]=With[{k=NumOfUnprimedSlotsDyad[dyad]@lhs,l=NumOfPrimedSlotsDyad[dyad]@lhs,comps=xAct`SpinFrames`AllDyadComponentsSymmetric[lhs,dyad]},
xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[lhs,dyad];
Thread[comps->Expand[Expand[comps//.$GHPExtraRules]/.GHPComponentRules[\[Phi],dyad]//.$GHPExtraRules]]];


ToGHPComponents[expr_,dyad_]:=ToGHPComponents1[expr,dyad]


(* ::Input::Initialization:: *)
ToGHPComponents1[0==expr_,dyad_]:=(0==#&)/@ToGHPComponents1[expr,dyad]
ToGHPComponents1[expr_==0,dyad_]:=(#==0&)/@ToGHPComponents1[expr,dyad]
ToGHPComponents1[expr_Equal,dyad_]:=Thread[ToGHPComponents1[#,dyad]&/@expr]
ToGHPComponents1[expr_Plus,dyad_]:=ToGHPComponents1[#,dyad]&/@expr
ToGHPComponents1[expr_TensorPlus,dyad_]:=ToGHPComponents1[#,dyad]&/@Plus@@expr
ToGHPComponents1[Times[c_?ScalarQ,more__],dyad_]:=Expand[c*ToGHPComponents1[Times[more],dyad]]


(* ::Input::Initialization:: *)
ToGHPComponents1[\[Phi]_?xTensorQ,dyad_]:=(xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[\[Phi],dyad];xAct`SpinFrames`AllDyadComponentsSymmetric[\[Phi],dyad])/.GHPComponentRules[\[Phi],dyad]


(* ::Input::Initialization:: *)
ToGHPComponents1[\[Phi]_?xTensorQ[inds___],dyad_]:=(xAct`SpinFrames`SetGHPWeightAndFormatSymmetric[\[Phi],dyad];xAct`SpinFrames`AllDyadComponentsSymmetric[\[Phi],dyad])/.GHPComponentRules[\[Phi],dyad]


(* ::Input::Initialization:: *)
ExtractSymMultCoeff[expr_,spinor_,spin_:$DefaultSpinBundle]:=Module[{sortedexpr=expr//.Join[SortSymMult[#=!=spinor&],CommuteSymMultRuleOut[spinor]],extractedcoeffs},
extractedcoeffs=ExtractSymMultCoeffFromTerm[#,spinor,spin]&/@TensorPlusToList[sortedexpr];
TensorPlus@@SortBy[extractedcoeffs,TensorPlusSortFunc]]


(* ::Input::Initialization:: *)
ExtractSymMultCoeffFromTerm[SymMult[spinor_,__][expr_],spinor_,spin_]:=expr;
ExtractSymMultCoeffFromTerm[MultScal[a_,SymMult[spinor_,__][expr_]],spinor_,spin_]:=MultScal[a,expr];
ExtractSymMultCoeffFromTerm[expr_,spinor_,spin_]:=ZeroTensor@SymmetricOfValence[Sequence@@(ValenceNumbersOfSpinor[spinor,spin]-ValenceNumbersOfSpinor[expr,spin]),spin]/;FreeQ[expr,spinor]


(* ::Input::Initialization:: *)
RemoveMultscalCoeff[MultScal[a_,expr_]]:=expr;
RemoveMultscalCoeff[expr_]:=expr;


(* ::Input::Initialization:: *)
ExtractCoeffsIndexFree[expr_?xTensorQ,field_?xTensorQ,spin_:$DefaultSpinBundle]:=Module[{sortedexpr,exprlist,TestField,ExprValence,FieldValence,TestFieldFieldIrrDecParts,contractedExpr,extractedcoeffs},
(* Make sure that "field" is innermost.*)
sortedexpr=expr//.Join[SortSymMult[Not@FreeQ[#,field]&],{CommuteSymMultRuleIn[field]}];
exprlist=Select[TensorPlusToList[sortedexpr],Not@FreeQ[#,field]&];
ExprValence=ValenceNumbersOfSpinor[expr,spin];
FieldValence=ValenceNumbersOfSpinor[field,spin];
(* Constructing a test field to contract away all free indices in expr. *)
TestField=SymSpinor[Sequence@@ExprValence,"T",spin];
(* Constructing all irreducible products of "field" and the test field. *)
TestFieldFieldIrrDecParts=Flatten@Table[SymMult[TestField,i,j]@field,{i,0,Min[ExprValence[[1]],FieldValence[[1]]]},{j,0,Min[ExprValence[[2]],FieldValence[[2]]]}];
(* Contracting the expression with the test field to get a scalar expression. *)
contractedExpr=TensorPlus@@SymMult[TestField,Sequence@@ExprValence]/@exprlist;
(* Commute the thestfield inside. *)
contractedExpr=contractedExpr//.CommuteSymMultRuleIn[TestField];
(* Return a list of extracted coefficients. *)
ExtractSymMultCoeff[contractedExpr,#,spin]&/@TestFieldFieldIrrDecParts
]


(* ::Input::Initialization:: *)
ExtractCoeffsIndexFree[Equal[lhs_?xTensorQ,rhs_?xTensorQ],field_?xTensorQ]:=DeleteCases[(ZeroTensorOfTensor[#]==#)&/@ExtractCoeffsIndexFree[TensorMinus[rhs,lhs],field],True]


(* ::Input::Initialization:: *)
ExtractCoeffsIndexFree[l_List,f_]:=ExtractCoeffsIndexFree[#,f]&/@l


(* ::Input::Initialization:: *)
ExtractCoeffsIndexFree[expr_,l_List]:=ExtractCoeffsIndexFree[expr,#]&/@l


(* ::Input::Initialization:: *)
End[];
EndPackage[];



