(* Subroutine for program Ampcalculator *)

sestimebegin=SessionTime[];

If[$VersionNumber<4.,Tr[ject_]:=Sum[ject[[ii,ii]],{ii,1,First[Dimensions[ject]]}]]

AllParticles={(Subscript[ϕ, 1]),(Subscript[ϕ, 2]),(Subscript[ϕ, 3]),(Subscript[ϕ, 4]),(Subscript[ϕ, 5]),(Subscript[ϕ, 6]),(Subscript[ϕ, 7]),(Subscript[ϕ, 8]),γ,(Subscript[W, 1]),(Subscript[W, 2])};

erroroutput[ject_]:=Do[resnb=CreateWindow[WindowSize->{800,700}];NotebookWrite[resnb,ject];]

(* strangenesscheck[ject] and its subfunctions nummeson[ject] and indizes[ject] check if the weak part of the generating functional is needed. Apart from that, the program stops if the process is of order 2 or higher in (Subscript[G, F]). *)

strangenesscheck[ject_]:=Do[stlist={0,0,0,1,-1,1,-1,0};wblist={0,0,0,0,0,0,0,0,1,1};stp=Sum[stlist[[nummeson[ject][[ii]]]],{ii,Length[nummeson[ject]]}];wbp=Sum[wblist[[indizes[ject][[ii]]]],{ii,Length[indizes[ject]]}];If[wbp>1,erroroutput["Program terminated: more than one W boson is not allowed."];Exit[]];If[Abs[stp]>1,erroroutput["Program terminated: processes with |ΔS|>1 are not allowed."];Exit[]];If[wbp==1,stcheck=0,If[Abs[stp]==0,stcheck=0,stcheck=1]]]

nummeson[ject_]:=Table[intern[[ii,2]],{ii,Length[intern=DeleteCases[DeleteCases[ject,(Subscript[W, ii_])],γ]]}]

indizes[ject_]:=Table[intern[[ii,2]],{ii,Length[intern=DeleteCases[ReplaceAll[ject,{(Subscript[W, 1])->(Subscript[W, 9]),(Subscript[W, 2])->(Subscript[W, 10])}],γ]]}]

(* chargecheck[ject] stops the program if charge is not conserved. *)

chargecheck[ject_]:=Do[qlist={1,-1,0,1,-1,0,0,0,1,-1};charge=Sum[qlist[[indizes[ject][[ii]]]],{ii,Length[indizes[ject]]}];If[charge≠0,erroroutput["Program terminated: charge is not conserved."];Exit[]]]

(* particlecheck[ject] stops the program if the list ExternalParticles has been chosen incorrectly. *)

particlecheck[ject_]:=Do[If[¬MemberQ[AllParticles,ject[[ii]]],erroroutput["Program terminated: the list ExternalParticles has been chosen incorrectly."];Exit[]],{ii,Length[ject]}]

(* removemomenta[ject] removes the momentum arguments from the list of external particles. *)

removemomenta[ject_]:=Do[expar={};Do[expar=Append[expar,Part[ject,ii][[0]]],{ii,Length[ject]}]]

expar1=ExternalParticles;

(* expar2 is the list of the external particles in the form {(Subscript[ϕ, 1])[(Subscript[p, 1])],...}.  *)

expar2=ReplaceAll[expar1,{(SubPlus[π])->(Subscript[ϕ, 1]),(SubMinus[π])->(Subscript[ϕ, 2]),(Subscript[π, 0])->(Subscript[ϕ, 3]),(SubPlus[K])->(Subscript[ϕ, 4]),(SubMinus[K])->(Subscript[ϕ, 5]),(Subscript[K, 0])->(Subscript[ϕ, 6]),(Subscript[(OverBar[K]), 0])->(Subscript[ϕ, 7]),(Subscript[η, 8])->(Subscript[ϕ, 8]),(SubPlus[W])->(Subscript[W, 1]),(SubMinus[W])->(Subscript[W, 2])}];

(* removemomenta[expar2] generates the list expar that is of the form {(Subscript[ϕ, 1]),...}.  *)

removemomenta[expar2]

chargecheck[expar];

strangenesscheck[expar];

particlecheck[expar];

protitle="External particles and definition of scalar products:";

contitle1=ExternalParticles;

contitle2=Table[replist1[[ii]]->replist2[[ii]],{ii,Length[replist1]}];

sestimetitle="Runtime:  ";

sesunit= " hours";

sesstyle="Section";

output[ject_]:=Do[sestimeend=SessionTime[];sestime=N[((sestimeend-sestimebegin)/3600),4];resnb=CreateWindow[WindowSize->{800,700}];NotebookWrite[resnb,Cell[protitle,sesstyle]];NotebookWrite[resnb,Cell[RowBox[{ToBoxes[contitle1]}],"Input"]];NotebookWrite[resnb,Cell[RowBox[{ToBoxes[contitle2]}],"Input"]];NotebookWrite[resnb,Cell[title,sesstyle]];NotebookWrite[resnb,Cell[RowBox[{ToBoxes[result]}],"Input"]];NotebookWrite[resnb,Cell[TextData[{sestimetitle,ToString[sestime],sesunit}],"Input"]];]

If[Length[expar]==6,phi6=1,phi6=0];

(* If the lists replist1 and replist2 have not been defined by the user, they are set equal to empty lists. *)

If[Length[replist1]==0||Length[replist2]==0,replist1={};replist2={}]

SortsOfParticles=11;

(* In the next line the warning message Part::pspec is turned off. *)

Off[Part::pspec];Off[General::spell1];Off[General::spell];Off[Simplify::time];Off[Power::infy];Off[Coefficient::ivar]

(* conj[ject] performs a complex conjugation. *)

conj[ject_]:=ReplaceAll[ject,{Complex[m_,n_]->Complex[m,-n]}]

(* de[ject,var,ind] transforms the derivatives of the fields occurring in ject into special functions so that a vector index ind can be introduced. *)

de[ject_,var_,ind_]:=ReplaceAll[(∂_var)ject,{(∂_var)(Subscript[ϕ, cc_])[var]->(Subscript[Dϕ, cc])[var,ind],(Subscript[Dϕ, cc_]^((1,0)))[var,ind1_]->(Subscript[DDϕ, cc])[var,ind1,ind],Derivative[0,1][Subscript][_,_]->0 ,Derivative[1,0][A][var,ind1_]->DA[var,ind,ind1],Derivative[1,0,0][DA][var,ind1_,ind2_]->DDA[var,ind,ind1,ind2],Derivative[1,0][(Subscript[W, dd_])][var,ind1_]->(Subscript[DW, dd])[var,ind,ind1],Derivative[1,0,0][(Subscript[DW, dd_])][var,ind1_,ind2_]->(Subscript[DDW, dd])[var,ind,ind1,ind2]}]

simp[ject_]:=If[phi6==0,Simplify[ject],ject]

simp2[ject_]:=If[norecurs==0,Simplify[ject],ReplaceAll[If[stcheck==0,simpfast[ject],Coefficient[ject,(Subscript[G, 8])]*(Subscript[G, 8])+Coefficient[ject,(Subscript[G, 27])]*(Subscript[G, 27])],{(Subscript[F, π])->F}]]

(* getord[ject, ord] and getord[ject, ord, par] evaluate ject at order ord in f, which is the number of external particles (a photon counts as two particles), or in par, which is a variable. *)

getord[ject_,ord_]:=Do[If[ord<0,ordn=0,ordn=ord];Return[simp[ReplaceAll[D[ject,{f,ordn}],{f->0}]]*((f^ordn)/(ordn !))]]

getord[ject_,ord_,par_]:=simp[ReplaceAll[D[ject,{par,ord}],{par->0}]]*((par^ord)/(ord !))

(* getordns[ject, ord] is the same as getord[ject, ord] without Simplify. *)

getordns[ject_,ord_]:=ReplaceAll[D[ject,{f,ord}],{f->0}]*((f^ord)/(ord !))

getordns[ject_,ord_,par_]:=ReplaceAll[D[ject,{par,ord}],{par->0}]*((par^ord)/(ord !))

(* ex2ord[ject, ord] and ex2ord[ject, ord, par] supply ject up to order ord in f or in par. *)

ex2ord[ject_,ord_]:=Do[If[ord<0,ordn=0,ordn=ord];Return[Sum[getord[ject,ii],{ii,0,ordn}]]]

ex2ord[ject_,ord_,par_]:=Do[If[ord<0,ordn=0,ordn=ord];Return[Sum[getord[ject,ii,par],{ii,0,ordn}]]]

(* ex2ordns[ject, ord] is the same as ex2ord[ject, ord] without Simplify. *)

ex2ordns[ject_,ord_]:=Sum[getordns[ject,ii],{ii,0,ord}]

ex2ordns[ject_,ord_,par_]:=Do[If[ord<0,ordn=0,ordn=ord];Return[Sum[getordns[ject,ii,par],{ii,0,ordn}]]]

ex2ordns2[ject_,ord_,par_]:=Do[Do[zxc=0;Do[jj=-ii;If[Position[ject,Power[par,jj]]≠{},zxc=1],{ii,1,8}]];If[zxc==0,Return[ex2ordns[ject,ord,par]],tject=ReplaceAll[ject,{par->(1/par)}];Return[ReplaceAll[Sum[getordns[tject,ii,par],{ii,0,ord}],{par->(1/par)}]]]]

(* Tr2ord[ject, ord] evaluates the trace of ject up to order ord in f. *)

Tr2ord[ject_,ord_]:=ex2ord[Tr[ject],ord]

(* funcchange[ject] transforms the fields so that they can be used as arguments of a function. *)

funcchange[ject_]:=ReplaceAll[ject,{ϕ[1]->ϕ1,ϕ[2]->ϕ2,ϕ[3]->ϕ3,ϕ[4]->ϕ4,ϕ[5]->ϕ5,ϕ[6]->ϕ6,ϕ[7]->ϕ7,ϕ[8]->ϕ8,W[1]->W1,W[2]->W2}]

(* dephi[ject,i] is needed for derivatives with respect to (Subscript[ϕ, i]). *)

dephi[ject_,i_]:=ReplaceAll[ReplaceAll[ReplaceAll[Expand[D[ReplaceAll[ject,{(Subscript[Dϕ, cc_])[n_,m_]->ab[(Subscript[ϕ, cc])[n],n]}],(Subscript[ϕ, i])[var]]],{(tt_)*Derivative[1,0][ab][m_,n_]->-ab[tt,n]}],{ab->D}],{(Subscript[ϕ, cc_]^′)[var_]->(Subscript[Dϕ, cc])[var,(Subscript[α, var])],(Subscript[ϕ, cc_]^″)[var]->(Subscript[DDϕ, cc])[var,(Subscript[α, var]),(Subscript[α, var])],Derivative[1,0][A][var_,ind_]->DA[var,(Subscript[α, var]),ind],Derivative[1,0,0][DA][var_,inde_,ind_]->DA[var,inde,(Subscript[α, var]),ind],Derivative[1,0][(Subscript[W, dd_])][var_,ind_]->(Subscript[DW, dd])[var,(Subscript[α, var]),ind],Derivative[1,0,0][(Subscript[DW, dd_])][var_,inde_,ind_]->(Subscript[DDW, dd])[var,(Subscript[α, var]),inde,ind],Derivative[0,1][Subscript][_,_]->0}]

(* wkorr[ject] extracts the squared momentum of a W boson. *)

wkorr[ject_]:=ReplaceAll[ject,{(Subscript[DDW, cc_]^((1,0,0,0)))[var_,ind1_,ind2_,ind2_]->-sp[wmo,wmo](Subscript[DW, cc])[var,ind2,ind2]}]

(* trans[ject] introduces the functions (Subscript[e, cc])[var] and their derivatives that are used in the following for the Fourier transformation. *)

trans[ject_]:=ReplaceAll[ject,{(Subscript[ϕ, cc_])[var_]->ϕ[cc] (Subscript[e, cc])[var],(Subscript[Dϕ, cc_])[var_,ind_]->ϕ[cc] (Subscript[De, cc])[var,ind],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->ϕ[cc] (Subscript[DDe, cc])[var,ind1,ind2],A[var_,ind_]->(Subscript[e, 9])[var,ind] γ,DA[var_,ind1_,ind2_]->(Subscript[De, 9])[var,ind1,ind2] γ,DDA[var_,ind1_,ind2_,ind3_]->(Subscript[DDe, 9])[var,ind1,ind2,ind3] γ,(Subscript[W, dd_])[var_,ind_]->(Subscript[e, dd+9])[var,ind] W[dd],(Subscript[DW, dd_])[var_,ind1_,ind2_]->(Subscript[De, dd+9])[var,ind1,ind2] W[dd],(Subscript[DDW, dd_])[var_,ind1_,ind2_,ind3_]->(Subscript[DDe, dd+9])[var,ind1,ind2,ind3] W[dd]}]

(* numbers[ject] extracts the indices of the fields. *)

numbers[ject_]:=Table[Complement[ject,{γ,(Subscript[W, 1]),(Subscript[W, 2])}][[ii,2]],{ii,Length[Complement[ject,{γ,(Subscript[W, 1]),(Subscript[W, 2])}]]}]

(* setzero[ject] deletes the fields that do not occur as external fields.  *)

setzero[ject_]:=Do[tempject=ject;Do[tempject=ReplaceAll[tempject,{Subscript[ϕ,nonexnumbers[[bb]]][x_]->0,Subscript[Dϕ,nonexnumbers[[bb]]][a_,b_]->0,Subscript[DDϕ,nonexnumbers[[bb]]][a_,b_,c_]->0}],{bb,Length[nonexnumbers]}];If[MemberQ[nonexpar,γ],tempject=ReplaceAll[tempject,{A[__]->0,DA[__]->0,DDA[__]->0}]];Do[If[MemberQ[nonexpar,(Subscript[W, ii])],tempject=ReplaceAll[tempject,{(Subscript[W, ii])[__]->0,(Subscript[DW, ii])[__]->0,(Subscript[DDW, ii])[__]->0}]],{ii,8}];Return[tempject]]

(* eom[ject,ord] inserts the solution of the classical equation of motion. *)

eom[ject_,ord_]:=Do[Return[ex2ord[setzero[ReplaceAll[ject,{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]+(f^2) (Subscript[pg, cc])[var] (Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]+(f^2) (Subscript[Dpg, cc])[var,ind](Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]+(f^2) (Subscript[DDpg, cc])[var,ind1,ind2](Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext]}]],ord]]]

(* eomp2[ject,ord,ftr] is the same as eom[ject,ord] with an additional factor ftr and the inclusion of the physical decay constants. *)

eomp2[ject_,ord_,ftr_]:=Do[Return[ex2ord[setzero[ReplaceAll[ject,{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]fp[[cc]]+ftr*(f^2) (Subscript[pg, cc])[var](((fp[[cc]]))^2)(Subscript[Φ, cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],ext],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]fp[[cc]]+ftr*(f^2) (Subscript[Dpg, cc])[var,ind](((fp[[cc]]))^2)(Subscript[Φ, cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],ext],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]fp[[cc]]+ftr*(f^2) (Subscript[DDpg, cc])[var,ind1,ind2](((fp[[cc]]))^2)(Subscript[Φ, cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],ext]}]],ord]]]

(* At order (p^2) the decay constant in the chiral limit F has to be replaced by the physical decay constant ((Subscript[F, π]), (Subscript[F, K]) or (Subscript[F, η])) for each pseudoscalar meson by multiplying with the corresponding factor fp[[cc]]. *)

fp={F/(Subscript[F, π]),F/(Subscript[F, π]),F/(Subscript[F, π]),F/(Subscript[F, K]),F/(Subscript[F, K]),F/(Subscript[F, K]),F/(Subscript[F, K]),F/(Subscript[F, η])};

(* diff[ject] and its subfunctions deriv[ject,list] and listtrans[ject]perform the functional differentiation with respect to the external fields. *)

diff[ject_]:=deriv[trans[ject],listtrans[expar]]

deriv[ject_,list_]:=ReplaceAll[Fold[D,ject,list],{ϕ[cc_]->0,γ->0,W[cc_]->0}]

listtrans[ject_]:=ReplaceAll[ject,{(Subscript[ϕ, cc_])->ϕ[cc],(Subscript[W, cc_])->W[cc]}]

(* fourier[ject] and its subfunctions change, fourphi and fourmultiphi perform a Fourier transformation. *)

fourier[ject_]:=Do[object=Expand[ject];elist={(Subscript[e, aaa])[var_],(Subscript[De, aaa])[var_,ind_],(Subscript[DDe, aaa])[var_,ind1_,ind2_],(Subscript[e, aaa])[var_,ind1_],(Subscript[De, aaa])[var_,ind1_,ind2_],(Subscript[DDe, aaa])[var_,ind1_,ind2_,ind3_]};clist={1,2,3,4,5,6};out=0;Do[If[Length[1+object]>2,arb=object[[ii]],arb=object];Do[count=1;se=If[aaa<9,1,4];Do[tl=Position[arb,((elist[[ll]])^(expo_))];fourmultiphi[arb,tl,clist[[ll]]];tl=Position[arb,elist[[ll]]];fourphi[arb,tl,clist[[ll]]],{ll,se,se+2}],{aaa,1,12}];out+=arb,{ii,1,Length[1+object]-1}];Return[out]]

change[num_,1,vari_,inde_,inde2_,inde3_]:=(E^(I (Subscript[(p[num]), sub]) vari))

change[num_,2,vari_,inde_,inde2_,inde3_]:=-I (Subscript[(Subscript[(p[num]), sub]), inde])(E^(I (Subscript[(p[num]), sub]) vari))

change[num_,3,vari_,inde_,inde2_,inde3_]:=(-1) (Subscript[(Subscript[(p[num]), sub]), inde]) (Subscript[(Subscript[(p[num]), sub]), inde2])(E^(I (Subscript[(p[num]), sub]) vari))

change[num_,4,vari_,inde_,inde2_,inde3_]:=(E^(I (Subscript[(p[num]), sub]) vari))(Subscript[(ϵ[(Subscript[(p[num]), sub])]), inde3])

change[num_,5,vari_,inde_,inde2_,inde3_]:=-I (Subscript[(Subscript[(p[num]), sub]), inde])(E^(I (Subscript[(p[num]), sub]) vari))(Subscript[(ϵ[(Subscript[(p[num]), sub])]), inde3])

change[num_,6,vari_,inde_,inde2_,inde3_]:=(-1) (Subscript[(Subscript[(p[num]), sub]), inde]) (Subscript[(Subscript[(p[num]), sub]), inde2])(E^(I (Subscript[(p[num]), sub]) vari))(Subscript[(ϵ[(Subscript[(p[num]), sub])]), inde3])

fourphi[barb_,tl_,gate_]:=Do[Do[sub=barb[[Part[tl,ii]]][[0,2]];vari[ii]=barb[[Part[tl,ii]]][[1]];If[Length[barb[[Part[tl,ii]]]]>1,inde[ii]=barb[[Part[tl,ii]]][[2]]];If[Length[barb[[Part[tl,ii]]]]>2,inde2[ii]=barb[[Part[tl,ii]]][[3]]];If[gate>3,inde3[ii]=barb[[Part[tl,ii]]][[-1]]],{ii,1,Length[tl]}];If[tl≠{},new=Table[change[count+ll-1,gate,vari[ll],inde[ll],inde2[ll],inde3[ll]],{ll,1,Length[tl]}];rlist=Table[barb[[tl[[nn]]]]->new[[nn]],{nn,Length[tl]}];arb=ReplaceAll[barb,rlist];count+=Length[tl]]]

fourmultiphi[barb_,tl_,gate_]:=Do[Clear[temp];Do[sub=barb[[Part[tl,ii]]][[1,0,2]];vari[ii]=barb[[Part[tl,ii]]][[1,1]];If[Length[barb[[Part[tl,ii]]][[1]]]>1,inde[ii]=barb[[Part[tl,ii]]][[1,2]]];If[Length[barb[[Part[tl,ii]]][[1]]]>2,inde2[ii]=barb[[Part[tl,ii]]][[1,3]]];If[gate>3,inde3[ii]=barb[[Part[tl,ii]]][[1,-1]]];expo[ii]=barb[[Part[tl,ii]]][[2]];temp[ii]=1;Do[temp[ii]*=change[count+kk-1,gate,vari[ii],inde[ii],inde2[ii],inde3[ii]],{kk,1,expo[ii]}];count+=expo[ii],{ii,1,Length[tl]}];If[tl≠{},new=Table[temp[ll],{ll,Length[tl]}];rlist=Table[barb[[tl[[nn]]]]->new[[nn]],{nn,Length[tl]}];arb=ReplaceAll[barb,rlist]]]

(* extract[ject,func] and its subfunctions changealph, countalph and zero extract the momenta that have to be inserted in the constituent functions of the generating functional. *)

extract[ject_,func_]:=Do[exobject=Expand[ject];erge=0;Do[If[Length[1+exobject]>2,arb=exobject[[ii]],arb=exobject];zero[arb];tl=Position[arb,x];If[tl!={},Do[xx-=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]+1,-1]],{jj,1,Length[tl]}]];tl=Position[arb,y];If[tl!={},Do[yy+=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]+1,-1]],{jj,1,Length[tl]}]];arb=ReplaceAll[arb,{x->0,y->0,z->0}];tl=Position[arb,μ];If[tl!={},uu=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];tl=Position[arb,ν];If[tl!={},vv=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];tl=Position[arb,λ];If[tl!={},ww=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];arb=ReplaceAll[arb,{((Subscript[(Subscript[(p[var_]), ind1_]), ind2_])^2)->sp[(Subscript[(p[var]), ind1]),(Subscript[(p[var]), ind1])],(Subscript[ϵ, ind1_]^2)->sp[ϵ,ϵ],(Subscript[any_, α[7]])->epst[ξ,any],(Subscript[any_, α[8]])->epst[ρ,any],(Subscript[any_, α[9]])->epst[σ,any],(Subscript[any_, α[10]])->epst[τ,any]}];countalph[arb];Do[tl=Position[arb,Part[alphlist,kk]];If[tl!={},Do[alph[jj]=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]-1,-1]],{jj,1,Length[tl]}];aa[kk]=sp[alph[1],alph[2]]],{kk,1,numa}];arb2=ReplaceAll[arb,{(Subscript[(Subscript[(p[var_]), ind1_]), ind2_])->1,(Subscript[(ϵ[xx_]), ind_])->1,epst[ind_,any_]->(Subscript[any, ind])}];
		factor=arb2;Do[factor=factor*aa[kk],{kk,1,numa}];erge+=factor*func[xx,yy,uu,vv,ww,0,0,0],{ii,1,Length[1+exobject]-1}];Return[erge]]

changealph[ject_]:=ReplaceAll[ject,{(Subscript[α, x])->α[1],(Subscript[α, y])->α[2],(Subscript[α, z])->α[3],(Subscript[α, ext])->α[4],(Subscript[α, int])->α[5],(Subscript[β, x])->α[6],(Subscript[ξ, x])->α[7],(Subscript[ρ, x])->α[8],(Subscript[σ, x])->α[9],(Subscript[τ, x])->α[10],(Subscript[ζ, x])->α[11],(Subscript[ζ, y])->α[12],(Subscript[ζ, z])->α[13]}]

countalph[ject_]:=Do[numa=0;alphlist={};αα=Array[α,13];
Do[tl=Position[ject,αα[[kkk]]];
If[tl≠{},numa=numa+1;
alphlist=Append[alphlist,αα[[kkk]]]],{kkk,1,13}]]

zero[ject_]:=Do[eext=0;xx=0;yy=0; uu=0;vv=0;ww=0;Do[aa[kkk]=0,{kkk, 1, 20}]]

(* extract[ject,func,ms1,ms2,ms3] is the version of extract[]with different mass arguments ms1, ms2 and ms3 for the particles propagating in the loop. *)

extract[ject_,func_,ms1_,ms2_,ms3_]:=Do[exobject=Expand[ject];erge=0;Do[If[Length[1+exobject]>2,arb=exobject[[ii]],arb=exobject];zero[arb];countalph[arb];tl=Position[arb,x];If[tl≠{},Do[xx-=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]+1,-1]],{jj,1,Length[tl]}]];tl=Position[arb,y];If[tl≠{},Do[yy+=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]+1,-1]],{jj,1,Length[tl]}]];arb=ReplaceAll[arb,{x->0,y->0,z->0}];tl=Position[arb,μ];If[tl≠{},uu=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];tl=Position[arb,ν];If[tl≠{},vv=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];tl=Position[arb,λ];If[tl≠{},ww=Extract[arb,ReplacePart[tl[[1]],Last[tl[[1]]]-1,-1]]];arb=ReplaceAll[arb,{((Subscript[(Subscript[(p[var_]), ind1_]), ind2_])^2)->sp[(Subscript[(p[var]), ind1]),(Subscript[(p[var]), ind1])],(Subscript[ϵ, ind1_]^2)->sp[ϵ,ϵ]}];countalph[arb];Do[tl=Position[arb,Part[alphlist,kk]];If[tl≠{},Do[alph[jj]=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]-1,-1]],{jj,1,Length[tl]}];aa[kk]=sp[alph[1],alph[2]]],{kk,1,numa}];arb=ReplaceAll[arb,{(Subscript[(Subscript[(p[var_]), ind1_]), ind2_])->1,(Subscript[(ϵ[xx_]), ind_])->1}];
		factor=arb;Do[factor=factor*aa[kk],{kk,1,numa}];erge+=factor*func[xx,yy,uu,vv,ww,ms1,ms2,ms3],{ii,1,Length[1+exobject]-1}];Return[erge]]

(* momenta[ject] adds the extracted momenta and takes the minus signs out of the index brackets. *)

momenta[ject_]:=If[Length[1+ject]>2,Sum[ReplaceAll[(Subscript[(Part[Expand[ject],jj]), ind]),{(Subscript[(-moment_), ind])->-(Subscript[(moment), ind])}],{jj,1,Length[ject]}],ReplaceAll[(Subscript[(Expand[ject]), ind]),{(Subscript[(-moment_), ind])->-(Subscript[(moment), ind])}]]

(* momenta2[ject] acts like momenta[ject] in case of a product of extracted momenta. *)

momenta2[ject_]:=ReplaceAll[momenta[ject],{ind->ind1}]*ReplaceAll[momenta[ject],{ind->ind2}]

(* extracttree[ject] inserts the propagators that are outside of the loop. *)

extracttree[ject_]:=Do[erge=0;exobject=Expand[ject];Do[Do[If[Length[1+exobject]>2,arb=exobject[[ii]],arb=exobject];tl=Position[arb,ext];eext=0;If[tl≠{},Do[eext-=Extract[arb,ReplacePart[tl[[jj]],Last[tl[[jj]]]+1,-1]],{jj,1,Length[tl]}]]];arb=ReplaceAll[arb,{ext->0}];
erge+=Simplify[ReplaceAll[arb,{(Subscript[pg, par_])[var_]->-((E^(-I var eext))/(sp[eext,eext]-((Subscript[M, par])^2))),(Subscript[Dpg, par_])[var_,ind_]->-I*momenta[eext]*((E^(-I var eext))/(sp[eext,eext]-((Subscript[M, par])^2))),(Subscript[DDpg, par_])[var_,ind1_,ind2_]->momenta2[eext]*((E^(-I var eext))/(sp[eext,eext]-((Subscript[M, par])^2)))}]],{ii,1,Length[1+exobject]-1}];Return[erge]]

(* remnonex[ject,func] generates matrices where the non-external meson fields are removed. *)

remnonex[ject_,func_]:=Do[result=ject;Do[If[Not[MemberQ[expar,(Subscript[ϕ, ii])]],result=ReplaceAll[result,{(Subscript[ϕ, ii])[var]->0}]],{ii,1,8}];func=ReplaceAll[Function[var,rtemp],{rtemp->result}]]

(* remnonex[ject,func] is used to remove all except one non-external meson fields. *)

remnonex[ject_,func_,xxlist_]:=Do[result=ject;Do[If[Not[MemberQ[xxlist,(Subscript[ϕ, ii])]],result=ReplaceAll[result,{(Subscript[ϕ, ii])[var]->0}],If[MemberQ[nonexpar,(Subscript[ϕ, ii])],result=ReplaceAll[result,{(Subscript[ϕ, ii])[var]->(Subscript[ϕ, ii])[var]*Length[nonexmeson]}]]],{ii,1,8}];func=ReplaceAll[Function[var,rtemp],{rtemp->result}]]

(* parmomenta[ject] generates a list, called list1,  of the external momenta in the standard form (Subscript[p[j], i]) and another list, called list2, of the external momenta in the user-defined form. *)

parmomenta[ject_]:=Do[list1={};list2={};Do[mult[jj]=1;equalparmom[jj]={},{jj,SortsOfParticles}];object=ReplaceAll[ject,{γ->(Subscript[γ, 9]),(Subscript[W, dd_])->(Subscript[W, dd+9])}];Do[sort=Part[object,ii][[0,2]];list1=Append[list1,Subscript[p[mult[sort]],Part[object,ii][[0,2]]]];list2=Append[list2,Part[object,ii][[1]]];mult[sort]+=1;equalparmom[sort]=Append[equalparmom[sort],list2[[ii]]],{ii,Length[object]}]]

(* changemomenta[ject] replaces the standard momenta by the user-defined momenta. *)

changemomenta[ject_]:=ReplaceAll[ject,Table[list1[[ii]]->list2[[ii]],{ii,Length[list1]}]]

(* sym[ject] takes indistinguishable particles into account. *)

sym[ject_]:=Do[temp=ject;Do[If[Length[equalparmom[kk]]>1,permu=Permutations[equalparmom[kk]];temp=(1/(Length[permu[[1]]]!))(temp+Sum[ReplaceAll[temp,Table[permu[[1,ii]]->permu[[jj,ii]],{ii,1,Length[permu[[1]]]}]],{jj,2,Length[permu]}])],{kk,SortsOfParticles}];Return[temp]]

(* makemasslist[ject] generates a list of the particle masses. *)

makemasslist[ject_]:=Do[masslist={};expargam={};expargam=ReplaceAll[ject,{γ->(Subscript[γ, 9]),(Subscript[W, dd_])->(Subscript[W, dd+9])}];Do[masslist=Append[masslist,expargam[[ii,2]]],{ii,Length[expargam]}]]

(* isolim[ject] assigns the isospin limit values to the masses. *)

isolim[ject_]:=ReplaceAll[ject,{(Subscript[M, 1])->(Subscript[M, π]),(Subscript[M, 2])->(Subscript[M, π]),(Subscript[M, 3])->(Subscript[M, π]),(Subscript[M, 4])->(Subscript[M, K]),(Subscript[M, 5])->(Subscript[M, K]),(Subscript[M, 6])->(Subscript[M, K]),(Subscript[M, 7])->(Subscript[M, K]),(Subscript[M, 8])->(Subscript[M, η])}]

(* mattable[smat,amat] calculates the matrices that occur in the one-loop functional. *)

mattable[samat_]:=Do[mat2[samat[[ll]]]=matrix2[samat[[ll]],stp];mat4[samat[[ll]]]=If[particleorder>4,matrix4[samat[[ll]],stp],mat2[samat[[ll]]]],{ll,1,Length[samat]}]

(* getamp[ject,func,ΦΦ] extracts the amplitude. *)

getamp[ject_,func_,ΦΦ_]:=Do[out1=ReplaceAll[diff[ject],{Φ->ΦΦ}];out2=changealph[fourier[out1]];out3=extracttree[out2];parmomenta[expar2];out4=simp[sym[changemomenta[ReplaceAll[extract[out3,func],{delta[__]->1,f->1}]]]];Return[out4]]

(* An additional argument φφ that is used for including bilinear terms can be added. *)

getamp[ject_,func_,ΦΦ_,φφ_]:=Do[out1=ReplaceAll[diff[ject],{Φ->ΦΦ,φ->φφ}];
out2=changealph[fourier[out1]];out3=extracttree[out2];
parmomenta[expar2];
out4=simp[sym[changemomenta[ReplaceAll[extract[out3,func],{delta[__]->1,f->1}]]]];
Return[out4]]

(* getamp[ject,func,ΦΦ,ms1,ms2,ms3] and getamp[ject,func,ΦΦ,φφ,ms1,ms2,ms3] extract the loop amplitude. *)

getamp[ject_,func_,ΦΦ_,ms1_,ms2_,ms3_]:=Do[out1=ReplaceAll[diff[ject],{Φ->ΦΦ}];out2=changealph[fourier[out1]];out3=extracttree[out2];parmomenta[expar2];out4=simp[sym[changemomenta[ReplaceAll[extract[out3,func,ms1,ms2,ms3],{delta[__]->1,f->1}]]]];Return[out4]]

getamp[ject_,func_,ΦΦ_,φφ_,ms1_,ms2_,ms3_]:=Do[out1=ReplaceAll[diff[ject],{Φ->ΦΦ,φ->φφ}];
out2=changealph[fourier[out1]];out3=extracttree[out2];
parmomenta[expar2];
out4=simp[sym[changemomenta[ReplaceAll[extract[out3,func,ms1,ms2,ms3],{delta[__]->1,f->1}]]]];Return[out4]]

(* second[ject] and third[ject] change the spatial arguments and indices of ject. *)

second[ject_]:=ReplaceAll[ject,{x->y,μ->ν}]

third[ject_,ll_,terml_]:=ReplaceAll[ject,If[MemberQ[smatrixlist,terml[[ll,1]]] && MemberQ[smatrixlist,terml[[ll,2]]] && MemberQ[amatrixlist,terml[[ll,3]]],{x->z},{x->z,μ->λ}]]

(* special[ject] changes μ into the summation index (Subscript[α, int]). *)

special[ject_]:=ReplaceAll[ject,{μ->(Subscript[α, int])}]

(* *)

(* transform[ject] and its subfunctions csp, dsp, scaleout, subscaleout, and msp transform the scalar products into their standard form. *)

transform[ject_]:=Do[slist2=Sort[list2];epslist2=slist2;Do[epslist2=Append[epslist2,ϵ[slist2[[ii]]]],{ii,Length[slist2]}];Return[isolim[ReplaceAll[ReplaceAll[scalout[ReplaceAll[ReplaceAll[ReplaceAll[ject,{sp->dsp}],{sp[-x_,-y_]->sp[x,y],sp[-x_,y_]->-sp[x,y],sp[x_,-y_]->-sp[x,y]}],{sp[x_,0]->0,sp[0,x_]->0,sp[0,0]->0}]],{sp->csp}],{sp->msp}]]]]

csp[x_,y_]:=If[Position[epslist2,sp[x,y][[1]]][[1,1]]>Position[epslist2,sp[x,y][[2]]][[1,1]],sp[y,x],sp[x,y]]

dsp[x_,y_]:=Distribute[sp[x,y]]

subscalout[ject_]:=ReplaceAll[ject,{sp[x_*y_,z_]->x*sp[y,z],sp[x_,y_*z_]->y*sp[x,z]}]

scalout[ject_]:=subscalout[subscalout[ject]]

msp[x_,y_]:=If[x==y,massind=masslist[[Position[list2,x][[1,1]]]];If[massind<9,Return[((Subscript[M, massind])^2)],sp[x,x]],,sp[x,y]]

(* sppoint[ject] changes sp[x,y] into x.y . *)

sppoint[ject_]:=ReplaceAll[ject,{sp->Dot}]

(* replacements[ject] replaces the standard scalar products by the user-defined expressions. *)

replacements[ject_]:=ReplaceAll[ReplaceAll[addlep[ject],Table[ReplaceAll[replist1[[ii]],{sp->csp}]->replist2[[ii]],{ii,Length[replist1]}]],{sp[x_,x_]->(x^2),((x_)^μ)->(Subscript[x, μ])}]

(* final[ject] renormalizes the coupling constants and transforms the third argument of the C function. *)

final[ject_]:=ReplaceAll[ReplaceAll[ReplaceAll[ject,{(Subscript[L, ii_])->((Subscript[L, ii])^r),(Subscript[N, ii_])->((Subscript[N, ii])^r),(Subscript[R, ii_])->((Subscript[R, ii])^r),(Subscript[(OverHat[N]), ii_])->((Subscript[(OverHat[N]), ii])^r),(Subscript[(OverHat[R]), ii_])->((Subscript[(OverHat[R]), ii])^r)}],{C[x_,y_,z_,ma_,mb_,mc_]->C[x,y,Simplify[(1/2)(x+y-z)],ma,mb,mc]}],{(0^r)->0,((OverHat[0])^r)->0}]

(* epscorr[ject] removes a possible minus sign in the argument of ϵ. *)

epscorr[ject_]:=ReplaceAll[ject,{ϵ[-q_]->ϵ[q]}]

(* particleordercount[ject] determines the order of the external particles. *)

particleordercount[ject_]:=Do[countinglist=ReplaceAll[ject,{ϕ->tei,γ->(Subscript[tei, 9]),(Subscript[W, 1])->(Subscript[tei, 10]),(Subscript[W, 2])->(Subscript[tei, 11])}];valuelist={1,1,1,1,1,1,1,1,2,1,1};temporlist=Table[countinglist[[jj,2]],{jj,Length[countinglist]}];particleorder=Sum[valuelist[[temporlist[[ii]]]],{ii,Length[temporlist]}]]

(* *)

(* insertbilin[ject, ord] is used to insert the bilinear terms that come from parts which are bilinear in the meson fields. *)

insertbilin[ject_,ord_]:=Sum[ex2ordbi[ReplaceAll[Expand[ject][[jj]],{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]+wp (Subscript[φ, cc])[var],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]+wp dete[(Subscript[φ, cc])[var],var,ind],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]+wp  dete[dete[(Subscript[φ, cc])[var],var,ind1],var,ind2]}],ord],{jj,Length[Expand[ject]]}]

(* fieldtransform[ject, name1, name2] and its subfunction bilin[ject] generate the bilinar terms. *)

fieldtransform[ject_,name1_,name2_]:=Do[ap={2,1,3,5,4,7,6,8};bilinmat=bilin[ject];Do[erg=Sum[If[MS[[ap[[ii]],ap[[ii]]]]==MS[[jj,jj]],0,,1/(MS[[ap[[ii]],ap[[ii]]]]-MS[[jj,jj]])ReplaceAll[bilinmat[[ii,jj]],{(p^2)->MS[[jj,jj]]}]]*(Subscript[ϕ, jj])[var],{jj,8}];(Subscript[name1, ap[[ii]]])=ReplaceAll[Function[{var},eer],{eer->erg}];(Subscript[name2, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->funcchange[trans[erg]]}];(Subscript[(OverBar[name2]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind},eer],{eer->funcchange[trans[de[erg,var,ind]]]}];(Subscript[(Overscript[name2, (OverBar[_])]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind1,ind2},eer],{eer->funcchange[trans[de[de[erg,var,ind1],var,ind2]]]}],{ii,8}]]

fieldtransform2[ject_,name1_,name2_]:=Do[ap={2,1,3,5,4,7,6,8};bilinmat=bilin[ject];Do[erg=Sum[If[MS[[ap[[ii]],ap[[ii]]]]==MS[[jj,jj]],0,,1/(MSP[[ap[[ii]],ap[[ii]]]]-MSP[[jj,jj]])ReplaceAll[bilinmat[[ii,jj]],{(p^2)->MSP[[jj,jj]]}]]*(Subscript[ϕ, jj])[var],{jj,8}];(Subscript[name1, ap[[ii]]])=ReplaceAll[Function[{var},eer],{eer->erg}];(Subscript[name2, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->funcchange[trans[erg]]}];(Subscript[(OverBar[name2]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind},eer],{eer->funcchange[trans[de[erg,var,ind]]]}];(Subscript[(Overscript[name2, (OverBar[_])]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind1,ind2},eer],{eer->funcchange[trans[de[de[erg,var,ind1],var,ind2]]]}],{ii,8}]]

MSP={{(Subscript[Mp, π]^2),0,0,0,0,0,0,0},{0,(Subscript[Mp, π]^2),0,0,0,0,0,0},{0,0,(Subscript[Mp, π]^2),0,0,0,0,0},{0,0,0,(Subscript[Mp, K]^2),0,0,0,0},{0,0,0,0,(Subscript[Mp, K]^2),0,0,0},{0,0,0,0,0,(Subscript[Mp, K]^2),0,0},{0,0,0,0,0,0,(Subscript[Mp, K]^2),0},{0,0,0,0,0,0,0,(Subscript[Mp, η]^2)}};

bilin[ject_]:=Do[expartemp=expar;expar2temp=expar2;ret=Table[expar={(Subscript[ϕ, ii]),(Subscript[ϕ, jj])};expar2={(Subscript[ϕ, ii])[p],(Subscript[ϕ, jj])[-p]};ReplaceAll[minpick[Simplify[getamp[ject,delta,Φ]]],{sp[p,p]->(p^2)}],{ii,8},{jj,8}];expar2=expar2temp;expar=expartemp;Return[ret]]

minpick[ject_]:=ReplaceAll[ject,{sp[-x_,-y_]->sp[x,y],sp[-x_,y_]->-sp[x,y],sp[x_,-y_]->-sp[x,y]}]

(* getordbi[ject, ord] supplies ject to order ord in f and to order one in the weak perturbation wp. *)

getordbi[ject_,ord_]:=Simplify[ReplaceAll[D[ject,{f,ord},wp],{f->0,wp->0}]]*((f^ord)/(ord !))

(* ex2ordbi[ject, ord] supplies ject up to order ord in f and to order one in wp. *)

ex2ordbi[ject_,ord_]:=Sum[getordbi[ject,ii],{ii,0,ord}]

(* addbilin[ject, ord] adds the bilinear terms. *)

addbilin[ject_,ord_]:=Sum[ex2ordbi[setzero[ReplaceAll[If[Length[1+ject]<2,Expand[ject][[jj]],ject],{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]+wp  (Subscript[φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]+wp  (Subscript[(OverBar[φ]), cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var,ind],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]+wp  (Subscript[(Overscript[φ, (OverBar[_])]), cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var,ind1,ind2]}]],ord],{jj,If[Length[1+ject]<2,Length[Expand[ject]],1]}]

(* addbilin2[ject, ord] includes also the decay constant factors fp[[cc]]. *)

addbilinp2[ject_,ord_]:=Sum[ex2ordbi[setzero[ReplaceAll[If[Length[1+ject]<2,Expand[ject][[jj]],ject],{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]fp[[cc]]+wp*(((fp[[cc]]))^2)* (Subscript[φ, cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],var],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]fp[[cc]]+wp *(((fp[[cc]]))^2)* (Subscript[(OverBar[φ]), cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],var,ind],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]fp[[cc]]+wp *(((fp[[cc]]))^2)* (Subscript[(Overscript[φ, (OverBar[_])]), cc])[ϕ[1]fp[[1]],ϕ[2]fp[[2]],ϕ[3]fp[[3]],ϕ[4]fp[[4]],ϕ[5]fp[[5]],ϕ[6]fp[[6]],ϕ[7]fp[[7]],ϕ[8]fp[[8]],γ,W[1],W[2],var,ind1,ind2]}]],ord],{jj,If[Length[1+ject]<2,Length[Expand[ject]],1]}]

(* eombilin[ject, ord] combines eom[ject, ord] and addbilin[ject, ord]. *)

eombilin[ject_,ord_]:=Sum[ex2ordbi[setzero[ReplaceAll[If[Length[1+ject]<2,Expand[ject][[jj]],ject],{(Subscript[ϕ, cc_])[var_]->(Subscript[ϕ, cc])[var]+wp  (Subscript[φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var]+(f^2) (Subscript[pg, cc])[var] (Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext],(Subscript[Dϕ, cc_])[var_,ind_]->(Subscript[Dϕ, cc])[var,ind]+wp  (Subscript[(OverBar[φ]), cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var,ind]+(f^2) (Subscript[Dpg, cc])[var,ind](Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext],(Subscript[DDϕ, cc_])[var_,ind1_,ind2_]->(Subscript[DDϕ, cc])[var,ind1,ind2]+wp  (Subscript[(Overscript[φ, (OverBar[_])]), cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],var,ind1,ind2]+(f^2) (Subscript[DDpg, cc])[var,ind1,ind2](Subscript[Φ, cc])[ϕ[1],ϕ[2],ϕ[3],ϕ[4],ϕ[5],ϕ[6],ϕ[7],ϕ[8],γ,W[1],W[2],ext]}]],ord],{jj,If[Length[1+ject]<2,Length[Expand[ject]],1]}]

(* fieldtransformmono[name] generates the terms that are proportional to one pseudoscalar meson field *)

fieldtransformmono[name_]:=Do[ap={2,1,3,5,4,7,6,8};monolinvec={0,0,0,0,0,monophi6,monophi7,0};Do[erg=monolinvec[[ii]];(Subscript[name, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->funcchange[trans[erg]]}];(Subscript[(OverBar[name]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind},eer],{eer->funcchange[trans[de[erg,var,ind]]]}];(Subscript[(Overscript[name, (OverBar[_])]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind1,ind2},eer],{eer->funcchange[trans[de[de[erg,var,ind1],var,ind2]]]}],{ii,8}]]

(* momzero[ject] removes the terms containing derivates that act on pseudoscalar meson fields. *)

momzero[ject_]:=Simplify[ReplaceAll[ject,{(Subscript[Dϕ, i_])[__]->0,(Subscript[DDϕ, i_])[__]->0}]]

(* tadep[ject, ord] provides the tadpole terms. *)

tadep[ject_,ord_]:=(1/2)Sum[A[MS[[kk,kk]]]ject[x,0,kk,kk,stp,ord,u,U],{kk,8}]

(* imlist is a list of the momenta of the external particles. *)

imlist=Table[ExternalParticles[[ii,1]],{ii,Length[ExternalParticles]}];

(* momentumcons[ject] replaces the last momentum of imlist by the sum of the other external momenta times (-1). *)

momentumcons[ject_]:=Simplify[ReplaceAll[ject,{(Subscript[(imlist[[-1]]), i_])->Sum[-(Subscript[(imlist[[jj]]), i]),{jj,Length[imlist]-1}]}]]

(* eps[am,bm,cm,dm] and epsorder[ject] are used to simplify the terms containing the epsilon-tensor. *)

eps[am_,bm_,cm_,dm_]:=(Subscript[((seps=Sort[{am,bm,cm,dm}])[[1]]), ξ]) (Subscript[(seps[[2]]), ρ]) (Subscript[(seps[[3]]), σ])(Subscript[(seps[[4]]), τ])*Signature[{am,bm,cm,dm}]

epsorder[ject_]:=Simplify[ReplaceAll[ReplaceAll[Expand[momentumcons[Coefficient[ject,ϵ^(ξρστ)]]],{(Subscript[(am_), ξ]) (Subscript[(bm_), τ]) (Subscript[(cm_), σ]) (Subscript[(dm_), ρ])->epstemp[am,bm,cm,dm]}],{epstemp->eps}]]*ϵ^(ξρστ)+ReplaceAll[ject,{ϵ^(ξρστ)->0}]

(* The function addlep[ject] adds the leptonic matrix element to the amplitude. *)

addlep[ject_]:=If[wpar==1,ReplaceAll[ReplaceAll[Expand[ject],{sp[conmom_,ϵ[wmo]]->lμ (conmom^μ),(ϵ^(ξρστ)) (Subscript[(mtm1_), ξ]) (Subscript[(mtm2_), σ]) (Subscript[(mtm3_), ρ]) (Subscript[(ϵ[wmo]), τ])->(ϵ^(ξρσμ)) (Subscript[(mtm1), ξ]) (Subscript[(mtm2), σ]) (Subscript[(mtm3), ρ]) lμ,(ϵ^(ξρστ)) (Subscript[(mtm1_), ξ]) (Subscript[(mtm2_), σ]) (Subscript[(ϵ[wmo]), ρ]) (Subscript[(mtm3_), τ])->(ϵ^(ξμστ)) (Subscript[(mtm1), ξ]) (Subscript[(mtm2), σ]) (Subscript[(mtm3), τ]) lμ,(ϵ^(ξρστ)) (Subscript[(mtm1_), ξ]) (Subscript[(ϵ[wmo]), σ]) (Subscript[(mtm2_), ρ]) (Subscript[(mtm3_), τ])->(ϵ^(ξρμτ)) (Subscript[(mtm1), ξ]) (Subscript[(mtm2), ρ]) (Subscript[(mtm3), τ]) lμ,(ϵ^(ξρστ)) (Subscript[(ϵ[wmo]), ξ]) (Subscript[(mtm1_), σ]) (Subscript[(mtm2_), ρ]) (Subscript[(mtm3_), τ])->(ϵ^(μρστ)) (Subscript[(mtm1), σ]) (Subscript[(mtm2), ρ]) (Subscript[(mtm3), τ]) lμ}],{(wmo^μ)->-Sum[((imlistwow[[ii]])^μ),{ii,Length[imlistwow]}]}],ject]

(* The function normalization[ject] assigns the correct normalization to the amplitude. *)

normalization[ject_]:=pullout[If[wpar==1,ReplaceAll[ex2ordns[ex2ordns[ex2ordns[ex2ordns[ex2ordns[ex2ordns[getordns[ReplaceAll[ject,{Dot->sp,(Subscript[C, w])->(Subscript[G, F])}],1,lμ],5,e],1,(Subscript[G, F])],1,(Subscript[V, ud])],1,(Subscript[V, us])],1,(Subscript[(OverHat[V]), ud])],1,(Subscript[(OverHat[V]), us])],{lμ->(Subscript[lep, μ]),sp->Dot,(Subscript[x_, μ])->(x^μ)}],sppoint[ject]]]

pullout[ject_]:=ReplaceAll[ex2ordns[ex2ordns[subpullout[ReplaceAll[ject,{Dot->sp,Complex[ui1_,ui2_]->ui1+imag*ui2}]],5,e],1,imag],{sp->Dot,imag->I}]

subpullout[ject_]:=If[stcheck==1,getordns[ject,1,(Subscript[G, 8])]+getordns[ject,1,(Subscript[G, 27])]+getordns[ject,1,(Subscript[(OverHat[G]), 8])]+getordns[ject,1,(Subscript[(OverHat[G]), 27])],ject]

(* *)

particleordercount[expar];

If[(particleorder>4 &&(MemberQ[expar,(Subscript[ϕ, 4])]||MemberQ[expar,(Subscript[ϕ, 5])]||MemberQ[expar,(Subscript[ϕ, 6])]||MemberQ[expar,(Subscript[ϕ, 7])]||MemberQ[expar,(Subscript[ϕ, 8])])&&Count[expar,γ]<2)||phi6==1,norecurs=1]

If[MemberQ[expar,(Subscript[W, 1])],lep=(OverHat[l]),lep=l];

(* nonexpar is the list of the fields that do not occur as external particles. *)

nonexpar=Complement[AllParticles,expar];

(* nonexmeson is the list of the pseudoscalar mesons that do not occur as external particles. *)

nonexmeson=Complement[AllParticles,Union[expar,{γ,(Subscript[W, 1]),(Subscript[W, 2])}]];

(* nonexnumbers is the list of the numbers of the fields that do not occur as external particles. *)

nonexnumbers=numbers[nonexpar];

makemasslist[expar]

(* In the following the matrices U and u that appear in the chiral Lagrangian are defined up to different orders: *)

(Subscript[U, 0])[x_]:= (Subscript[δ, 3])

(Subscript[U, 1])[x_]:= (Subscript[δ, 3])+(I/F)Ψ[x]*f

(Subscript[U, 2])[x_]:= (Subscript[U, 1])[x]+(1/(2!)) ((((I/F)))^2)Ψ2[x]*(f^2)

(Subscript[U, 3])[x_]:= (Subscript[U, 2])[x]+(1/(3!)) ((((I/F)))^3)Ψ3[x]*(f^3)

(Subscript[U, 4])[x_]:=(Subscript[U, 3])[x]+(1/(4!)) ((((I/F)))^4)Ψ3[x].Ψ[x]*(f^4)

(Subscript[U, 5])[x_]:=(Subscript[U, 4])[x]+(1/(5!)) ((((I/F)))^5)Ψ3[x].Ψ2[x]*(f^5)

(Subscript[U, 6])[x_]:=(Subscript[U, 5])[x]+(1/(6!)) ((((I/F)))^6)Ψ3[x].Ψ3[x]*(f^6)

(Subscript[u, 0])[x_]:= (Subscript[δ, 3])

(Subscript[u, 1])[x_]:= (Subscript[δ, 3])+(I/(2F))Ψ[x]*f

(Subscript[u, 2])[x_]:= (Subscript[u, 1])[x]+(1/((2^2)*2!)) ((((I/F)))^2)Ψ2[x]*(f^2)

(Subscript[u, 3])[x_]:= (Subscript[u, 2])[x]+(1/((2^3)*3!)) ((((I/F)))^3)Ψ3[x]*(f^3)

(Subscript[u, 4])[x_]:= (Subscript[u, 3])[x]+(1/((2^4)*4!)) ((((I/F)))^4)Ψ3[x].Ψ[x]*(f^4)

(Subscript[u, 5])[x_]:=(Subscript[u, 4])[x]+(1/((2^5)*5!)) ((((I/F)))^5)Ψ3[x].Ψ2[x]*(f^5)

(Subscript[u, 6])[x_]:=(Subscript[u, 5])[x]+(1/((2^6)*6!)) ((((I/F)))^6)Ψ3[x].Ψ3[x]*(f^6)

(Subscript[U, -1])[x_]:= (Subscript[δ, 3])

(Subscript[u, -1])[x_]:= (Subscript[δ, 3])

(Subscript[U, -2])[x_]:= (Subscript[δ, 3])

(Subscript[u, -2])[x_]:= (Subscript[δ, 3])

(Subscript[U, -3])[x_]:= (Subscript[δ, 3])

(Subscript[u, -3])[x_]:= (Subscript[δ, 3])

(Subscript[(OverHat[name_]), jj_])[x_]:=conj[(Subscript[name, jj])[x]]

(* In the following a quantity with a hat is complex conjugated. *)

(Subscript[δ, 3])={{1,0,0},{0,1,0},{0,0, 1}};

(* Ψ[x] is the matrix containing the pseudoscalar meson fields, Ψ2[x] is the squared matrix and Ψ3[x] is Ψ[x].Ψ[x].Ψ[x] . *)

Ψ[x_]:={{(Subscript[ϕ, 3])[x]+(((Subscript[ϕ, 8])[x])/(Sqrt[3])),(Sqrt[2]) (Subscript[ϕ, 1])[x],(Sqrt[2]) (Subscript[ϕ, 4])[x]},{(Sqrt[2]) (Subscript[ϕ, 2])[x],-(Subscript[ϕ, 3])[x]+(((Subscript[ϕ, 8])[x])/(Sqrt[3])),(Sqrt[2]) (Subscript[ϕ, 6])[x]},{(Sqrt[2]) (Subscript[ϕ, 5])[x],(Sqrt[2]) (Subscript[ϕ, 7])[x],-((2 (Subscript[ϕ, 8])[x])/(Sqrt[3]))}}

Ψ2[x_]:={{2 ((Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+(Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+((((((Subscript[ϕ, 3])[x])/(Sqrt[2]))+(((Subscript[ϕ, 8])[x])/(Sqrt[6]))))^2)),2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 7])[x]+2 (Sqrt[(2/3)]) (Subscript[ϕ, 1])[x] (Subscript[ϕ, 8])[x],2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 6])[x]+(1/3) (Sqrt[2]) (Subscript[ϕ, 4])[x] (3 (Subscript[ϕ, 3])[x]-(Sqrt[3]) (Subscript[ϕ, 8])[x])},{2 (Subscript[ϕ, 5])[x] (Subscript[ϕ, 6])[x]+2 (Sqrt[(2/3)]) (Subscript[ϕ, 2])[x] (Subscript[ϕ, 8])[x],2 ((Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+(Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+((((((Subscript[ϕ, 3])[x])/(Sqrt[2]))-(((Subscript[ϕ, 8])[x])/(Sqrt[6]))))^2)),2 (Subscript[ϕ, 2])[x] (Subscript[ϕ, 4])[x]-(1/3) (Sqrt[2]) (Subscript[ϕ, 6])[x] (3 (Subscript[ϕ, 3])[x]+(Sqrt[3]) (Subscript[ϕ, 8])[x])},{2 (Subscript[ϕ, 2])[x] (Subscript[ϕ, 7])[x]+(1/3) (Sqrt[2]) (Subscript[ϕ, 5])[x] (3 (Subscript[ϕ, 3])[x]-(Sqrt[3]) (Subscript[ϕ, 8])[x]),2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 5])[x]-(1/3) (Sqrt[2]) (Subscript[ϕ, 7])[x] (3 (Subscript[ϕ, 3])[x]+(Sqrt[3]) (Subscript[ϕ, 8])[x]),2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(4/3) (((Subscript[ϕ, 8])[x])^2)}}

Ψ3[x_]:={{2 (Sqrt[2]) (Subscript[ϕ, 2])[x] (Subscript[ϕ, 4])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^3)+(Sqrt[3]) (((Subscript[ϕ, 3])[x])^2) (Subscript[ϕ, 8])[x]+((((Subscript[ϕ, 8])[x])^3)/(3 (Sqrt[3])))+(Subscript[ϕ, 3])[x] (4 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+(((Subscript[ϕ, 8])[x])^2))+2 (Subscript[ϕ, 1])[x] ((Sqrt[2]) (Subscript[ϕ, 5])[x] (Subscript[ϕ, 6])[x]+(Subscript[ϕ, 2])[x] ((Subscript[ϕ, 3])[x]+(Sqrt[3]) (Subscript[ϕ, 8])[x])),(Sqrt[2]) (Subscript[ϕ, 1])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2)),(Sqrt[2]) (Subscript[ϕ, 4])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2))},{(Sqrt[2]) (Subscript[ϕ, 2])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2)),2 (Sqrt[2]) (Subscript[ϕ, 2])[x] (Subscript[ϕ, 4])[x] (Subscript[ϕ, 7])[x]-(((Subscript[ϕ, 3])[x])^3)+(Sqrt[3]) (((Subscript[ϕ, 3])[x])^2) (Subscript[ϕ, 8])[x]+((((Subscript[ϕ, 8])[x])^3)/(3 (Sqrt[3])))-(Subscript[ϕ, 3])[x] (4 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 8])[x])^2))+(Subscript[ϕ, 1])[x] (2 (Sqrt[2]) (Subscript[ϕ, 5])[x] (Subscript[ϕ, 6])[x]-2 (Subscript[ϕ, 2])[x] ((Subscript[ϕ, 3])[x]-(Sqrt[3]) (Subscript[ϕ, 8])[x])),(Sqrt[2]) (Subscript[ϕ, 6])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2))},{(Sqrt[2]) (Subscript[ϕ, 5])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2)),(Sqrt[2]) (Subscript[ϕ, 7])[x] (2 (Subscript[ϕ, 1])[x] (Subscript[ϕ, 2])[x]+2 (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]+2 (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]+(((Subscript[ϕ, 3])[x])^2)+(((Subscript[ϕ, 8])[x])^2)),2 (Sqrt[2]) (Subscript[ϕ, 1])[x] (Subscript[ϕ, 5])[x] (Subscript[ϕ, 6])[x]+2 (Sqrt[2]) (Subscript[ϕ, 2])[x] (Subscript[ϕ, 4])[x] (Subscript[ϕ, 7])[x]+2 ((Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x]-(Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x]) (Subscript[ϕ, 3])[x]-2 (Sqrt[3]) (Subscript[ϕ, 4])[x] (Subscript[ϕ, 5])[x] (Subscript[ϕ, 8])[x]-2 (Sqrt[3]) (Subscript[ϕ, 6])[x] (Subscript[ϕ, 7])[x] (Subscript[ϕ, 8])[x]-((8 (((Subscript[ϕ, 8])[x])^3))/(3 (Sqrt[3])))}}

(* ph[var,ind] is the matrix that contains the photon field. *)

ph[var_,ind_]:=-e A[var,ind]*{{(2/3),0,0},{0,-(1/3),0},{0,0,-(1/3)}}

(* wbos[var,ind] is the matrix that contains the W boson fields. The constant (Subscript[C, w]) is defined later. *)

wbos[var_,ind_]:=-(Sqrt[2])(Subscript[C, w])((Subscript[W, 1])[var,ind]*{{0,(Subscript[V, ud]),(Subscript[V, us])},{0,0,0},{0,0,0}}+(Subscript[W, 2])[var,ind]*{{0,0,0},{(Subscript[(OverHat[V]), ud]),0,0},{(Subscript[(OverHat[V]), us]),0,0}})

(* If no external W bosons occur, wbos[var,ind] is set equal to zero. *)

If[!MemberQ[expar,(Subscript[W, cc_])],wbos[var_,ind_]:=0;wpar=0,wpar=1];

(* wmo is the momentum of the W boson. *)

If[wpar==1,wmo=expar2[[Position[expar,(Subscript[W, cc_])][[1,1]],1]]];

(* imlistwow is the list of the momenta without the momentum of the W boson. *)

imlistwow=Complement[imlist,{wmo}];

(* set0[ject] and dot0[ject] allow to remove products with zero. *)

set0[ject_]:=ReplaceAll[ject,{Dot[__,0,__]->0,Dot[__,0]->0,Dot[0,__]->0}]

dot0[ject_]:=set0[set0[ject]]

(* r[var,ind] and l[var,ind] are external fields. *)

r[var_,ind_]:=ph[var,ind]*(f^2)

l[var_,ind_]:=ph[var,ind]*(f^2)+wbos[var,ind]*f

(* FL[var,ind1,ind2] and FR[var,ind1,ind2] are tensors that appear in the chiral Lagrangian. *)

FL[var_,ind1_,ind2_]:=de[l[var,(Subscript[ind2, var])],var,(Subscript[ind1, var])]-de[l[var,(Subscript[ind1, var])],var,(Subscript[ind2, var])]-I kom[l[var,(Subscript[ind1, var])],l[var,(Subscript[ind2, var])]]

FR[var_,ind1_,ind2_]:=de[r[var,(Subscript[ind2, var])],var,(Subscript[ind1, var])]-de[r[var,(Subscript[ind1, var])],var,(Subscript[ind2, var])]-I kom[r[var,(Subscript[ind1, var])],r[var,(Subscript[ind2, var])]]

(* MM is proportional to the quark mass matrix. *)

MM:={{((Subscript[M, π])^2),0,0},{0,((Subscript[M, π])^2),0},{0,0,2((Subscript[M, K])^2)-((Subscript[M, π])^2)}}

χ[var_]:=MM

(OverHat[χ])[var_]:=MM

(* MS is the squared mass matrix of the pseudoscalar mesons. *)

MS=DiagonalMatrix[{((Subscript[M, π])^2),((Subscript[M, π])^2),((Subscript[M, π])^2),((Subscript[M, K])^2),((Subscript[M, K])^2),((Subscript[M, K])^2),((Subscript[M, K])^2),((Subscript[M, η])^2)}];

(* kom[m1, m2] is the commutator. *)

kom[m1_,m2_]:=m1.m2-m2.m1

(* akom[m1, m2] is the anticommutator. *)

akom[m1_,m2_]:=m1.m2+m2.m1

(* With remnonex the matrices (Subscript[Ur, ii]) and (Subscript[ur, ii]) are generated that contain only external fields. *)

Remove[Ur]

Do[remnonex[(Subscript[U, ii])[var],(Subscript[Ur, ii])],{ii,-3,6}]

Remove[ur]

Do[remnonex[(Subscript[u, ii])[var],(Subscript[ur, ii])],{ii,-3,6}]

(* (Subscript[Uw[jj], ii]) and (Subscript[ur[jj], ii]) contain the external fields and one additional pseudoscalar meson field. *)

Remove[Uw]

Do[Do[remnonex[(Subscript[U, ii])[var],(Subscript[Uw[jj], ii]),Union[expar,{nonexmeson[[jj]]}]],{ii,-3,6}],{jj,1,Length[nonexmeson]}]

Remove[uw]

Do[Do[remnonex[(Subscript[u, ii])[var],(Subscript[uw[jj], ii]),Union[expar,{nonexmeson[[jj]]}]],{ii,-3,6}],{jj,1,Length[nonexmeson]}]

(* *)

(* kov[ject,var,ind] is the covariant derivative. *)

kov[ject_,var_,ind_]:=de[ject,var,ind]-I r[var,ind].ject+I ject.l[var,ind]

(* kovc[ject,var,ind] is the covariant derivative with an interchange of l[var,ind] and r[var,ind]. *)

kovc[ject_,var_,ind_]:=de[ject,var,ind]-I l[var,ind].ject+I ject.r[var,ind]

(* ublist is the list of the complex conjugated matrices that contain the pseudoscalar meson fields. *)

ublist=Union[{(OverHat[U]),(OverHat[Ur])},Table[(OverHat[Uw[ii]]),{ii,8}]];

(* In the following the Lagrangians are defined: *)

d[ind_,UU_,ord_]:=If[MemberQ[ublist,UU],kovc[(Subscript[UU, ord])[var],var,(Subscript[ind, var])],kov[(Subscript[UU, ord])[var],var,(Subscript[ind, var])]]

d[ind_,UU_,ord_,var_]:=If[MemberQ[ublist,UU],kovc[(Subscript[UU, ord])[var],var,(Subscript[ind, var])],kov[(Subscript[UU, ord])[var],var,ind]]

mt[U_,ord_]:=(Subscript[(OverHat[U]), ord])[var].χ[var]+(OverHat[χ])[var].(Subscript[U, ord])[var]

mtmin[U_,ord_]:=(Subscript[(OverHat[U]), ord])[var].χ[var]-(OverHat[χ])[var].(Subscript[U, ord])[var]

L[U_,ind_,ord_]:=I(Subscript[(OverHat[U]), ord])[var].d[ind,U,ord]

(* lagp2[var, ord, u, U] is the strong Lagrangian of order (p^2) *)

lagp2[pointvar_,ord_,u_,U_]:=ReplaceAll[((F^2)/4)(Tr2ord[d[α,U,ord-1].d[α,(OverHat[U]),ord-1]+(Subscript[(OverHat[U]), ord])[var].χ[var]+(OverHat[χ])[var].(Subscript[U, ord])[var],ord]-ReplaceAll[Tr2ord[d[α,U,1].d[α,(OverHat[U]),1]+(Subscript[(OverHat[U]), 2])[var].χ[var]+(OverHat[χ])[var].(Subscript[U, 2])[var],2],{(Subscript[W, cc_])[__]->0}]),{var->pointvar}]

(* lagwp2[pointvar, ord, fromord, u, U] is the nonleptonic weak Lagrangian of order (p^2) *)

lagwp2[pointvar_,ord_,fromord_,u_,U_]:=Do[lgtemp=(F^4)ReplaceAll[If[stp==1,ex2ord[(Subscript[G, 8])(L[U,α,ord-1].L[U,α,ord-1])[[2,3]],ord]+ex2ord[(Subscript[G, 27])(L[U,α,ord-1][[2,3]] L[U,α,ord-1][[1,1]]+(2/3)L[U,α,ord-1][[2,1]] L[U,α,ord-1][[1,3]]),ord],ex2ord[(Subscript[(OverHat[G]), 8])(L[U,α,ord-1].L[U,α,ord-1])[[3,2]],ord]+ex2ord[(Subscript[(OverHat[G]), 27])(L[U,α,ord-1][[3,2]] L[U,α,ord-1][[1,1]]+(2/3)L[U,α,ord-1][[1,2]] L[U,α,ord-1][[3,1]]),ord]],{var->pointvar}];lgge=lgtemp-ex2ord[lgtemp,fromord-1];Return[lgge]]

(* lagp4[var, ord, u, U] is the strong Lagrangian of order (p^4). In lagp4[pointvar,ord,u,U] the terms producing two-point functions are subtracted. Their contributions to the amplitude are already included in the definition of the physical decay constants (Subscript[F, π]),  (Subscript[F, K]) and (Subscript[F, η]). *)

lagp4[pointvar_,ord_,u_,U_]:=ReplaceAll[(Subscript[L, 1]) ex2ord[ Tr2ord[d[α,(OverHat[U]),ord-3].d[α,U,ord-3],ord-2]*Tr2ord[d[β,(OverHat[U]),ord-3].d[β,U,ord-3],ord-2],ord]+(Subscript[L, 2])ex2ord[Tr2ord[d[α,(OverHat[U]),ord-3].d[β,U,ord-3],ord-2] *Tr2ord[d[α,(OverHat[U]),ord-3].d[β,U,ord-3],ord-2],ord]+(Subscript[L, 3]) Tr2ord[d[α,(OverHat[U]),ord-3].d[α,U,ord-3].d[β,(OverHat[U]),ord-3].d[β,U,ord-3],ord]+(Subscript[L, 4]) Tr2ord[d[α,(OverHat[U]),ord-1].d[α,U,ord-1],ord]* Tr2ord[mt[U,ord-2],ord]+(Subscript[L, 5]) Tr2ord[d[α,(OverHat[U]),ord-1].d[α,U,ord-1].mt[U,ord-2],ord]+(Subscript[L, 6]) ex2ord[((Tr2ord[mt[U,ord],ord])^2),ord]+(Subscript[L, 7]) ex2ord[((Tr2ord[mtmin[U,ord],ord])^2),ord]+(Subscript[L, 8]) Tr2ord[(OverHat[χ])[var].(Subscript[U, ord])[var].(OverHat[χ])[var].(Subscript[U, ord])[var]+χ[var].(Subscript[(OverHat[U]), ord])[var].χ[var].(Subscript[(OverHat[U]), ord])[var],ord]-I (Subscript[L, 9]) Tr2ord[FR[var,α,β].d[α,U,ord-2].d[β,(OverHat[U]),ord-2]+FL[var,α,β].d[α,(OverHat[U]),ord-2].d[β,U,ord-2],ord]+(Subscript[L, 10]) Tr2ord[(Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var].FL[var,α,β],ord]-ReplaceAll[ex2ord[(Subscript[L, 4]) Tr2ord[d[α,(OverHat[U]),1].d[α,U,1],2]* Tr2ord[mt[U,1],2]+(Subscript[L, 5]) Tr2ord[d[α,(OverHat[U]),1].d[α,U,1].mt[U,1],2]+(Subscript[L, 6]) ((Tr2ord[mt[U,2],2])^2)+(Subscript[L, 7]) ex2ord[((Tr2ord[mtmin[U,2],2])^2),2]+(Subscript[L, 8]) Tr2ord[(OverHat[χ])[var].(Subscript[U, 2])[var].(OverHat[χ])[var].(Subscript[U, 2])[var]+χ[var].(Subscript[(OverHat[U]), 2])[var].χ[var].(Subscript[(OverHat[U]), 2])[var],2],2],{(Subscript[W, cc_])[__]->0}],{var->pointvar}]

(* wzw[pointvar,ord,u,U] contains the terms coming from the anomaly functional. Only terms containing external gauge fields are included *)

wzw[pointvar_,ord_,u_,U_]:=-ϵ^(ξρστ)(I/(16 (π^2)))ReplaceAll[Tr2ord[(Subscript[U, ord-4])[var].l[var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), ord-4])[var].r[var,(Subscript[τ, var])],ord]+(1/4)Tr2ord[(Subscript[U, ord-4])[var].l[var,(Subscript[ξ, var])].(Subscript[(OverHat[U]), ord-4])[var].r[var,(Subscript[ρ, var])].(Subscript[U, ord-4])[var].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), ord-4])[var].r[var,(Subscript[τ, var])],ord]+I Tr2ord[(Subscript[U, ord-3])[var].de[l[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), ord-3])[var].r[var,(Subscript[τ, var])],ord]+I Tr2ord[de[r[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].(Subscript[U, ord-3])[var].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), ord-3])[var].r[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), ord-4])[var].r[var,(Subscript[σ, var])].(Subscript[U, ord-4])[var].l[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[(OverHat[U]), ord-3])[var].de[(Subscript[U, ord-2])[var],var,(Subscript[ξ, var])].(Subscript[(OverHat[U]), ord-3])[var].de[r[var,(Subscript[σ, var])],var,(Subscript[ρ, var])].(Subscript[U, ord-3])[var].l[var,(Subscript[τ, var])],ord]-Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ξ, var])].(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), ord-4])[var].r[var,(Subscript[σ, var])].(Subscript[U, ord-4])[var].l[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[(OverHat[U]), ord-3])[var].de[(Subscript[U, ord-2])[var],var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].de[l[var,(Subscript[τ, var])],var,(Subscript[σ, var])],ord]+Tr2ord[(Subscript[(OverHat[U]), ord-3])[var].de[(Subscript[U, ord-2])[var],var,(Subscript[ξ, var])].de[l[var,(Subscript[σ, var])],var,(Subscript[ρ, var])].l[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].l[var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]+(1/2)Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ξ, var])].(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), ord-4])[var].de[(Subscript[U, ord-3])[var],var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]-Tr2ord[l[var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].l[var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]-(1/4)Tr2ord[l[var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].l[var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]-I Tr2ord[de[l[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].l[var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]-I Tr2ord[de[r[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].l[var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]-(-I Tr2ord[de[l[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].r[var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]-I Tr2ord[de[r[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].r[var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]-(1/4) Tr2ord[r[var,(Subscript[ξ, var])].l[var,(Subscript[ρ, var])].r[var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]-Tr2ord[r[var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].r[var,(Subscript[σ, var])].l[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[U, -3+ord])[var].de[(Subscript[(OverHat[U]), -2+ord])[var],var,(Subscript[ξ, var])].de[r[var,(Subscript[σ, var])],var,(Subscript[ρ, var])].r[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[U, -3+ord])[var].de[(Subscript[(OverHat[U]), -2+ord])[var],var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].de[r[var,(Subscript[τ, var])],var,(Subscript[σ, var])],ord]+I Tr2ord[de[l[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].(Subscript[(OverHat[U]), -3+ord])[var].r[var,(Subscript[σ, var])].(Subscript[U, -3+ord])[var].l[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].r[var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]+I Tr2ord[(Subscript[(OverHat[U]), -3+ord])[var].de[r[var,(Subscript[ρ, var])],var,(Subscript[ξ, var])].r[var,(Subscript[σ, var])].(Subscript[U, -3+ord])[var].l[var,(Subscript[τ, var])],ord]+(1/2) Tr2ord[(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[U, -3+ord])[var].de[(Subscript[(OverHat[U]), -2+ord])[var],var,(Subscript[ξ, var])].(Subscript[U, -3+ord])[var].de[l[var,(Subscript[σ, var])],var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), -3+ord])[var].r[var,(Subscript[τ, var])],ord]+Tr2ord[(Subscript[(OverHat[U]), -4+ord])[var].r[var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].r[var,(Subscript[σ, var])].(Subscript[U, -4+ord])[var].l[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ξ, var])].r[var,(Subscript[ρ, var])].(Subscript[U, -4+ord])[var].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), -4+ord])[var].r[var,(Subscript[τ, var])],ord]-I Tr2ord[(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ξ, var])].(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ρ, var])].(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[σ, var])].r[var,(Subscript[τ, var])],ord]-Tr2ord[(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ξ, var])].(Subscript[U, -4+ord])[var].de[(Subscript[(OverHat[U]), -3+ord])[var],var,(Subscript[ρ, var])].(Subscript[U, -4+ord])[var].l[var,(Subscript[σ, var])].(Subscript[(OverHat[U]), -4+ord])[var].r[var,(Subscript[τ, var])],ord]+(1/4) Tr2ord[(Subscript[(OverHat[U]), -4+ord])[var].r[var,(Subscript[ξ, var])].(Subscript[U, -4+ord])[var].l[var,(Subscript[ρ, var])].(Subscript[(OverHat[U]), -4+ord])[var].r[var,(Subscript[σ, var])].(Subscript[U, -4+ord])[var].l[var,(Subscript[τ, var])],ord]),{var->pointvar}]

(* Q[2,3] and Q[3,2] are used to project out the 23-and 32-components of a matrix. *)

Q[2,3]={{0,0,0},{0,0,0},{0,1,0}};

Q[3,2]={{0,0,0},{0,0,1},{0,0,0}};

(* If ΔS=1, Q[2,3] is used. In case of ΔS=-1, Q[3,2] is used. *)

proj=If[stp==1,Q[2,3],Q[3,2]];

(* The sign of ΔS determines whether the coupling constants are complex conjugated. *)

GG=If[stp==1,G,(OverHat[G])];

NN=If[stp==1,N,(OverHat[N])];

RR=If[stp==1,R,(OverHat[R])];

(* lagwp48[pointvar, ord, fromord, u, U] is the octet part of the nonleptonic weak Lagrangian of order (p^4) *)

lagwp48[pointvar_,ord_,fromord_,u_,U_]:=Do[lgtemp=(Subscript[GG, 8])*(F^2)ReplaceAll[(Subscript[NN, 1])Tr2ord[proj.d[α,(OverHat[U]),ord-3].d[α,U,ord-3].d[β,(OverHat[U]),ord-3].d[β,U,ord-3],ord]+(Subscript[NN, 2])Tr2ord[proj.d[α,(OverHat[U]),ord-3].d[β,U,ord-3].d[β,(OverHat[U]),ord-3].d[α,U,ord-3],ord]+(Subscript[NN, 3])ex2ord[Tr2ord[proj.d[α,(OverHat[U]),ord-3].d[β,U,ord-3],ord-2]Tr2ord[d[α,(OverHat[U]),ord-3].d[β,U,ord-3],ord-2],ord]+(Subscript[NN, 4])ex2ord[Tr2ord[proj.d[α,(OverHat[U]),ord-3].(Subscript[U, ord-4])[var],ord-3]Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].d[α,U,ord-3].d[β,(OverHat[U]),ord-3].d[β,U,ord-3],ord-1],ord]+(Subscript[NN, 5])Tr2ord[proj.akom[mt[U,ord-2],d[α,(OverHat[U]),ord-1].d[α,U,ord-1]],ord]+(Subscript[NN, 6])ex2ord[Tr2ord[proj.d[α,(OverHat[U]),ord-1].(Subscript[U, ord-2])[var],ord-1]Tr2ord[(Subscript[(OverHat[U]), ord-2])[var].d[α,U,ord-1].mt[U,ord-2],ord-1],ord]+(Subscript[NN, 7])ex2ord[Tr2ord[proj.mt[U,ord-2],ord-2]Tr2ord[d[α,(OverHat[U]),ord-1].d[α,U,ord-1],ord],ord]+(Subscript[NN, 8])ex2ord[Tr2ord[proj.d[α,(OverHat[U]),ord-1].d[α,U,ord-1],ord]Tr2ord[mt[U,ord-2],ord-2],ord]+(Subscript[NN, 9])Tr2ord[proj.kom[mtmin[U,ord-2],d[α,(OverHat[U]),ord-1].d[α,U,ord-1]],ord]+(Subscript[NN, 10])Tr2ord[proj.mt[U,ord].mt[U,ord],ord]+(Subscript[NN, 11])ex2ord[Tr2ord[proj.mt[U,ord],ord]Tr2ord[mt[U,ord],ord],ord]+(Subscript[NN, 12])Tr2ord[proj.mtmin[U,ord].mtmin[U,ord],ord]+(Subscript[NN, 13])ex2ord[Tr2ord[proj.mtmin[U,ord],ord]Tr2ord[mtmin[U,ord],ord],ord]+(Subscript[NN, 14]) I Tr2ord[proj.akom[FL[var,α,β]+(Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var],d[α,(OverHat[U]),ord-2].d[β,U,ord-2]],ord]+(Subscript[NN, 15]) I Tr2ord[proj.d[α,(OverHat[U]),ord-2].((Subscript[U, ord-3])[var].FL[var,α,β].(Subscript[(OverHat[U]), ord-3])[var]+FR[var,α,β]).d[β,U,ord-2],ord]+(Subscript[NN, 16]) I Tr2ord[proj.akom[FL[var,α,β]-(Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var],d[α,(OverHat[U]),ord-2].d[β,U,ord-2]],ord]+(Subscript[NN, 17]) I Tr2ord[proj.d[α,(OverHat[U]),ord-2].((Subscript[U, ord-3])[var].FL[var,α,β].(Subscript[(OverHat[U]), ord-3])[var]-FR[var,α,β]).d[β,U,ord-2],ord]+2(Subscript[NN, 18])Tr2ord[proj.(FL[var,α,β].(Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]+(Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var].FL[var,α,β]),ord]+(Subscript[NN, 28]) I ϵ^(ξρστ)ex2ord[Tr2ord[proj.d[ξ,(OverHat[U]),ord-3].(Subscript[U, ord-4])[var],ord-3]Tr2ord[(Subscript[(OverHat[U]), ord-4])[var].d[ρ,U,ord-3].d[σ,(OverHat[U]),ord-3].d[τ,U,ord-3],ord-1],ord]+2(Subscript[NN, 29])ϵ^(ξρστ)Tr2ord[proj.kom[(Subscript[(OverHat[U]), ord-3])[var].FR[var,σ,τ].(Subscript[U, ord-3])[var],d[ξ,(OverHat[U]),ord-2].d[ρ,U,ord-2]],ord]+(Subscript[NN, 30])ϵ^(ξρστ)ex2ord[Tr2ord[proj.(Subscript[(OverHat[U]), ord-3])[var].d[ξ,U,ord-2],ord-2]Tr2ord[(FL[var,σ,τ]+(Subscript[(OverHat[U]), ord-3])[var].FR[var,σ,τ].(Subscript[U, ord-3])[var]).d[ρ,(OverHat[U]),ord-2].(Subscript[U, ord-3])[var],ord-1],ord]+(Subscript[NN, 31])ϵ^(ξρστ)ex2ord[Tr2ord[proj.(Subscript[(OverHat[U]), ord-3])[var].d[ξ,U,ord-2],ord-2]Tr2ord[(FL[var,σ,τ]-(Subscript[(OverHat[U]), ord-3])[var].FR[var,σ,τ].(Subscript[U, ord-3])[var]).d[ρ,(OverHat[U]),ord-2].(Subscript[U, ord-3])[var],ord-1],ord],{var->pointvar}];lgge=lgtemp-ex2ord[lgtemp,fromord-1];Return[lgge]]

sublagwp427[pointvar_,ord_,fromord_,u_,U_,i_,j_,k_,l_]:=Do[lgtemp=ReplaceAll[ex2ord[(Subscript[RR, 1])(L[U,α,ord-3].L[U,α,ord-3])[[i,j]] (L[U,β,ord-3].L[U,β,ord-3])[[k,l]],ord] +ex2ord[(Subscript[RR, 2])(L[U,α,ord-3].L[U,β,ord-3])[[i,j]] (L[U,α,ord-3].L[U,β,ord-3])[[k,l]],ord] +ex2ord[(Subscript[RR, 3])(L[U,α,ord-3].L[U,β,ord-3])[[i,j]] (L[U,β,ord-3].L[U,α,ord-3])[[k,l]],ord] +ex2ord[(Subscript[RR, 4])(L[U,β,ord-3].L[U,α,ord-3].L[U,β,ord-3])[[k,l]] L[U,α,ord-3][[i,j]] ,ord]+ex2ord[(Subscript[RR, 5])L[U,α,ord-3][[i,j]] (akom[L[U,α,ord-3],L[U,β,ord-3].L[U,β,ord-3]])[[k,l]],ord] +ex2ord[(Subscript[RR, 6])Tr[L[U,α,ord-3].L[U,α,ord-3]]*L[U,β,ord-3][[i,j]] L[U,β,ord-3][[k,l]],ord] +ex2ord[(Subscript[RR, 7])(L[U,α,ord-1].L[U,α,ord-1])[[k,l]] mt[U,ord][[i,j]] ,ord]+ex2ord[(Subscript[RR, 8])L[U,α,ord-1][[k,l]] (akom[mt[U,ord],L[U,α,ord-1]])[[i,j]] ,ord]-ex2ord[I (Subscript[RR, 9])L[U,α,ord-1][[k,l]] I (kom[mtmin[U,ord],L[U,α,ord-1]])[[i,j]] ,ord]+ex2ord[(Subscript[RR, 10])Tr[mt[U,ord]]*L[U,α,ord-1][[i,j]] L[U,α,ord-1][[k,l]] ,ord]+ex2ord[(Subscript[RR, 11])mt[U,ord][[i,j]] mt[U,ord][[k,l]] ,ord]+ex2ord[(Subscript[RR, 12])(I mtmin[U,ord])[[i,j]] (I mtmin[U,ord])[[k,l]] ,ord]+ex2ord[I (Subscript[RR, 13])((Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var]+FL[var,α,β])[[i,j]] (L[U,α,ord-2].L[U,β,ord-2])[[k,l]] ,ord]+ex2ord[I (Subscript[RR, 14])(L[U,α,ord-2].L[U,β,ord-2])[[k,l]] ((Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var]-FL[var,α,β])[[i,j]] ,ord]+ex2ord[I (Subscript[RR, 15])(kom[(Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var]+FL[var,α,β],L[U,α,ord-2]])[[i,j]] L[U,β,ord-2][[k,l]],ord] +ex2ord[I (Subscript[RR, 16])(kom[((Subscript[(OverHat[U]), ord-3])[var].FR[var,α,β].(Subscript[U, ord-3])[var]-FL[var,α,β]),L[U,α,ord-2]])[[i,j]] L[U,β,ord-2][[k,l]] ,ord]+ex2ord[(Subscript[RR, 17])((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]+FL[var,α,β])[[i,j]] ((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]+FL[var,α,β])[[k,l]],ord] +ex2ord[(Subscript[RR, 18])((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]-FL[var,α,β])[[k,l]] ((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]+FL[var,α,β])[[i,j]] ,ord]+ex2ord[(Subscript[RR, 19])((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]-FL[var,α,β])[[i,j]] ((Subscript[(OverHat[U]), ord-2])[var].FR[var,α,β].(Subscript[U, ord-2])[var]-FL[var,α,β])[[k,l]],ord] +ex2ord[I (Subscript[RR, 20])ϵ^(ξρστ) (L[U,ξ,ord-3].L[U,ρ,ord-3])[[i,j]] (L[U,σ,ord-3].L[U,τ,ord-3])[[k,l]],ord] +ex2ord[I (Subscript[RR, 21])ϵ^(ξρστ) L[U,ξ,ord-3][[i,j]] (L[U,ρ,ord-3].L[U,σ,ord-3].L[U,τ,ord-3])[[k,l]] ,ord]+ex2ord[(Subscript[RR, 22])ϵ^(ξρστ) L[U,ξ,ord-2][[i,j]] (L[U,ρ,ord-2].((Subscript[(OverHat[U]), ord-3])[var].FR[var,σ,τ].(Subscript[U, ord-3])[var]+FL[var,σ,τ]))[[k,l]],ord] +ex2ord[(Subscript[RR, 23])ϵ^(ξρστ) L[U,ξ,ord-2][[i,j]] (L[U,ρ,ord-2].((Subscript[(OverHat[U]), ord-3])[var].FR[var,σ,τ].(Subscript[U, ord-3])[var]-FL[var,σ,τ]))[[k,l]] ,ord],{var->pointvar}];lgge=lgtemp-ex2ord[lgtemp,fromord-1];Return[lgge]]

(* lagwp427[pointvar, ord, fromord, u, U] is the 27-plet part of the nonleptonic weak Lagrangian of order (p^4). *)

lagwp427[pointvar_,ord_,fromord_,u_,U_]:=If[stp==1,(Subscript[G, 27])*(F^2)((1/3)(sublagwp427[pointvar,ord,fromord,u,U,2,1,1,3]+sublagwp427[pointvar,ord,fromord,u,U,1,3,2,1]+sublagwp427[pointvar,ord,fromord,u,U,2,3,1,1]+sublagwp427[pointvar,ord,fromord,u,U,1,1,2,3])-(1/6)(sublagwp427[pointvar,ord,fromord,u,U,2,2,2,3]+sublagwp427[pointvar,ord,fromord,u,U,2,3,2,2]+sublagwp427[pointvar,ord,fromord,u,U,2,3,3,3]+sublagwp427[pointvar,ord,fromord,u,U,3,3,2,3])),(Subscript[(OverHat[G]), 27])*(F^2)((1/3)(sublagwp427[pointvar,ord,fromord,u,U,1,2,3,1]+sublagwp427[pointvar,ord,fromord,u,U,3,1,1,2]+sublagwp427[pointvar,ord,fromord,u,U,3,2,1,1]+sublagwp427[pointvar,ord,fromord,u,U,1,1,3,2])-(1/6)(sublagwp427[pointvar,ord,fromord,u,U,2,2,3,2]+sublagwp427[pointvar,ord,fromord,u,U,3,2,2,2]+sublagwp427[pointvar,ord,fromord,u,U,3,2,3,3]+sublagwp427[pointvar,ord,fromord,u,U,3,3,3,2]))]

lagwp4[var_,ord_,fromord_,u_,U_]:=lagwp48[var,ord,fromord,u,U]+lagwp427[var,ord,fromord,u,U]

(* In the following the lists of the coupling constants are generated. *)

Llist=Table[(Subscript[L, ii]),{ii,10}];

Nlist=Union[Table[If[stp==1,(Subscript[N, ii]),(Subscript[(OverHat[N]), ii])],{ii,18}],If[stp==1,{(Subscript[N, 28]),(Subscript[N, 29]),(Subscript[N, 30]),(Subscript[N, 31])},{(Subscript[(OverHat[N]), 28]),(Subscript[(OverHat[N]), 29]),(Subscript[(OverHat[N]), 30]),(Subscript[(OverHat[N]), 31])}]];

Rlist=Union[Table[If[stp==1,(Subscript[R, ii]),(Subscript[(OverHat[R]), ii])],{ii,23}]];

constlist=If[stcheck==0,Llist,Union[Llist,Nlist,Rlist]];

(* Gell-Mann matrices: *)

(Subscript[λ, 1])={{0,1,0},{1,0,0},{0,0,0}}; (Subscript[λ, 2])={{0,-I,0},{I,0,0},{0,0,0}};  (Subscript[λ, 3])={{1,0,0},{0,-1,0},{0,0,0}};  (Subscript[λ, 4])={{0,0,1},{0,0,0},{1,0,0}};  (Subscript[λ, 5])={{0,0,-I},{0,0,0},{I,0,0}};  (Subscript[λ, 6])={{0,0,0},{0,0,1},{0,1,0}};  (Subscript[λ, 7])={{0,0,0},{0,0,-I},{0,I,0}};(Subscript[λ, 8])=(1/(Sqrt[3])){{1,0,0},{0,1,0},{0,0,-2}};

(* In the following expressions the solution of the classical equation of motion for each field is calculated. *)

Do[ap={2,1,3,5,4,7,6,8};Do[erg=funcchange[trans[Simplify[setzero[dephi[(1/(f^4))*ReplaceAll[lagp2[var,4,u,U],{(Subscript[W, cc_])[__]->0,(Subscript[DW, cc_])[__]->0}],ii]]]]];(Subscript[Φs, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->erg}],{ii,1,8}]]

If[wpar==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=ReplaceAll[Simplify[setzero[dephi[(1/(f^4))*lagp2[var,2,u,U],ii]]],{α->ζ}]*(((-1))/(sp[wmo,wmo]-((Subscript[M, ap[[ii]]])^2)));(Subscript[ϕsp2, ap[[ii]]])=ReplaceAll[Function[{var},eer],{eer->erg}];(Subscript[φsp2, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->funcchange[trans[erg]]}];(Subscript[(OverBar[φsp2]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind},eer],{eer->funcchange[trans[de[erg,var,ind]]]}]],{ii,1,8}]]

If[wpar==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=ReplaceAll[Simplify[setzero[dephi[(1/(f^4))*getordns[lagp4[var,2,u,U],1,(Subscript[C, w])],ii]]],{α->ζ}]*(((-1))/(sp[wmo,wmo]-((Subscript[M, ap[[ii]]])^2)));(Subscript[φsp4, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->funcchange[trans[erg]]}];(Subscript[(OverBar[φsp4]), ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var,ind},eer],{eer->funcchange[trans[de[erg,var,ind]]]}]],{ii,1,8}]]

If[wpar==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=funcchange[trans[wkorr[Simplify[setzero[dephi[(1/(f^4))*(lagp2[var,4,u,U]-lagp2[var,2,u,U]+If[wpar==1,getordns[ReplaceAll[ReplaceAll[insertbilin[lagp2[var,4,u,U],4],{φ->ϕsp2}],dete->de],1,(Subscript[C, w])],0]),ii]]]]]];(Subscript[Φsa, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->erg}],{ii,1,8}]]]

If[stcheck==1,fieldtransform[lagwp2[x,2,2,u,U],ϕwp2,φp2]]

If[stcheck==1,fieldtransform2[lagwp2[x,2,2,u,U],ϕwp2p,φp2p]]

If[stcheck==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=funcchange[trans[Simplify[setzero[dephi[(1/(f^4))*(lagwp2[var,4,3,u,U]+ReplaceAll[ReplaceAll[insertbilin[lagp2[var,4,u,U],4],{φ->ϕwp2p}],dete->de]),ii]]]]]+ReplaceAll[(Subscript[φp2p, ap[[ii]]])[[2]],{ϕ1->1,ϕ2->1,ϕ3->1,ϕ4->1,ϕ5->1,ϕ6->1,ϕ7->1,ϕ8->1,γ->1,W1->1,W2->1,(Subscript[e, i_])[var_]->(Subscript[Φs, i])[ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var]}];(Subscript[Φwp, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->erg}];(Subscript[Φw, ap[[ii]]])=ReplaceAll[(Subscript[Φwp, ap[[ii]]]),{Mp->M}],{ii,1,8}]]]

(* In some cases the trees coming from the solution of the classical equation of motion are set to zero and the following identities are needed. *)

(Subscript[0, cc_])[__]:=0

(Subscript[(OverBar[0]), cc_])[__]:=0

(Subscript[(Overscript[0, (OverBar[_])]), cc_])[__]:=0

(* *)

(* strongp2 is the order (p^2) part coming from the strong Lagrangian, where trees or bilinear terms have to be added. *)

(* version of April 2012 *)(*strongp2=dot0[dot0[eomp2[lagp2[x,4,u,U]-lagp2[x,2,u,U],6,If[stcheck==0,(1/2),1]]+Sum[addbilinp2[lagp2[x,6,uw[jj],Uw[jj]],6],{jj,Length[nonexmeson]}]/Length[nonexmeson]]];*)
(*new version of August 2012 *)
strongp2=dot0[dot0[eomp2[lagp2[x,4,u,U]-lagp2[x,2,u,U],6,If[stcheck==0,(1/2),1]]+Sum[addbilinp2[lagp2[x,6,uw[jj],Uw[jj]],6],{jj,Length[nonexmeson]}]/Length[nonexmeson]+eomp2[Simplify[lagp2[x,2,u,U]],2,1]]];

(* lagtreep2 is the complete order (p^2) part. *)

lagtreep2=strongp2+ReplaceAll[If[stcheck==0,getord[setzero[lagp2[x,6,ur,Ur]],6]+getord[setzero[lagp2[x,5,ur,Ur]],5],setzero[lagwp2[x,6,3,ur,Ur]]],Flatten[Table[{(Subscript[ϕ, ii])[x_]->(Subscript[ϕ, ii])[x] fp[[ii]],(Subscript[Dϕ, ii])[x_,ind_]->(Subscript[Dϕ, ii])[x,ind] fp[[ii]],(Subscript[DDϕ, ii])[x_,ind1_,ind2_]->(Subscript[DDϕ, ii])[x,ind1,ind2] fp[[ii]]},{ii,8}]]];

treep2temp2=getamp[lagtreep2,delta,If[stcheck==0,If[wpar==0,Φs,Φsa],Φwp],If[stcheck==0,If[wpar==0,0,φsp2],φp2p]];

treep2temp=ReplaceAll[treep2temp2,{Mp->M}];

(* treep2 is the tree amplitude of O((p^2)). *)

treep2=If[onlyloops==0,replacements[transform[ReplaceAll[treep2temp,{(F^(ii_))->((Subscript[F, π])^ii),((Subscript[F, K])^(ii_))->((Subscript[F, π])^ii),((Subscript[F, η])^(ii_))->((Subscript[F, π])^ii)}]]],0];

treeresult2temp=If[phi6==0,Simplify[treep2],ex2ordns2[Collect[FactorTerms[treep2],Power[__,-1]],4,(Subscript[F, π])]];

treeresult2=If[onlyloops==0,normalization[treeresult2temp],"Not evaluated"];

If[onlytreep2==1,result=treeresult2;title="Tree level amplitude of O(p^2):";treeresult4="Not evaluated";treeresult2and4="Not evaluated";loopresult="Not evaluated";output[result]];

If[onlytreep2==1,Quit[]]

etakorr[ject_]:=If[nogmo==1,ject,ReplaceAll[ject,{(Subscript[M, η])->(Sqrt[((4 (Subscript[M, K]^2)-(Subscript[M, π]^2))/3)])}],ReplaceAll[ject,{(Subscript[M, η])->(Sqrt[((4 (Subscript[M, K]^2)-(Subscript[M, π]^2))/3)])}]]

(* renomp4 is the order (p^4) amplitude due to wave function and mass renormalization. *)

renomp4=replacements[transform[ReplaceAll[getord[ReplaceAll[ReplaceAll[etakorr[treep2temp2],{(Subscript[M, π]^(ii_))->((Subscript[M, π])^ii)(1+(mo^2)(ii/2)(((A[((Subscript[M, π])^2)])/(2(F^2)))-(1/3)*((A[((Subscript[M, η])^2)])/(2(F^2)))-((8((Subscript[M, π])^2))/(F^2))(2(Subscript[L, 8])+2(Subscript[L, 6])-(Subscript[L, 4])-(Subscript[L, 5]))-((16((Subscript[M, K])^2))/(F^2))(2(Subscript[L, 6])-(Subscript[L, 4])))),(Subscript[M, K]^(ii_))->((Subscript[M, K])^ii)(1+(mo^2)(ii/2)((2/3)*((A[((Subscript[M, η])^2)])/(2(F^2)))-((16((Subscript[M, K])^2))/(F^2))(2(Subscript[L, 6])+(Subscript[L, 8])-(Subscript[L, 4])-(1/2)(Subscript[L, 5]))-((8((Subscript[M, π])^2))/(F^2))(2(Subscript[L, 6])-(Subscript[L, 4]))))}],{(F^(ii_))->(F^ii)(1+(mo^2)ii(-((A[((Subscript[M, π])^2)])/(F^2))-((A[((Subscript[M, K])^2)])/(2(F^2)))-((4((Subscript[M, π])^2))/(F^2))((Subscript[L, 5])+(Subscript[L, 4]))-((8((Subscript[M, K])^2))/(F^2))(Subscript[L, 4]))),((Subscript[F, K])^(ii_))->(F^ii)(1+(1/(8 (F^2)))((mo^2) ii(2 A[(Subscript[M, K]^2)]-5 A[(Subscript[M, π]^2)]+3 A[(Subscript[M, η]^2)]+32 (Subscript[L, 5]) (Subscript[M, K]^2)-32 (Subscript[L, 5]) (Subscript[M, π]^2)))),((Subscript[F, η])^(ii_))->(F^ii)(1+(1/(3 (F^2)))((mo^2) ii(3 A[(Subscript[M, K]^2)]-3 A[(Subscript[M, π]^2)]+16 (Subscript[L, 5]) ((Subscript[M, K]^2)-(Subscript[M, π]^2))))),(Subscript[F, π])->F,(Subscript[F, K])->F(1+(1/(8 (F^2)))((mo^2) (2 A[(Subscript[M, K]^2)]-5 A[(Subscript[M, π]^2)]+3 A[(Subscript[M, η]^2)]+32 (Subscript[L, 5]) (Subscript[M, K]^2)-32 (Subscript[L, 5]) (Subscript[M, π]^2)))),(Subscript[F, η])->F(1+(1/(3 (F^2)))((mo^2)(3 A[(Subscript[M, K]^2)]-3 A[(Subscript[M, π]^2)]+16 (Subscript[L, 5]) ((Subscript[M, K]^2)-(Subscript[M, π]^2)))))}],2,mo],{mo->1,Mp->M}]]];

If[onlyloops==1,(Subscript[L, i_])=0];

(* lagtreep4a, lagtreep4b and lagtreep4c are the order (p^4) tree parts of the functional. *)

lagtreep4a=If[onlyloops==0,If[stcheck==0 && wpar==0,eom[lagp4[x,4,u,U]+wzw[x,4,u,U],6],If[stcheck==0 ,getord[setzero[wzw[x,5,ur,Ur]],5]+getord[setzero[wzw[x,6,ur,Ur]],6]+getord[setzero[lagp4[x,5,ur,Ur]],5]+getord[setzero[lagp4[x,6,ur,Ur]],6]]+Sum[getord[eombilin[lagp4[x,4,uw[jj],Uw[jj]],6],6]/Length[nonexmeson],{jj,Length[nonexmeson]}]+If[stcheck==1,Sum[Sum[eom[Coefficient[lagwp48[x,4,3,uw[jj],Uw[jj]],Nlist[[ii]]]*Nlist[[ii]],6]/Length[nonexmeson],{ii,Length[Nlist]}]+Sum[eom[Coefficient[lagwp427[x,4,3,uw[jj],Uw[jj]],Rlist[[ii]]]*Rlist[[ii]],6]/Length[nonexmeson],{ii,Length[Rlist]}],{jj,Length[nonexmeson]}]+setzero[lagwp4[x,6,5,ur,Ur]],0]],0];

lagtreep4b=If[onlyloops==0,If[stcheck==1 || wpar==1,Sum[eom[Sum[Coefficient[lagp4[x,4,uw[jj],Uw[jj]],(Subscript[L, ii])]*(Subscript[L, ii]),{ii,10}],6]/Length[nonexmeson],{jj,Length[nonexmeson]}]+setzero[Sum[(addbilin[lagp4[x,particleorder,uw[jj],Uw[jj]],particleorder]+addbilin[wzw[x,particleorder,uw[jj],Uw[jj]],particleorder]+eom[wzw[x,particleorder,uw[jj],Uw[jj]],particleorder+2])/Length[nonexmeson],{jj,Length[nonexmeson]}]],0],0];

lagtreep4c=If[onlyloops==0,If[wpar==1,Sum[getord[eombilin[lagp2[x,4,uw[jj],Uw[jj]],6],6]/Length[nonexmeson],{jj,Length[nonexmeson]}],0],0];

If[stcheck==1,fieldtransform[Sum[Simplify[Coefficient[Sum[lagwp48[x,2,2,uw[jj],Uw[jj]]/Length[nonexmeson],{jj,Length[nonexmeson]}],Nlist[[ii]]]]*Nlist[[ii]],{ii,Length[Nlist]}]+Sum[Simplify[Coefficient[Sum[lagwp427[x,2,2,uw[jj],Uw[jj]]/Length[nonexmeson],{jj,Length[nonexmeson]}],Rlist[[ii]]]]*Rlist[[ii]],{ii,Length[Rlist]}],ϕwp4t,φp4t]]

If[stcheck==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=funcchange[trans[Simplify[setzero[dephi[(1/(f^4))*(ReplaceAll[ReplaceAll[insertbilin[lagp2[var,4,u,U],4],{φ->ϕwp4t}],dete->de]),ii]]]]];(Subscript[Φwp4t, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->erg}],{ii,1,8}]]]

treebilin=ReplaceAll[If[stcheck==1,replacements[transform[getamp[strongp2,delta,Φwp4t,φp4t]]],0],{(Subscript[F, π])->F,(Subscript[F, K])->F,(Subscript[F, η])->F}];

(* In the following the tree graphs coming from terms proportional to one pseudoscalar meson field are calculated. *)

monophi6=If[stcheck==1,Simplify[Coefficient[lagwp4[x,1,1,u,U],(Subscript[ϕ, 6])[x]]/(Subscript[M, K]^2)],0];

monophi7=If[stcheck==1,Simplify[Coefficient[lagwp4[x,1,1,u,U],(Subscript[ϕ, 7])[x]]/(Subscript[M, K]^2)],0];

If[stcheck==1,fieldtransformmono[φmt]]

monopret=If[stcheck==1,Sum[eombilin[lagp2[x,6,uw[jj],Uw[jj]],6]/Length[nonexmeson],{jj,Length[nonexmeson]}],0];

monotrees=If[stcheck==1,Simplify[replacements[transform[getamp[monopret,delta,Φs,φmt]]]],0];

(* treep4 is the tree amplitude from the order (p^4) Lagrangian. *)

treep4=If[onlyloops==0,etakorr[ReplaceAll[Sum[simp[replacements[transform[epsorder[getamp[Coefficient[lagtreep4a,constlist[[ii]]]*constlist[[ii]],delta,Φs,If[wpar==1,φsp2,φp2]]+If[stcheck==1|| wpar==1,getamp[Coefficient[lagtreep4b,constlist[[ii]]]*constlist[[ii]],delta,If[wpar==1,Φsa,Φw],If[wpar==1,φsp2,φp2]]+getamp[Coefficient[lagtreep4c,constlist[[ii]]]*constlist[[ii]],delta,Φs,φsp4],0]]]]],{ii,Length[constlist]}]+simp[replacements[transform[epsorder[getamp[ReplaceAll[lagtreep4a+lagtreep4b,Table[constlist[[ii]]->0,{ii,Length[constlist]}]],delta,If[stcheck==0,If[wpar==1,Φs,0],Φw],If[stcheck==0,If[wpar==1,φsp2,0],φp2]]]]]]+monotrees+treebilin,{F->(Subscript[F, π])}]],0];

collect[ject_]:=Do[worklist=Union[constlist,{(OverBar[A])[__],(OverBar[B])[__],C[__]}];tgp=ject;Do[tgp=ReplaceAll[tgp,{worklist[[jj]]->0}],{jj,Length[worklist]}];Return[Sum[simp[etakorr[Coefficient[ject,constlist[[jj]]]]]*constlist[[jj]],{jj,Length[constlist]}]+simp[tgp]]]

simpfast[ject_]:=If[phi6==0,Simplify[ject],ReplaceAll[ex2ordns2[Collect[FactorTerms[ReplaceAll[ject,{(Subscript[F, π])->F}]],Power[__,-1]],8,F],{F->(Subscript[F, π])}]]

treeresult4temp=final[collect[simpfast[treep4+ReplaceAll[renomp4,{A[_]->0,F->(Subscript[F, π])}]]]];

treeresult4=normalization[treeresult4temp];

treeresult2and4=normalization[treeresult2temp+treeresult4temp];

If[onlytreep2and4==1,result=treeresult2and4;title="Complete tree level amplitude to O(p^4):";loopresult="Not evaluated";output[result]];

If[onlytreep2and4==1,Quit[]]

(* *)

(* The matrices (Subscript[(OverHat[λ]), i]) are needed to express ϕ in the spherical basis of SU(3). *)

(Subscript[(OverHat[λ]), 1])=(1/(Sqrt[2]))((Subscript[λ, 1])+I (Subscript[λ, 2]));

(Subscript[(OverHat[λ]), 2])=(1/(Sqrt[2]))((Subscript[λ, 1])-I (Subscript[λ, 2]));

(Subscript[(OverHat[λ]), 3])=(Subscript[λ, 3]);

(Subscript[(OverHat[λ]), 4])=(1/(Sqrt[2]))((Subscript[λ, 4])+I (Subscript[λ, 5]));

(Subscript[(OverHat[λ]), 5])=(1/(Sqrt[2]))((Subscript[λ, 4])-I (Subscript[λ, 5]));

(Subscript[(OverHat[λ]), 6])=(1/(Sqrt[2]))((Subscript[λ, 6])+I (Subscript[λ, 7]));

(Subscript[(OverHat[λ]), 7])=(1/(Sqrt[2]))((Subscript[λ, 6])-I (Subscript[λ, 7]));

(Subscript[(OverHat[λ]), 8])=(Subscript[λ, 8]);

hc[ject_]:=Transpose[Conjugate[ject]]

(* gam[var,ind,aa,bb,cp,ord,u,U] and sig[var,ind,aa,bb,cp,ord,u,U] are the components of the matrices that occur in the strong part of the generating functional. *)

gam[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=-(1/2)Tr2ord[kom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].gamma[var,ind,ord,u,U],ord]

gamma[var_,ind_,ord_,u_,U_]:=(1/2)kom[(Subscript[(OverHat[u]), ord-1])[var],de[(Subscript[u, ord])[var],var,ind]]-(1/2)I (Subscript[(OverHat[u]), ord])[var].r[var,ind].(Subscript[u, ord])[var]-(1/2)I (Subscript[u, ord])[var].l[var,ind].(Subscript[(OverHat[u]), ord])[var]

sighat[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/2)Tr2ord[(kom[(Subscript[(OverHat[λ]), aa]),((1/2)(Subscript[(OverHat[u]), ord-2])[var].d[(Subscript[α, var]),U,ord-1,var].(Subscript[(OverHat[u]), ord-2])[var])].kom[hc[(Subscript[(OverHat[λ]), bb])],((1/2)(Subscript[(OverHat[u]), ord-2])[var].d[(Subscript[α, var]),U,ord-1,var].(Subscript[(OverHat[u]), ord-2])[var])]),ord]+(1/4)Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].sigma[var,ord,u,U],ord]

sig[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=sighat[var,ind,aa,bb,cp,ord,u,U]-(1/4)Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].sigma[var,0,u,U],ord]

sigma[var_,ord_,u_,U_]:=(1/2)((Subscript[(OverHat[u]), ord])[var].χ[var].(Subscript[(OverHat[u]), ord])[var]+(Subscript[u, ord])[var].(OverHat[χ])[var].(Subscript[u, ord])[var])

Y[var_,ind_,ord_,u_,U_]:=(1/2)(Subscript[(OverHat[u]), ord-1])[var].d[ind,U,ord,var].(Subscript[(OverHat[u]), ord-1])[var]

(* In the following the projectors needed for the nonleptonic weak loop part are defined. *)

k[1]=(4/(F^2))(Subscript[G, 8]) (F^4);P[1]=Q[2,3];k[-1]=(4/(F^2))(Subscript[(OverHat[G]), 8]) (F^4);P[-1]=Q[3,2];

PR[1]:={{(1/2)(Subscript[G, 27]) (F^4)*Q[2,3],Q[1,1]},{(1/3) (Subscript[G, 27]) (F^4)*Q[2,1],Q[1,3]}}

PR[-1]:={{(1/2)(Subscript[(OverHat[G]), 27]) (F^4)*Q[3,2],Q[1,1]},{(1/3) (Subscript[(OverHat[G]), 27]) (F^4)*Q[1,2],Q[3,1]}}

R[cp_]:=Join[PR[cp],Table[Reverse[PR[cp][[ii]]],{ii,1,2}]]

Q[1,1]={{1,0,0},{0,0,0},{0,0,0}};

Q[1,2]={{0,0,0},{1,0,0},{0,0,0}};

Q[3,1]={{0,0,1},{0,0,0},{0,0,0}};

Q[2,1]={{0,1,0},{0,0,0},{0,0,0}};

Q[1,3]={{0,0,0},{0,0,0},{1,0,0}};

(* Tr2ord[ject, ord] evaluates the trace of a matrix up to order ord in the external fields. *)

Tr2ord[ject_,ord_]:=ex2ord[dot0[Tr[ject]],ord]

q2ord2[ject1_,ject2_,cp_,ord_]:=ex2ord[(4/(F^2))*Sum[ex2ord[Tr2ord[ject1,ord]Tr2ord[ject2,ord],ord],{si,4}],ord]

ex2ord2[ject_,ord_]:=ex2ord[dot0[ject],ord]

(* gamw[var,ind,aa,bb,cp,ord,u,U], sigw1[var,ind,aa,bb,cp,ord,u,U],, sigw2[var,ind,aa,bb,cp,ord,u,U], sigw2[var,ind,aa,bb,cp,ord,u,U], sigw4[var,ind,aa,bb,cp,ord,u,U] are the components of the matrices that occur in the weak part of the generating functional. NP, alph, dalph,alphsub and sigrotmat are their subfunctions. *)

NP[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/4)k[cp]Tr2ord[kom[(Subscript[u, ord-1])[var].P[cp].(Subscript[(OverHat[u]), ord-1])[var],Y[var,ind,ord,u,U]].akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]],ord]+(1/2)q2ord2[(Subscript[u, ord-1])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-1])[var].(Subscript[(OverHat[λ]), aa]),kom[(Subscript[u, ord-1])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord-1])[var],Y[var,ind,ord,u,U]].hc[(Subscript[(OverHat[λ]), bb])],cp,ord]+(1/2)q2ord2[(Subscript[u, ord-1])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-1])[var].hc[(Subscript[(OverHat[λ]), bb])],kom[(Subscript[u, ord-1])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord-1])[var],Y[var,ind,ord,u,U]].(Subscript[(OverHat[λ]), aa]),cp,ord]

alph[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/4)k[cp]Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].(Subscript[u, ord])[var].P[cp].(Subscript[(OverHat[u]), ord])[var],ord]+(1/2)q2ord2[(Subscript[(OverHat[λ]), aa]).(Subscript[u, ord])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord])[var],hc[(Subscript[(OverHat[λ]), bb])].(Subscript[u, ord])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord])[var],cp,ord]

dalph[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=de[alph[var,0,aa,bb,cp,ord,u,U],var,ind]

alphsub[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/4)k[cp]Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].(Subscript[u, ord])[var].P[cp].(Subscript[(OverHat[u]), ord])[var],ord]+(1/2)q2ord2[(Subscript[(OverHat[λ]), aa]).(Subscript[u, ord])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord])[var],hc[(Subscript[(OverHat[λ]), bb])].(Subscript[u, ord])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord])[var],cp,ord]-((1/4)k[cp]Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].(Subscript[u, 0])[var].P[cp].(Subscript[(OverHat[u]), 0])[var],0]+(1/2)q2ord2[(Subscript[(OverHat[λ]), aa]).(Subscript[u, 0])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), 0])[var],hc[(Subscript[(OverHat[λ]), bb])].(Subscript[u, 0])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), 0])[var],cp,0])

rotmat=Table[ReplaceAll[Coefficient[(Subscript[φp2, ii])[ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var],(Subscript[e, jj])[var]],{ϕ1->1,ϕ2->1,ϕ3->1,ϕ4->1,ϕ5->1,ϕ6->1,ϕ7->1,ϕ8->1}],{ii,8},{jj,8}];

trotmat=Transpose[rotmat];

rotmat2=Table[Simplify[(1/2)(rotmat[[ii,jj]]-rotmat[[ap[[jj]],ap[[ii]]]])],{ii,8},{jj,8}];

trotmat2=Transpose[rotmat2];

sigrotmat[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=Simplify[Sum[sig[var,ind,aa,ii,cp,ord,u,U]*rotmat[[ap[[ii]],ap[[bb]]]],{ii,8}]]+Simplify[Sum[trotmat[[aa,ii]]*sig[var,ind,ii,bb,cp,ord,u,U],{ii,8}]]

gamrotmat[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=Simplify[Sum[gam[var,ind,aa,ii,cp,ord,u,U]*rotmat2[[ap[[ii]],ap[[bb]]]],{ii,8}]]+Simplify[Sum[trotmat2[[aa,ii]]*gam[var,ind,ii,bb,cp,ord,u,U],{ii,8}]]

gamw[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=-(1/2)(-(1/4)k[cp]Tr2ord[akom[(Subscript[u, ord-1])[var].P[cp].(Subscript[(OverHat[u]), ord-1])[var],Y[var,ind,ord,u,U]].kom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]],ord]-(1/2)k[cp]Tr2ord[(Subscript[u, ord-1])[var].P[cp].(Subscript[(OverHat[u]), ord-1])[var].(Subscript[(OverHat[λ]), aa]).Y[var,ind,ord,u,U].hc[(Subscript[(OverHat[λ]), bb])]-(Subscript[u, ord-1])[var].P[cp].(Subscript[(OverHat[u]), ord-1])[var].hc[(Subscript[(OverHat[λ]), bb])].Y[var,ind,ord,u,U].(Subscript[(OverHat[λ]), aa]),ord]-q2ord2[kom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].(Subscript[u, ord])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord])[var],(Subscript[u, ord])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord])[var].Y[var,ind,ord,u,U],cp,ord]-(1/2)q2ord2[(Subscript[u, ord])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord])[var].(Subscript[(OverHat[λ]), aa]),kom[(Subscript[u, ord])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord])[var],Y[var,ind,ord,u,U]].hc[(Subscript[(OverHat[λ]), bb])],cp,ord]+(1/2)q2ord2[(Subscript[u, ord])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord])[var].hc[(Subscript[(OverHat[λ]), bb])],kom[(Subscript[u, ord])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord])[var],Y[var,ind,ord,u,U]].(Subscript[(OverHat[λ]), aa]),cp,ord])+gamrotmat[var,ind,aa,bb,cp,ord,u,U]

sigw1[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=-(1/4)k[cp]Tr2ord[akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]].akom[Y[var,(Subscript[α, var]),ord-1,u,U].Y[var,(Subscript[α, var]),ord-1,u,U],(Subscript[u, ord-2])[var].P[cp].(Subscript[(OverHat[u]), ord-2])[var]],ord]+(1/4)k[cp]Tr2ord[(Subscript[u, ord-2])[var].P[cp].(Subscript[(OverHat[u]), ord-2])[var].(Subscript[(OverHat[λ]), aa]).Y[var,(Subscript[α, var]),ord-1,u,U].Y[var,(Subscript[α, var]),ord-1,u,U].hc[(Subscript[(OverHat[λ]), bb])]+(Subscript[u, ord-2])[var].P[cp].(Subscript[(OverHat[u]), ord-2])[var].hc[(Subscript[(OverHat[λ]), bb])].Y[var,(Subscript[α, var]),ord-1,u,U].Y[var,(Subscript[α, var]),ord-1,u,U].(Subscript[(OverHat[λ]), aa]),ord]+(1/4)k[cp]Tr2ord[akom[(Subscript[u, ord-2])[var].P[cp].(Subscript[(OverHat[u]), ord-2])[var],Y[var,(Subscript[α, var]),ord-1,u,U]].((Subscript[(OverHat[λ]), aa]).Y[var,(Subscript[α, var]),ord-1,u,U].hc[(Subscript[(OverHat[λ]), bb])]+hc[(Subscript[(OverHat[λ]), bb])].Y[var,(Subscript[α, var]),ord-1,u,U].(Subscript[(OverHat[λ]), aa])),ord]-(1/4)k[cp]Tr2ord[Y[var,(Subscript[α, var]),ord-1,u,U].(Subscript[u, ord-2])[var].P[cp].(Subscript[(OverHat[u]), ord-2])[var].Y[var,(Subscript[α, var]),ord-1,u,U].akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]],ord]-(1/2)q2ord2[kom[ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-2])[var],ord-2],ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1]].(Subscript[(OverHat[λ]), aa]),kom[ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord-2])[var],ord-2],ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1]].hc[(Subscript[(OverHat[λ]), bb])],cp,ord]-(1/2)q2ord2[ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-2])[var],ord-2].ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1],akom[ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord-2])[var],ord-2],ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1]].akom[(Subscript[(OverHat[λ]), aa]),hc[(Subscript[(OverHat[λ]), bb])]],cp,ord]+q2ord2[ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-2])[var],ord-2].(Subscript[(OverHat[λ]), aa]).ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1].hc[(Subscript[(OverHat[λ]), bb])]+ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,1]].(Subscript[(OverHat[u]), ord-2])[var],ord-2].hc[(Subscript[(OverHat[λ]), bb])].ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1].(Subscript[(OverHat[λ]), aa]),ex2ord2[(Subscript[u, ord-2])[var].R[cp][[si,2]].(Subscript[(OverHat[u]), ord-2])[var],ord-2].ex2ord2[Y[var,(Subscript[α, var]),ord-1,u,U],ord-1],cp,ord]+sigrotmat[var,ind,aa,bb,cp,ord,u,U]

sigw2[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/2)(de[NP[var,(Subscript[α, var]),aa,bb,cp,ord,u,U],var,(Subscript[α, var])]+Sum[gam[var,(Subscript[α, var]),aa,ii,cp,ord-1,u,U]*NP[var,(Subscript[α, var]),ii,bb,cp,ord-2,u,U]-NP[var,(Subscript[α, var]),aa,ii,cp,ord-2,u,U]*gam[var,(Subscript[α, var]),ii,bb,cp,ord-1,u,U],{ii,1,8}])

sigw3[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=(1/2)(-ex2ord[de[dalph[var,(Subscript[α, var]),aa,bb,cp,ord,u,U]-mp2[alph,gam,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord-2,ord]+mp2[gam,alph,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord,ord-2],var,(Subscript[α, var])],ord]+ex2ord[mp2[dalph,gam,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord-2,ord]-mp2[gam,dalph,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord,ord-2],ord]+ex2ord[2mp3[gam,alph,gam,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord-2,ord-2,ord-2]-mp3[gam,gam,alph,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord-2,ord-2,ord-2]-mp3[alph,gam,gam,var,(Subscript[α, var]),aa,bb,cp,ord,u,U,ord-2,ord-2,ord-2],ord])

sigw4[var_,ind_,aa_,bb_,cp_,ord_,u_,U_]:=-(1/2)(ex2ord[Sum[sig[var,(Subscript[α, var]),aa,ii,cp,ord,u,U]*alphsub[var,(Subscript[α, var]),ii,bb,cp,ord-2,u,U]+alphsub[var,(Subscript[α, var]),aa,ii,cp,ord-2,u,U]*sig[var,(Subscript[α, var]),ii,bb,cp,ord,u,U],{ii,1,8}],ord]+ex2ord[Sum[sighat[var,(Subscript[α, var]),aa,ii,cp,0,u,U]*alphsub[var,(Subscript[α, var]),ii,bb,cp,ord,u,U]+alphsub[var,(Subscript[α, var]),aa,ii,cp,ord,u,U]*sighat[var,(Subscript[α, var]),ii,bb,cp,0,u,U],{ii,1,8}],ord])

(* mp2[ject1, ject2, var, ind, aa, bb, cp, ord, u, U, ord1, ord2] and mp3[ject1, ject2, var, ind, aa, bb, cp, ord, u, U, ord1, ord2] perform the matrix product of two/three matrices. *)

mp2[ject1_,ject2_,var_,ind_,aa_,bb_,cp_,ord_,u_,U_,ord1_,ord2_]:=ex2ord[Sum[ex2ord[ject1[var,(Subscript[α, var]),aa,ii,cp,ord1,u,U]*ject2[var,(Subscript[α, var]),ii,bb,cp,ord2,u,U],ord],{ii,1,8}],ord]

mp3[ject1_,ject2_,ject3_,var_,ind_,aa_,bb_,cp_,ord_,u_,U_,ord1_,ord2_,ord3_]:=ex2ord[Sum[ex2ord[ject1[var,(Subscript[α, var]),aa,ii,cp,ord1,u,U]*ject2[var,(Subscript[α, var]),ii,jj,cp,ord2,u,U]*ject3[var,(Subscript[α, var]),jj,bb,cp,ord3,u,U],ord],{ii,1,8},{jj,1,8}],ord]

(* ord3n4[ject], ord5n6[ject] and ord3n6[ject] provide the parts of the matrices with the corresponding number of external fields. *)

ord3n4[ject_]:=getordns[ject,3]+getordns[ject,4]

ord5n6[ject_]:=getordns[ject,5]+getordns[ject,6]

ord3to6[ject_]:=getordns[ject,6]-getordns[ject,2]

(* *)

(* matrix2[ject,select] and matrix4[ject,select] are used to table the matrices up to order (f^4). *)

matrix2[ject_,cp_]:=Do[tempmat=Table[0,{ii,1,8},{jj,1,8}];Do[tempmat[[ii,jj]]=If[Or[(stcheck==0&& wpar==0),MemberQ[matrixlistw,ject]],setzero[ject[x,μ,ii,jj,cp,1,ur,Ur]+getord[ject[x,μ,ii,jj,cp,2,ur,Ur],2]],setzero[ject[x,μ,ii,jj,cp,2,ur,Ur]]+setzero[Sum[addbilin[ject[x,μ,ii,jj,cp,2,uw[kk],Uw[kk]],2]/Length[nonexmeson],{kk,Length[nonexmeson]}]]],{ii,1,8},{jj,1,8}];Return[tempmat];]

matrix4[ject_,cp_]:=Do[tempmat=Table[0,{ii,1,8},{jj,1,8}];
Do[tempmat[[ii,jj]]=If[Or[(stcheck==0&& wpar==0),MemberQ[matrixlistw,ject]],Sum[eom[ject[x,μ,ii,jj,cp,1,uw[nn],Uw[nn]]+getord[ject[x,μ,ii,jj,cp,2,uw[nn],Uw[nn]],2],4]/Length[nonexmeson],{nn,Length[nonexmeson]}]+getordns[setzero[ject[x,μ,ii,jj,cp,3,ur,Ur]],3]+getordns[setzero[ject[x,μ,ii,jj,cp,4,ur,Ur]],4],getord[setzero[ject[x,μ,ii,jj,cp,4,ur,Ur]],4]+setzero[Sum[addbilin[ex2ord[ject[x,μ,ii,jj,0,4,uw[kk],Uw[kk]],4],4]/Length[nonexmeson],{kk,Length[nonexmeson]}]]+Sum[eom[getord[ject[x,μ,ii,jj,0,2,uw[kk],Uw[kk]],2]/Length[nonexmeson],4],{kk,Length[nonexmeson]}]+Sum[eombilin[getord[ject[x,μ,ii,jj,0,2,uw[kk],Uw[kk]],2],4]/Length[nonexmeson],{kk,Length[nonexmeson]}]],{ii,1,8},{jj,1,8}];Return[tempmat];]

(* oneweak[ject] ensures in case of weak processes that there is only one weak perturbation. *)

oneweak[ject_]:=If[stcheck==1,ReplaceAll[ex2ordns[ject,1,wp],{wp->0}],ject]

(* termtable2[na,terml,funcl,factorl], termtable3[na,terml,funcl,factorl] and their subfunctions calculate the products of the matrices and perform the functional differentation. *)

termtable2[na_,terml_,funcl_,factorl_,ΦΦ_,φ_]:=Do[subtermtable2[na[[ii]],terml,funcl,factorl,ii,ΦΦ,φ],{ii,Length[terml]}]

subtermtable2[namel_,terml_,funcl_,factorl_,ll_,ΦΦ_,φ_]:=Do[namel[ii,jj]=If[Length[terml[[ll]]]==2,factorl[[ll]]*getamp[oneweak[If[particleorder<5,(1/2),1]*(mat4[terml[[ll,1]]][[jj,ii]]*second[mat2[terml[[ll,2]]][[ii,jj]]]+mat2[terml[[ll,1]]][[jj,ii]]*second[mat4[terml[[ll,2]]][[ii,jj]]])],funcl[[ll]],ΦΦ,φ,((Subscript[M, ii])^2),((Subscript[M, jj])^2),0],getamp[oneweak[Sum[factorl[[ll]]*special[mat2[terml[[ll,1]]][[jj,kk]]]*special[mat2[terml[[ll,2]]][[kk,ii]]]*second[mat2[terml[[ll,3]]][[ii,jj]]],{kk,1,8}]],funcl[[ll]],ΦΦ,φ,((Subscript[M, ii])^2),((Subscript[M, jj])^2),0]],{ii,1,8},{jj,1,8}]

termtable3[na_,terml_,funcl_,factorl_,φ_]:=Do[subtermtable3[na[[ii]],terml,funcl,factorl,ii,φ],{ii,Length[terml]}]

subtermtable3[namel_,terml_,funcl_,factorl_,lll_,φ_]:=Do[namel[iii,jjj,kkk]=factorl[[lll]]*getamp[oneweak[mat2[terml[[lll,1]]][[jjj,iii]]*second[mat2[terml[[lll,2]]][[iii,kkk]]]*third[mat2[terml[[lll,3]]][[kkk,jjj]],lll,terml]],funcl[[lll]],0,φ,((Subscript[M, iii])^2),((Subscript[M, jjj])^2),((Subscript[M, kkk])^2)],{iii,1,8},{jjj,1,8},{kkk,1,8}]

(* *)

(* smatrixlist is the list of the originally symmetric matrices. *)

smatrixlist=If[stcheck==0,{sig},{sig,sigw1,sigw2,sigw3,sigw4}];

(* amatrixlist is the list of the originally antisymmetric matrices. *)

amatrixlist=If[stcheck==0,{gam},{gam,gamw}];

matrixlist=Union[smatrixlist,amatrixlist];

(* smatrixlistw and amatrixlistw are the lists of the matrices due to nonleptonic weak interactions. *)

smatrixlistw={sigw1,sigw2,sigw3,sigw4};

amatrixlistw={gamw};

matrixlistw=Union[smatrixlistw,amatrixlistw];

(* Now the matrices occurring in the generating functional are calculated. *)

mattable[matrixlist]

(* *)

Remove[term,terma,termb]

(* The lists namelist2, termlist2, funclist2 and factorlist2 contain the information that is needed to calculate the two-propagator part of the loop amplitude. *)

namelist2={Table[(Subscript[term, ii]),{ii,1,4}],If[stcheck==1,Table[(Subscript[term, ii]),{ii,5,20}],{}]};

namelist2a={If[stcheck==1 || wpar==1,Table[(Subscript[terma, ii]),{ii,1,4}],{}]};

namelist2b={If[stcheck==1 || wpar==1,Table[(Subscript[termb, ii]),{ii,1,4}],{}]};

termlist2={{{gam,gam},{gam,sig},{sig,sig},{gam,gam,sig}},If[stcheck==1,{{gamw,gam},{gamw,sig},{gam,sigw1},{gam,sigw2},{gam,sigw3},{gam,sigw4},{sigw1,sig},{sigw2,sig},{sigw3,sig},{sigw4,sig},{gam,gam,sigw1},{gam,gam,sigw2},{gam,gam,sigw3},{gam,gam,sigw4},{gamw,gam,sig},{gam,gamw,sig}},{}]};

funclist2={{TG1,TG2,TG3,TG3},If[stcheck==1,{TG1,TG2,TG2,TG2,TG2,TG2,TG3,TG3,TG3,TG3,TG3,TG3,TG3,TG3,TG3,TG3},{}]};

factorlist2={{1,1,1,2},If[stcheck==1,{2,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2},{}]};

(* With termtable2 the different parts of the amplitude with two propagators in the loop, called term1[i,j], term2[i,j], ..., are calculated. *)

If[stcheck==0 ,termtable2[namelist2[[1]],termlist2[[1]],funclist2[[1]],factorlist2[[1]],If[wpar==1,Φsa,Φs],If[wpar==1,φsp2,0]]]

If[stcheck==1 ,termtable2[namelist2a[[1]],termlist2[[1]],funclist2[[1]],factorlist2[[1]],Φw,0];termtable2[namelist2b[[1]],termlist2[[1]],funclist2[[1]],factorlist2[[1]],Φs,φp2];Do[Do[(Subscript[term, ii])[jj,kk]=(Subscript[terma, ii])[jj,kk]+(Subscript[termb, ii])[jj,kk],{ii,4}],{jj,8},{kk,8}]]

If[stcheck==1,termtable2[namelist2[[2]],termlist2[[2]],funclist2[[2]],factorlist2[[2]],Φs,0]]

(* The terms are added up: *)

Do[(Subscript[expression, ll])=replacements[transform[epscorr[Sum[(Subscript[term, ll])[ii,jj],{ii,1,8},{jj,1,8}]]]],{ll,If[stcheck==0,4,20]}];

(* The part tadpoles2 is calculated by use of (Subscript[term, 1])[i,j] or (Subscript[term, 5])[i,j]: *)

tadpoles2=replacements[transform[Sum[ReplaceAll[If[stcheck==0,(Subscript[term, 1])[ii,jj],(Subscript[term, 1])[ii,jj]+(Subscript[term, 5])[ii,jj]],{TG1[x_,y_,u_,v_,w_,ms1_,ms2_,0]->(1/4)(A[ms1]+A[ms2])*sp[u,v]}],{ii,1,8},{jj,1,8}]]];

(* The lists namelist3, termlist3, funclist3 and factorlist3 contain the information that is needed to calculate the two-propagator part of the loop amplitude. *)

namelist3={Table[(Subscript[term, ii]),{ii,If[stcheck==0,5,21],If[stcheck==0,8,24]}],If[stcheck==1,Table[(Subscript[term, ii]),{ii,25,39}],{}]};

termlist3={{{gam,gam,gam},{gam,gam,sig},{sig,sig,gam},{sig,sig,sig}},{{gamw,gam,gam},{gamw,gam,sig},{gam,gam,sigw1},{gam,gam,sigw2},{gam,gam,sigw3},{gam,gam,sigw4},{sigw1,sig,gam},{sigw2,sig,gam},{sigw3,sig,gam},{sigw4,sig,gam},{sig,sig,gamw},{sigw1,sig,sig},{sigw2,sig,sig},{sigw3,sig,sig},{sigw4,sig,sig}}};

funclist3={{TG4,TG5,TG6,TG7},{TG4,TG5,TG5,TG5,TG5,TG5,TG6,TG6,TG6,TG6,TG6,TG7,TG7,TG7,TG7}};

factorlist3={{1,1,1,1},{3,2,1,1,1,1,2,2,2,2,1,3,3,3,3}};

(* With termtable3 the different parts of the amplitude with three propagators in the loop are calculated. *)

If[particleorder>4,termtable3[namelist3[[1]],termlist3[[1]],funclist3[[1]],factorlist3[[1]],If[stcheck==0,0,If[wpar==1,φsp2,φp2]]]]

If[particleorder>4,If[stcheck==1,termtable3[namelist3[[2]],termlist3[[2]],funclist3[[2]],factorlist3[[2]],0]]]

(* The terms are added up: *)

Do[(Subscript[expression, ll])=If[particleorder>4,replacements[transform[epscorr[Sum[(Subscript[term, ll])[ii,jj,kk],{ii,1,8},{jj,1,8},{kk,1,8}]]]],0],{ll,If[stcheck==0,5,21],If[stcheck==0,8,39]}];

(* *)

(* In the following the trees of order (p^4) are generated. *)

If[stcheck==1,fieldtransform[tadep[sigw1,2]+tadep[sigw2,2]+tadep[sigw3,2]+tadep[sigw4,2],ϕwp4l,φp4l]]

If[stcheck==1,Do[ap={2,1,3,5,4,7,6,8};Do[erg=funcchange[trans[Simplify[setzero[dephi[(1/(f^4))*(ReplaceAll[ReplaceAll[insertbilin[lagp2[var,4,u,U],4],{φ->ϕwp4l}],dete->de]),ii]]]]];(Subscript[Φwp4l, ap[[ii]]])=ReplaceAll[Function[{ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6,ϕ7,ϕ8,γ,W1,W2,var},eer],{eer->erg}],{ii,1,8}]]]

(* tadbilin contains the amplitude that originates from bilinear tadpole terms. *)

tadbilin=ReplaceAll[If[stcheck==1,replacements[transform[getamp[strongp2,delta,Φwp4l,φp4l]]],0],{(Subscript[F, π])->F,(Subscript[F, K])->F,(Subscript[F, η])->F}];

(* In the following the part of the amplitude that includes contributions linear in a pseudoscalar meson field are calculated. *)

monotpi=If[stcheck==1,Simplify[momzero[Sum[sigw1[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw2[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw3[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw4[x,(Subscript[α, x]),ii,ii,stp,1,u,U],{ii,3}]*(1/2)A[(Subscript[M, π]^2)]]],0];

monotkaon=If[stcheck==1,Simplify[momzero[Sum[sigw1[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw2[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw3[x,(Subscript[α, x]),ii,ii,stp,1,u,U]+sigw4[x,(Subscript[α, x]),ii,ii,stp,1,u,U],{ii,4,7}]*(1/2)A[(Subscript[M, K]^2)]]],0];

monoteta=If[stcheck==1,Simplify[momzero[(sigw1[x,(Subscript[α, x]),8,8,stp,1,u,U]+sigw2[x,(Subscript[α, x]),8,8,stp,1,u,U]+sigw3[x,(Subscript[α, x]),8,8,stp,1,u,U]+sigw4[x,(Subscript[α, x]),8,8,stp,1,u,U])*(1/2)A[(Subscript[M, η]^2)]]],0];

monophi6=If[stcheck==1,Simplify[Coefficient[monotpi+monotkaon+monoteta,(Subscript[ϕ, 6])[x]]/(Subscript[M, K]^2)],0];

monophi7=If[stcheck==1,Simplify[Coefficient[monotpi+monotkaon+monoteta,(Subscript[ϕ, 7])[x]]/(Subscript[M, K]^2)],0];

If[stcheck==1,fieldtransformmono[φml]]

monoprel=If[stcheck==1,Sum[eombilin[lagp2[x,6,uw[jj],Uw[jj]],6]/Length[nonexmeson],{jj,Length[nonexmeson]}],0];

monotadpoles=If[stcheck==1,Simplify[replacements[transform[getamp[monoprel,delta,Φs,φml]]]],0];

(* In the following the tadpole part of the amplitude is generated. *)

(Subscript[tadpoles1a, 1])=If[stcheck==1||wpar==1,replacements[transform[Sum[getamp[setzero[Sum[eombilin[getord[sig[x,0,ii,ii,0,4,uw[jj],Uw[jj]],4]/Length[nonexmeson],6],{jj,Length[nonexmeson]}]],delta,Φs,If[wpar==1,φsp2,φp2]]*(1/2)A[((Subscript[M, ii])^2)],{ii,8}]]],0];

(Subscript[tadpoles1a, 2])=If[stcheck==1,replacements[transform[Sum[Sum[getamp[eom[ord3n4[smatrixlistw[[1]][x,0,ii,ii,stp,4,uw[jj],Uw[jj]]],6]/Length[nonexmeson],delta,Φs,φp2],{jj,Length[nonexmeson]}]*(1/2)A[((Subscript[M, ii])^2)],{ii,8}]]],0];

(Subscript[tadpoles1a, 3])=If[stcheck==1,replacements[transform[Sum[Sum[getamp[eom[ord3n4[smatrixlistw[[2]][x,0,ii,ii,stp,4,uw[jj],Uw[jj]]],6]/Length[nonexmeson],delta,Φs,φp2],{jj,Length[nonexmeson]}]*(1/2)A[((Subscript[M, ii])^2)],{ii,8}]]],0];

(Subscript[tadpoles1a, 4])=If[stcheck==1,replacements[transform[Sum[Sum[getamp[eom[ord3n4[smatrixlistw[[3]][x,0,ii,ii,stp,4,uw[jj],Uw[jj]]],6]/Length[nonexmeson],delta,Φs,φp2],{jj,Length[nonexmeson]}]*(1/2)A[((Subscript[M, ii])^2)],{ii,8}]]],0];

(Subscript[tadpoles1a, 5])=If[stcheck==1,replacements[transform[Sum[Sum[getamp[eom[ord3n4[smatrixlistw[[4]][x,0,ii,ii,stp,4,uw[jj],Uw[jj]]],6]/Length[nonexmeson],delta,Φs,φp2],{jj,Length[nonexmeson]}]*(1/2)A[((Subscript[M, ii])^2)],{ii,8}]]],0];

Do[(Subscript[tadpoles1a, 5+kk])=If[stcheck==1,replacements[transform[Sum[getamp[setzero[ord5n6[smatrixlistw[[kk]][x,0,ii,ii,stp,6,ur,Ur]*(1/2)A[((Subscript[M, ii])^2)]]],delta,Φs,φp2],{ii,8}]]],0],{kk,1,4}]

tadpoles1aa=If[stcheck==0 && wpar==0,replacements[transform[Sum[getamp[getord[setzero[sig[x,0,ii,ii,0,6,ur,Ur]],6]*(1/2)A[((Subscript[M, ii])^2)],delta,Φs],{ii,8}]]],tadbilin+monotadpoles+Sum[(Subscript[tadpoles1a, ii]),{ii,9}]+If[wpar==1,replacements[transform[Sum[getamp[ord5n6[setzero[sig[x,0,ii,ii,0,6,ur,Ur]]]*(1/2)A[((Subscript[M, ii])^2)],delta,0],{ii,8}]]],0]];

tadpoles1bb=If[stcheck==0 && wpar==0,replacements[transform[Sum[getamp[eom[getord[sig[x,0,ii,ii,0,4,u,U],4],6]*(1/2)A[((Subscript[M, ii])^2)],delta,Φs],{ii,8}]]],replacements[transform[Sum[getamp[Sum[(eom[ord3n4[sig[x,0,ii,ii,0,4,uw[jj],Uw[jj]]],6]+eom[getord[sig[x,0,ii,ii,0,2,uw[jj],Uw[jj]],1,(Subscript[C, w]) ],4])*(1/2)A[((Subscript[M, ii])^2)]/Length[nonexmeson],{jj,Length[nonexmeson]}]+Sum[addbilin[ord5n6[sig[x,0,ii,ii,0,6,uw[jj],Uw[jj]]],6]*(1/2)A[((Subscript[M, ii])^2)]/Length[nonexmeson],{jj,Length[nonexmeson]}],delta,If[wpar==1,Φsa,Φw],If[wpar==1,φsp2,φp2]],{ii,8}]]]];

tadpoles1=tadpoles1aa+tadpoles1bb;

tadpoles=tadpoles1+tadpoles2;

(* *)

If[norecurs==1,nonre=1,nonre=0,nonre=0];

(Subscript[expressiont, 5])=If[nonre==0,Expand[ReplaceAll[(Subscript[expression, 5]),{TG4->TGG4}]],0];

(Subscript[expressiont, 21])=If[nonre==0 && stcheck==1,Expand[ReplaceAll[(Subscript[expression, 21]),{TG4->TGG4}]],0];

(Subscript[expressiont, 25])=If[nonre==0 && stcheck==1,Expand[ReplaceAll[(Subscript[expression, 25]),{TG4->TGG4}]],0];

momentalist[ject_]:=Do[perlist=Flatten[Table[Sort[Union[Flatten[Table[Partition[Permutations[ject][[ii]],jj],{ii,Length[Permutations[ject]]}],1]]],{jj,Length[ject]}],1];Return[Union[Table[Sum[Sort[perlist[[ii]]][[jj]],{jj,Length[perlist[[ii]]]}],{ii,Length[perlist]}]]]]

momlist2=momentalist[imlist];

propmasslist={((Subscript[M, π])^2),((Subscript[M, K])^2),((Subscript[M, η])^2)};

bflist=If[nonre==0,Union[Flatten[Table[B[sp[momlist2[[ii]],momlist2[[ii]]],propmasslist[[jj]],propmasslist[[kk]]],{ii,Length[momlist2]},{jj,Length[propmasslist]},{kk,Length[propmasslist]}]]],{}];

momlist1=Table[Flatten[{momlist2[[ii]]/.{Plus->List}}],{ii,Length[momlist2]}];

(* *)

Cswitch12[ject_]:=Do[saved1=ject[[1]];saved2=ject[[2]];saved5=ject[[5]];saved6=ject[[6]];Return[replacelist[ject,{saved2,saved1,saved6,saved5},{1,2,5,6}]]]

Cswitch13[ject_]:=Do[saved1=ject[[1]];saved3=ject[[3]];saved4=ject[[4]];saved6=ject[[6]];Return[replacelist[ject,{saved3,saved1,saved6,saved4},{1,3,4,6}]]]

Cswitch[ject_,1]:=Cswitch12[ject]

Cswitch[ject_,2]:=Cswitch13[ject]

Cswitch[ject_,3]:=Cswitch13[Cswitch12[ject]]

Cswitch[ject_,4]:=Cswitch12[Cswitch13[ject]]

Cswitch[ject_,5]:=Cswitch12[Cswitch13[Cswitch12[ject]]]

replacelist[ject_,novel_,inpos_]:=ReplacePart[ject,novel,{{inpos[[1]]},{inpos[[2]]},{inpos[[3]]},{inpos[[4]]}},{{1},{2},{3},{4}}]

Cchange[ject_]:=ReplaceAll[ject,{C[a_,b_,c_,d_,e_,f_]->C[Simplify[a],Simplify[b],Simplify[a+b-2c],d,e,f]}]

Cchange[ject_]:=ReplaceAll[ReplaceAll[ject,{C[a_,b_,c_,d_,e_,f_]->tempfunc[C[a,b,a+b-2c,d,e,f]]}],{tempfunc->Simplify}]

Border[ject_]:=ReplaceAll[ReplaceAll[ReplaceAll[ject,{(OverBar[B])[ps_,(Subscript[M, K]^2),(Subscript[M, π]^2)]->(OverBar[B])[ps,(Subscript[M, π]^2),(Subscript[M, K]^2)],(OverBar[B])[ps_,(Subscript[M, η]^2),(Subscript[M, π]^2)]->(OverBar[B])[ps,(Subscript[M, π]^2),(Subscript[M, η]^2)],(OverBar[B])[ps_,(Subscript[M, η]^2),(Subscript[M, K]^2)]->(OverBar[B])[ps,(Subscript[M, K]^2),(Subscript[M, η]^2)]}],{(OverBar[B])[ps_,ms1_,ms2_]->tempfunc[(OverBar[B])[ps,ms1,ms2]]}],{tempfunc->Simplify}]

(* *)

Corder2[ject_]:=Do[proClist=ject;orderedClist={};switcher=0;Do[Do[If[MemberQ[orderedClist,Cswitch[proClist[[ii]],jj]],switcher=1],{jj,5}];If[switcher==0,orderedClist=Append[orderedClist,proClist[[ii]]],switcher=0],{ii,Length[proClist]}]]

ckorr[ject_]:=ReplaceAll[ReplaceAll[ject,{C[a_,b_,c_,d_,e_,f_]->tempfunc[C[a,b,(1/2)(a+b-c),d,e,f]]}],{tempfunc->Simplify}]

cflist=If[nonre==0,Delete[Union[Flatten[Table[If[momlist1[[ii]]⋂momlist1[[iii]]=={}&&ii<iii&&Length[momlist1[[ii]]]+Length[momlist1[[iii]]]<Length[imlist],C[sp[momlist2[[ii]],momlist2[[ii]]],sp[momlist2[[iii]],momlist2[[iii]]],sp[imp3=Sum[klim[[kkk]],{kkk,Length[klim=Complement[imlist,momlist1[[ii]],momlist1[[iii]]]]}],imp3],propmasslist[[jj]],propmasslist[[kk]],propmasslist[[ll]]]],{ii,1,Length[momlist2]},{iii,1,Length[momlist2]},{jj,Length[propmasslist]},{kk,Length[propmasslist]},{ll,Length[propmasslist]}]]],1],{}];

Corder2[cflist]

(Subscript[orderedClist, 1])=orderedClist;

Do[(Subscript[orderedClist, jj+1])=Table[Cswitch[orderedClist[[ii]],jj],{ii,Length[orderedClist]}],{jj,1,5}]

alist={(OverBar[A])[(Subscript[M, π]^2)],(OverBar[A])[(Subscript[M, K]^2)],(OverBar[A])[(Subscript[M, η]^2)]};

(OverBar[B])[0,ms1_,ms2_]:=0

blist=Complement[Union[Border[ReplaceAll[Simplify[replacements[transform[bflist]]],{B->(OverBar[B])}]]],{0}];

Remove[rtorderedClist]

Do[(Subscript[rtorderedClist, ii])=Simplify[replacements[transform[(Subscript[orderedClist, ii])]]],{ii,6}]

Do[Do[Do[If[(Subscript[rtorderedClist, jj])[[ii]]==(Subscript[rtorderedClist, kk])[[ii]],xft=jj;(Subscript[rtorderedClist, jj])=ReplacePart[(Subscript[rtorderedClist, jj]),equal,ii]],{kk,1,jj-1}],{jj,2,6}],{ii,1,Length[(Subscript[rtorderedClist, 1])]}]

(* *)

replacelist2[ject_,novel_,inpos_]:=ReplacePart[ject,novel,{{inpos[[1]]},{inpos[[2]]}},{{1},{2}}]

tgcorr[ject_]:=If[norecurs==1,ReplaceAll[argsym[ject],{TG1[x__]->TGpp[x,1],TG2[x__]->TGpp[x,2],TG3[x__]->TGpp[x,3],TG4[x__]->TGppp[x,4],TG5[x__]->TGppp[x,5],TG6[x__]->TGppp[x,6],TG7[x__]->TGppp[x,7]}]]

argsym[ject_]:=ReplaceAll[ject,{TG2[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]->TG2[x,y,u,u,w,ms1,ms2,ms3],TG3[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]->TG3[x,y,u,u,w,ms1,ms2,ms3],TG6[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]->TG6[x,y,u,u,w,ms1,ms2,ms3],TG7[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]->TG7[x,y,u,u,u,ms1,ms2,ms3]}]

Gswitch[ject_,posi1_,posi2_]:=Do[saveda=ject[[posi1]];savedb=ject[[posi2]];Return[replacelist2[ject,{savedb,saveda},{posi1,posi2}]]]

Gswitch2p[ject_,1]:=Gswitch[ject,3,4]

Gswitch2p[ject_,2]:=Gswitch[ject,6,7]

Gswitch3p[ject_]:=Gswitch[Gswitch[Gswitch[ject,1,2],3,4],6,8]

intg2pa[a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_,num_]:=If[{a6,a7}==Sort[{a6,a7}],intg2pb[a1,a2,a3,a4,a5,a6,a7,a8,num],Gswitch2p[intg2pb[a1,a2,a3,a4,a5,a6,a7,a8,num],2]*(((-1))^(num+1)),Gswitch2p[intg2pb[a1,a2,a3,a4,a5,a6,a7,a8,num],2]*(((-1))^(num+1))]

intg2pb[a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_,num_]:=If[{a3,a4}==Sort[{a3,a4}],TGpp[a1,a2,a3,a4,a5,a6,a7,a8,num],Gswitch2p[TGpp[a1,a2,a3,a4,a5,a6,a7,a8,num],1],Gswitch2p[TGpp[a1,a2,a3,a4,a5,a6,a7,a8,num],1]]

intg3p[a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_,num_]:=If[{a1,a2}==Sort[{a1,a2}],TGppp[a1,a2,a3,a4,a5,a6,a7,a8,num],Gswitch3p[TGppp[a1,a2,a3,a4,a5,a6,a7,a8,num]],Gswitch3p[TGppp[a1,a2,a3,a4,a5,a6,a7,a8,num]]]

Gminus[ject_]:=ReplaceAll[ject,{TGppp[x_,-y_,u_,v_,w_,ms1_,ms2_,ms3_,4]->-TGppp[-x,y,u,v,w,ms1,ms2,ms3,4],TGppp[x_,-y_,u_,v_,w_,ms1_,ms2_,ms3_,5]->TGppp[-x,y,u,v,w,ms1,ms2,ms3,5],TGppp[x_,-y_,u_,v_,w_,ms1_,ms2_,ms3_,6]->-TGppp[-x,y,u,v,w,ms1,ms2,ms3,6],TGppp[x_,-y_,u_,v_,w_,ms1_,ms2_,ms3_,7]->TGppp[-x,y,u,v,w,ms1,ms2,ms3,7],TGppp[x_,-y1_-y2_,u_,v_,w_,ms1_,ms2_,ms3_,4]->-TGppp[-x,y1+y2,u,v,w,ms1,ms2,ms3,4],TGppp[x_,-y1_-y2_,u_,v_,w_,ms1_,ms2_,ms3_,5]->TGppp[-x,y1+y2,u,v,w,ms1,ms2,ms3,5],TGppp[x_,-y1_-y2_,u_,v_,w_,ms1_,ms2_,ms3_,6]->-TGppp[-x,y1+y2,u,v,w,ms1,ms2,ms3,6],TGppp[x_,-y1_-y2_,u_,v_,w_,ms1_,ms2_,ms3_,7]->TGppp[-x,y1+y2,u,v,w,ms1,ms2,ms3,7]}]

Gorder[ject_]:=Gminus[ReplaceAll[tgcorr[ject],{TGpp[x__]->intg2pa[x],TGppp[x__]->intg3p[x]}]]

Gmake[ject_]:=ReplaceAll[ject,{TGpp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,1]->G1[x,ms1,ms2](Subscript[u, α])(Subscript[v, β]),TGpp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,2]->G2[x,ms1,ms2](Subscript[u, α]),TGpp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,3]->G3[x,ms1,ms2],TGppp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,4]->G4[x,y,ms1,ms2,ms3](Subscript[u, α])(Subscript[v, β])(Subscript[w, γ]),TGppp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,5]->G5[x,y,ms1,ms2,ms3](Subscript[u, α])(Subscript[v, β]),TGppp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,6]->G6[x,y,ms1,ms2,ms3](Subscript[u, α]),TGppp[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_,7]->G7[x,y,ms1,ms2,ms3]}]

Gind[ject_]:=ReplaceAll[ject,{G1->((Subscript[G, 1])^(αβ)),G2->((Subscript[G, 2])^α),G3->(Subscript[G, 3]),G4->((Subscript[G, 4])^(αβγ)),G5->((Subscript[G, 5])^(αβ)),G6->((Subscript[G, 6])^α),G7->(Subscript[G, 7])}]

(* *)

Glist=If[stcheck==0,{G1,G2,G3,G3,G4,G5,G6,G7},{G1,G2,G3,G3,G1,G2,G2,G2,G2,G2,G3,G3,G3,G3,G3,G3,G3,G3,G3,G3,G4,G5,G6,G7,G4,G5,G5,G5,G5,G5,G6,G6,G6,G6,G6,G7,G7,G7,G7}];

Gshortlist={G1,G2,G3,G4,G5,G6,G7};

Do[(Subscript[gexp, ii])=Gmake[Gorder[(Subscript[expression, ii])]];(Subscript[gexplist, ii])=Position[(Subscript[gexp, ii]),Glist[[ii]]];(Subscript[gexplist2, ii])=Table[Delete[(Subscript[gexplist, ii])[[jj]],-1],{jj,Length[(Subscript[gexplist, ii])]}];(Subscript[gexplist3, ii])=Union[Table[Extract[(Subscript[gexp, ii]),(Subscript[gexplist2, ii])[[jj]]],{jj,Length[(Subscript[gexplist2, ii])]}]],{ii,If[stcheck==0,8,39]}]

Do[(Subscript[simpgexp, ii])=If[stcheck==0,Sum[Simplify[Coefficient[(Subscript[gexp, ii]),(Subscript[gexplist3, ii])[[jj]]]]*(Subscript[gexplist3, ii])[[jj]],{jj,Length[(Subscript[gexplist3, ii])]}],Sum[Simplify[Coefficient[Coefficient[(Subscript[gexp, ii]),(Subscript[G, 8])],(Subscript[gexplist3, ii])[[jj]]]]*(Subscript[gexplist3, ii])[[jj]],{jj,Length[(Subscript[gexplist3, ii])]}]*(Subscript[G, 8])+Sum[Simplify[Coefficient[Coefficient[(Subscript[gexp, ii]),(Subscript[G, 27])],(Subscript[gexplist3, ii])[[jj]]]]*(Subscript[gexplist3, ii])[[jj]],{jj,Length[(Subscript[gexplist3, ii])]}]*(Subscript[G, 27])],{ii,If[stcheck==0,8,39]}]

sumexp=Sum[(Subscript[simpgexp, ii]),{ii,If[stcheck==0,8,39]}];

Do[(Subscript[gsumexplist, ii])=Position[sumexp,Gshortlist[[ii]]];(Subscript[gsumexplist2, ii])=Table[Delete[(Subscript[gsumexplist, ii])[[jj]],-1],{jj,Length[(Subscript[gsumexplist, ii])]}];(Subscript[gsumexplist3, ii])=Union[Table[Extract[sumexp,(Subscript[gsumexplist2, ii])[[jj]]],{jj,Length[(Subscript[gsumexplist2, ii])]}]],{ii,7}]

Do[(Subscript[fullexp, ii])=If[(Subscript[gsumexplist3, ii])≠{},Sum[Simplify[Coefficient[sumexp,(Subscript[gsumexplist3, ii])[[jj]]]]*(Subscript[gsumexplist3, ii])[[jj]],{jj,Length[(Subscript[gsumexplist3, ii])]}],0],{ii,1,7}]

endexptemp=Sum[Gind[(Subscript[fullexp, ii])],{ii,7}];

endexp=If[stcheck==0,endexptemp,If[nonre==1,ReplaceAll[getordns[ReplaceAll[endexptemp,{Dot[x_,y_]->sp[x,y]}],1,(Subscript[G, 8])],{sp[x_,y_]->Dot[x,y]}]+ReplaceAll[getordns[ReplaceAll[endexptemp,{Dot[x_,y_]->sp[x,y]}],1,(Subscript[G, 27])],{sp[x_,y_]->Dot[x,y]}],0]];

(* *)

(* In the following the loop integrals are defined. *)

TG1[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=a[sp[x,x],ms1,ms2]*sp[x,u]*sp[x,v]+b[sp[x,x],ms1,ms2]*sp[u,v]

TG2[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=I((1/2)*B[sp[x,x],ms1,ms2]-(Subscript[B, 11])[sp[x,x],ms1,ms2])*sp[x,u]

TG3[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=(1/4)*B[sp[x,x],ms1,ms2]

TG4[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=I(c[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[x,v]*sp[x,w]+c[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u]*sp[y,v]*sp[y,w]+d[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[x,v]*sp[y,w]+e[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[y,v]*sp[x,w]+f[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[y,u]*sp[x,v]*sp[x,w]+d[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u]*sp[y,v]*sp[x,w]+f[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u]*sp[x,v]*sp[y,w]+e[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[x,u]*sp[y,v]*sp[y,w]+g[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[v,w]+h[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,v]*sp[u,w]+i[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,w]*sp[u,v]+h[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u]*sp[v,w]+g[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,v]*sp[u,w]+i[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,w]*sp[u,v])

TG5[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=j[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[u,v]+k[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[x,v]+k[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u]*sp[y,v]+l[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]*sp[y,v]+m[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[y,u]*sp[x,v]

TG6[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=I(n[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]*sp[x,u]+n[sp[y,y],sp[x,x],sp[x,y],ms1,ms3,ms2]*sp[y,u])

TG7[x_,y_,u_,v_,w_,ms1_,ms2_,ms3_]:=(1/6)C[sp[x,x],sp[y,y],sp[x,y],ms1,ms2,ms3]

(* *)

a[ps_,ms1_,ms2_]:=-(1/4)B[ps,ms1,ms2]+(Subscript[B, 11])[ps,ms1,ms2]-(Subscript[B, 22])[ps,ms1,ms2]

b[ps_,ms1_,ms2_]:=-(Subscript[B, 20])[ps,ms1,ms2]

c[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/3)((Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]-4(Subscript[C, 22])[ps1,ps2,ps12,ms1,ms2,ms3]+4(Subscript[C, 33])[ps1,ps2,ps12,ms1,ms2,ms3])

d[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/3)((Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[C, 22])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[(OverHat[C]), 22])[ps1,ps2,ps12,ms1,ms2,ms3]+4(Subscript[C, 33])[ps1,ps2,ps12,ms1,ms2,ms3])

e[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/3)(-C[ps1,ps2,ps12,ms1,ms2,ms3]/2+2(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]+(Subscript[C, 11])[ps2,ps1,ps12,ms1,ms3,ms2]-2(Subscript[C, 22])[ps1,ps2,ps12,ms1,ms2,ms3]-4(Subscript[(OverHat[C]), 22])[ps1,ps2,ps12,ms1,ms2,ms3]+4(Subscript[C, 33])[ps1,ps2,ps12,ms1,ms2,ms3])

f[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-(2/3)((Subscript[(OverHat[C]), 22])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[(OverHat[C]), 33])[ps1,ps2,ps12,ms1,ms2,ms3])

g[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-(2/3)((Subscript[C, 20])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[C, 31])[ps1,ps2,ps12,ms1,ms2,ms3])

h[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(4/3)(Subscript[C, 31])[ps1,ps2,ps12,ms1,ms2,ms3]

i[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-(2/3)((Subscript[C, 20])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[C, 31])[ps1,ps2,ps12,ms1,ms2,ms3])

j[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-2(Subscript[C, 20])[ps1,ps2,ps12,ms1,ms2,ms3]

k[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]-2(Subscript[C, 22])[ps1,ps2,ps12,ms1,ms2,ms3]

l[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-C[ps1,ps2,ps12,ms1,ms2,ms3]/2+(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]+(Subscript[C, 11])[ps2,ps1,ps12,ms1,ms3,ms2]-2(Subscript[(OverHat[C]), 22])[ps1,ps2,ps12,ms1,ms2,ms3]

m[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=-2(Subscript[(OverHat[C]), 22])[ps1,ps2,ps12,ms1,ms2,ms3]

n[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=C[ps1,ps2,ps12,ms1,ms2,ms3]/2-(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]

(* *)

(Subscript[B, 11])=Function[{ps,ms1,ms2},If[replacements[transform[ps]]==0,If[ms1==ms2,(((A[ms1])/ms1)-(1/(16(π^2))))/2,,(((ms1^2)-(ms2^2))/(64 (((ms1-ms2))^2) (π^2)))+((ms1 A[ms1])/(2 (((ms1-ms2))^2)))+(((ms2-2 ms1) A[ms2])/(2 (((ms1-ms2))^2)))],,(1/(2ps))(-A[ms1]+A[ms2]+B[ps,ms1,ms2]*(ms1-ms2+ps))]];

(Subscript[B, 20])=Function[{ps,ms1,ms2},If[replacements[transform[ps]]==0,If[ms1==ms2,((A[ms1])/2),,(1/4)(A[ms2]+ms1 ((A[ms1]-A[ms2])/(ms1-ms2))+(1/(32(π^2)))(ms1+ms2))],,(1/6)(A[ms2]+2B[ps,ms1,ms2]*ms1+(3ms1+3ms2-ps)/(48(π^2))-(Subscript[B, 11])[ps,ms1,ms2]*(ms1-ms2+ps))] ];

(Subscript[B, 22])=Function[{ps,ms1,ms2},If[replacements[transform[ps]]==0,If[ms1==ms2,((A[ms])/(3 ms))-(1/(48 (π^2))),,((A[ms1] (ms1^2))/(3 (((ms1-ms2))^3)))-(((3 (ms1^2)-3 ms2 ms1+(ms2^2)) A[ms2])/(3 (((ms1-ms2))^3)))+((5 (ms1^3)-9 (ms2^2) ms1+4 (ms2^3))/(288 (((ms1-ms2))^3) (π^2)))],,(1/(3ps))(A[ms2]-B[ps,ms1,ms2]*ms1-(3ms1+3ms2-ps)/(96(π^2))+2(Subscript[B, 11])[ps,ms1,ms2]*(ms1-ms2+ps))]];

(Subscript[C, 11])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 11])[ps1,ps2,ps12,ms1,ms2,ms3]*ps2-(Subscript[H, 11])[ps2,ps1,ps12,ms1,ms3,ms2]*ps12)];

(Subscript[C, 20])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/4)(2C[ps1,ps2,ps12,ms1,ms2,ms3]*ms1+1/(16(π^2))+B[ps1-2ps12+ps2,ms2,ms3]-(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]*(ps1+ms1-ms2)-(Subscript[C, 11])[ps2,ps1,ps12,ms1,ms3,ms2]*(ps2+ms1-ms3))];

(Subscript[C, 22])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 21])[ps1,ps2,ps12,ms1,ms2,ms3]*ps2-(Subscript[H, 22])[ps2,ps1,ps12,ms1,ms3,ms2]*ps12)];

(Subscript[(OverHat[C]), 22])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 22])[ps2,ps1,ps12,ms1,ms3,ms2]*ps1-(Subscript[H, 21])[ps1,ps2,ps12,ms1,ms2,ms3]*ps12)];

(Subscript[C, 33])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 31])[ps1,ps2,ps12,ms1,ms2,ms3]*ps2-(Subscript[H, 32])[ps2,ps1,ps12,ms1,ms3,ms2]*ps12)];

(Subscript[(OverHat[C]), 33])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 32])[ps2,ps1,ps12,ms1,ms3,ms2]*ps1-(Subscript[H, 31])[ps1,ps2,ps12,ms1,ms2,ms3]*ps12)];

(Subscript[C, 31])=Function[{ps1,ps2,ps12,ms1,ms2,ms3},(1/(ps1*ps2-(ps12^2)))((Subscript[H, 30])[ps1,ps2,ps12,ms1,ms2,ms3]*ps2-(Subscript[H, 30])[ps2,ps1,ps12,ms1,ms3,ms2]*ps12)];

(* *)

(Subscript[H, 11])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/2)(-B[ps2,ms1,ms3]+B[ps1-2ps12+ps2,ms2,ms3]+C[ps1,ps2,ps12,ms1,ms2,ms3]*(ps1 + ms1 - ms2))

(Subscript[H, 21])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/2)(B[ps1-2ps12+ps2,ms2,ms3]+(Subscript[C, 11])[ps1,ps2,ps12,ms1,ms2,ms3]*(ps1 + ms1 - ms2)-(Subscript[B, 11])[ps1-2ps12+ps2,ms2,ms3])-(Subscript[C, 20])[ps1,ps2,ps12,ms1,ms2,ms3]

(Subscript[H, 22])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:= (1/2)((Subscript[C, 11])[ps2,ps1,ps12,ms1,ms3,ms2]*(ps1 + ms1 - ms2)-(Subscript[B, 11])[ps2,ms1,ms3]+(Subscript[B, 11])[ps1-2ps12+ps2,ms2,ms3])

(Subscript[H, 30])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/2)((Subscript[C, 20])[ps1,ps2,ps12,ms1,ms3,ms2]*(ps1 + ms1 - ms2)-(Subscript[B, 20])[ps2,ms1,ms3]+(Subscript[B, 20])[ps1-2ps12+ps2,ms2,ms3])

(Subscript[H, 31])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/2)(B[ps1-2ps12+ps2,ms2,ms3]+(Subscript[C, 22])[ps1,ps2,ps12,ms1,ms2,ms3]*(ps1 + ms1 - ms2)-2*(Subscript[B, 11])[ps1-2ps12+ps2,ms2,ms3]+(Subscript[B, 22])[ps1-2ps12+ps2,ms2,ms3])-2(Subscript[C, 31])[ps1,ps2,ps12,ms1,ms2,ms3]

(Subscript[H, 32])[ps1_,ps2_,ps12_,ms1_,ms2_,ms3_]:=(1/2)((Subscript[C, 22])[ps2,ps1,ps12,ms1,ms3,ms2]*(ps1 + ms1 - ms2)-(Subscript[B, 22])[ps2,ms1,ms3]+(Subscript[B, 22])[ps1-2ps12+ps2,ms2,ms3])

B[ps_,ms1_,ms2_]:=(OverBar[B])[ps,ms1,ms2]+If[ms1≠ms2,((A[ms1]-A[ms2])/(ms1-ms2)),((A[ms1])/ms1)-(1/(16(π^2))),((A[ms1]-A[ms2])/(ms1-ms2))]

A[ms_]:=(OverBar[A])[ms]

rtgexpression1=If[nonre==0,Sum[Cchange[Border[replacements[transform[(Subscript[expression, ii])]]]],{ii,If[stcheck==0,4,20]}],0];

loopresult1a=Sum[simp2[Coefficient[rtgexpression1+renomp4+tadpoles,alist[[ii]]]]*alist[[ii]],{ii,Length[alist]}];

loopresult1b=Sum[Simplify[Coefficient[rtgexpression1,blist[[ii]]]]*blist[[ii]],{ii,Length[blist]}];

loopresult1poly=simp2[ReplaceAll[rtgexpression1+renomp4+tadpoles,{(OverBar[A])[__]->0,(OverBar[B])[__]->0,C[__]->0,(Subscript[L, i_])->0}]];

rtgexpression2=If[nonre==0,Sum[Cchange[Border[replacements[transform[If[ii≠25,(Subscript[expression, ii]),0]]]]],{ii,If[stcheck==0,6,22],If[stcheck==0,8,39]}],0];

loopresult2a=Sum[Simplify[Coefficient[rtgexpression2,alist[[ii]]]]*alist[[ii]],{ii,Length[alist]}];

loopresult2b=Sum[Simplify[Coefficient[rtgexpression2,blist[[ii]]]]*blist[[ii]],{ii,Length[blist]}];

loopresult2poly=Simplify[ReplaceAll[rtgexpression2,{(OverBar[A])[__]->0,(OverBar[B])[__]->0,C[__]->0}]];

loopresult2c=Sum[Simplify[Sum[Coefficient[rtgexpression2,(Subscript[rtorderedClist, jj])[[ii]]],{jj,6}]]*(Subscript[rtorderedClist, 1])[[ii]],{ii,Length[(Subscript[rtorderedClist, 1])]}];

Do[Do[resTG4a[ii]=0,{ii,Length[alist]}];If[particleorder>4,Do[tempex=Cchange[Border[replacements[transform[ReplaceAll[If[stcheck==0,(Subscript[expressiont, 5])[[jj]],(Subscript[expressiont, 21])[[jj]]+(Subscript[expressiont, 25])[[jj]]],{TGG4->TG4}]]]]];Do[resTG4a[ii]=Simplify[Coefficient[resTG4a[ii]+tempex,alist[[ii]]]]*alist[[ii]],{ii,Length[alist]}],{jj,Length[(Subscript[expressiont, 5])]}]]]

loopresult3a=Sum[resTG4a[ii],{ii,Length[alist]}];

Do[Do[resTG4b[ii]=0,{ii,Length[blist]}];If[particleorder>4,Do[tempex=Cchange[Border[replacements[transform[ReplaceAll[If[stcheck==0,(Subscript[expressiont, 5])[[jj]],(Subscript[expressiont, 21])[[jj]]+(Subscript[expressiont, 25])[[jj]]],{TGG4->TG4}]]]]];Do[resTG4b[ii]=Simplify[Coefficient[resTG4b[ii]+tempex,blist[[ii]]]]*blist[[ii]],{ii,Length[blist]}],{jj,Length[(Subscript[expressiont, 5])]}]]]

loopresult3b=Sum[resTG4b[ii],{ii,Length[blist]}];

Do[Do[resTG4c[ii]=0,{ii,Length[(Subscript[rtorderedClist, 1])]}];If[particleorder>4,Do[tempex=Cchange[Border[replacements[transform[ReplaceAll[If[stcheck==0,(Subscript[expressiont, 5])[[jj]],(Subscript[expressiont, 21])[[jj]]+(Subscript[expressiont, 25])[[jj]]],{TGG4->TG4}]]]]];Do[resTG4c[ii]=Simplify[Sum[Coefficient[resTG4c[ii]+tempex,(Subscript[rtorderedClist, kk])[[ii]]],{kk,6}]]*(Subscript[rtorderedClist, 1])[[ii]],{ii,Length[(Subscript[rtorderedClist, 1])]}],{jj,Length[(Subscript[expressiont, 5])]}]]]

loopresult3c=Sum[resTG4c[ii],{ii,Length[(Subscript[rtorderedClist, 1])]}];

Do[resTG4poly=0;If[particleorder>4,Do[tempex=replacements[transform[ReplaceAll[ReplaceAll[If[stcheck==0,(Subscript[expressiont, 5])[[jj]],(Subscript[expressiont, 21])[[jj]]+(Subscript[expressiont, 25])[[jj]]],{TGG4->TG4}],{(OverBar[A])[__]->0,(OverBar[B])[__]->0,C[__]->0}]]];resTG4poly=Simplify[resTG4poly+Expand[tempex]],{jj,Length[(Subscript[expressiont, 5])]}]]]

loopresult3poly=resTG4poly;

loopresult=ReplaceAll[Sum[simp2[etakorr[Coefficient[loopresult1a+loopresult2a+loopresult3a,alist[[ii]]]]]*alist[[ii]],{ii,Length[alist]}]+Sum[Simplify[etakorr[Coefficient[loopresult1b+loopresult2b+loopresult3b,blist[[ii]]]]]*blist[[ii]],{ii,Length[blist]}]+Sum[Simplify[etakorr[Coefficient[loopresult2c+loopresult3c,(Subscript[rtorderedClist, 1])[[ii]]]]]*ckorr[(Subscript[rtorderedClist, 1])[[ii]]],{ii,Length[(Subscript[rtorderedClist, 1])]}]+simp2[etakorr[loopresult1poly+loopresult2poly+loopresult3poly]]+endexp,{F->(Subscript[F, π])}];

result=normalization[treeresult2temp+treeresult4temp+loopresult];

title=If[onlyloops==1,"Loop amplitude:","Complete amplitude to 0(p^4):","Complete amplitude to 0(p^4):"];

output[result]

NotebookClose[notesub];