gps[[Mu]_,[Nu]_,[Rho]_,[Sigma]_]:=(MTD[[Mu],[Rho]]-(FVD[k,[Mu]]FVD[k,[Rho]])/SPD[k,k])(MTD[[Nu],[Sigma]]-(FVD[k,[Nu]]FVD[k,[Sigma]])/SPD[k,k])+(MTD[[Mu],[Sigma]]-(FVD[k,[Mu]]FVD[k,[Sigma]])/SPD[k,k])(MTD[[Nu],[Rho]]-(FVD[k,[Nu]]FVD[k,[Rho]])/SPD[k,k])-2/3 (MTD[[Mu],[Nu]]-(FVD[k,[Mu]]FVD[k,[Nu]])/SPD[k,k])(MTD[[Rho],[Sigma]]-(FVD[k,[Rho]]FVD[k,[Sigma]])/SPD[k,k])
gp[1]=I*(GSD[p3]+m).GAD[[Gamma]].(GSD[p4]-m)
gpc[1]=-I*GAD[[Epsilon]];
gammapart[a_, b_] := TR[gp[a].gpc[b]]
Y[k_,[Mu]_,[Nu]_,[Alpha]_,[Beta]_,[Gamma]_]:=FVD[k,[Mu]](MTD[[Nu],[Beta]]MTD[[Alpha],[Gamma]]-MTD[[Nu],[Gamma]]MTD[[Alpha],[Beta]])+FVD[k,[Beta]](MTD[[Mu],[Alpha]]MTD[[Nu],[Gamma]]-1/2 MTD[[Mu],[Nu]]MTD[[Alpha],[Gamma]])-FVD[k,[Gamma]](MTD[[Mu],[Alpha]]MTD[[Nu],[Beta]]-1/2 MTD[[Mu],[Nu]]MTD[[Alpha],[Beta]])
Z[k1_,k2_,k3_,[Mu]_,[Nu]_,[Alpha]_,[Beta]_,[Gamma]_]:=Y[k1,[Mu],[Nu],[Alpha],[Beta],[Gamma]]+Y[k2,[Mu],[Nu],[Beta],[Gamma],[Alpha]]+Y[k3,[Nu],[Mu],[Gamma],[Alpha],[Beta]]+Y[k1,[Nu],[Mu],[Alpha],[Beta],[Gamma]]+Y[k2,[Nu],[Mu],[Beta],[Gamma],[Alpha]]+Y[k3,[Nu],[Mu],[Gamma],[Alpha],[Beta]]
rule1 = {SPD[a_ + b_, a_ + b_] -> Calc[SPD[a + b, a + b]]};
rule2 = {SPD[a_ + b_ + c_, a_ + b_ + c_] ->     Calc[SPD[a + b + c, a + b + c]]};
cp[1]=Contract[Z[-k1,-k2,k1+k2-k,[Mu],[Nu],[Alpha],[Beta],[Gamma]]*gps[[Mu],[Nu],[Rho],[Sigma]]] 1/SPD[k1+k2-k,k1+k2-k]/.rule2;
cpc[1]=Z[-k1,-k2,k1+k2-k,[Rho],[Sigma],[Alpha],[Beta],[Epsilon]](1/SPD[k1+k2-k,k1+k2-k]/.rule2);
momentumcombination = {SPD[k1, k2] -> s12/2,   SPD[k1, p3] -> (m^2 - s13)/2,   SPD[k1, p4] -> (m^2 - s14)/2,   SPD[k1, k] -> (s12 + s13 + s14 - 2 m^2)/2,   SPD[k2, p3] -> (m^2 - s23)/2,   SPD[k2, p4] -> (m^2 - s24)/2,   SPD[k2, k] -> (s12 + s23 + s24 - 2 m^2)/2,   SPD[p3, p4] -> (-s12 - 2 s13 - 2 s14 - s24 + 2 m^2 + gm^2)/2,   SPD[p3, k] -> (s12 + s14 + s24 - 2 m^2 - gm^2)/2,   SPD[p4, k] -> (s12 + s13 + s23 - 2 m^2 - gm^2)/2,   SPD[k1, k1] -> 0,   SPD[k2, k2] -> 0,   SPD[p3, p3] -> m^2,   SPD[p4, p4] -> m^2,   SPD[k, k] -> gm^2};
cqpart[i_, j_] := (cp[i] /. FCI[momentumcombination]) (cpc[j] /.FCI[momentumcombination])
a11 = SUNSimplify[SUNF[a, b, c]SUNF[a, b, g] SUNTrace[SUNT[c, g], Explicit -> True], SUNNToCACF -> False]
bingo[i_, j_] :=   Contract[cqpart[i, j] gammapart[i, j]] /. FCI[momentumcombination]
SetDirectory[NotebookDirectory[]]mango[1] = bingo[1, 1] + Sum[{2*bingo[1, i]}, {i, 2, 7}][[1]]
Export["mango1.m", mango[1]]bingos1 = Simplify[1/64*a11*mango[1]]
Export["bingos1.m", FortranForm[bingos1 /. D -> 4]]