diff --git a/feynamp/form/form.py b/feynamp/form/form.py
index ffd1c6d..c5fa804 100644
--- a/feynamp/form/form.py
+++ b/feynamp/form/form.py
@@ -17,7 +17,7 @@
 Tensors colorcorrelation;
 Tensors Metric(symmetric),df(symmetric),da(symmetric),Identity(symmetric);
 Function ProjM,ProjP,VF,xg,xgi,P,dg,dgi,xeg,xegi;
-CFunctions Den,Denom,P,Gamma,u,v,ubar,vbar,eps,epsstar,VC,VA,GammaId, GammaCollect, GammaIdCollect;
+CFunctions Den,Denom,P,Gamma,u,v,ubar,vbar,eps,epsstar,VC,VA,VPol,GammaId, GammaCollect, GammaIdCollect;
 Indices a,o,n,m,tm,tn,beta,b,m,betap,alphap,a,alpha,ind,delta,k,j,l,c,d,e;
 """
 
diff --git a/feynamp/form/lorentz.py b/feynamp/form/lorentz.py
index 4e7be6c..8d0a2f4 100644
--- a/feynamp/form/lorentz.py
+++ b/feynamp/form/lorentz.py
@@ -138,7 +138,11 @@ def get_orthogonal_polarisation_momentum(
 
 
 def get_polarisation_sums(fds: List[FeynmanDiagram], model: FeynModel):
-    pol_sums = ""
+    pol_sums = """
+    repeat;
+    id VPol(Polb?,Moma?) * VPol(Pold?,Moma?) = d_(Polb,Pold);
+    endrepeat;
+    """
     # TODO might want to loop over all fds?
     for fd in [fds[0]]:
         for l in fd.legs:
@@ -162,21 +166,21 @@ def get_polarisation_sums(fds: List[FeynmanDiagram], model: FeynModel):
 
 def get_polarisation_sum_massive(mom_a):
     pol_sum = f"""
-    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Pold?,{mom_a}) = -Metric(Mul,Muc) + (P(Mul,{mom_a})*P(Muc,{mom_a}))*Den({mom_a}.{mom_a});
+    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Polb?,{mom_a}) = -Metric(Mul,Muc) + (P(Mul,{mom_a})*P(Muc,{mom_a}))*Den({mom_a}.{mom_a});
     """
     return pol_sum
 
 
 def get_polarisation_sum_feynman(mom_a):
     pol_sum = f"""
-    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Pold?,{mom_a}) = -Metric(Mul,Muc);
+    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Polb?,{mom_a}) = -Metric(Mul,Muc);
     """
     return pol_sum
 
 
 def get_polarisation_sum_physical(mom_a, mom_b):
     pol_sum = f"""
-    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Pold?,{mom_a}) = -Metric(Muc,Mul) 
+    id epsstar(Muc?,Polb?,{mom_a}) * eps(Mul?,Polb?,{mom_a}) = -Metric(Muc,Mul) 
     + (P(Muc,{mom_a})*P(Mul,{mom_b}) +  P(Mul,{mom_a})*P(Muc,{mom_b}))*Den({mom_b}.{mom_a}) 
     - P(Muc,{mom_a})*P(Mul,{mom_a})*({mom_b}.{mom_b})*Den({mom_b}.{mom_a})*Den({mom_b}.{mom_a});
     """
diff --git a/feynamp/leg.py b/feynamp/leg.py
index 1194b97..1b9a98b 100644
--- a/feynamp/leg.py
+++ b/feynamp/leg.py
@@ -112,9 +112,9 @@ def get_leg_math(fd, leg, model):  # epsilons or u/v optionally also barred
 
     if p.spin == 3:
         if leg.is_incoming():
-            ret += f"eps(Mu{p.particle.id},Pol{p.particle.id},{mom})"
+            ret += f"eps(Mu{p.particle.id},Pol{p.particle.id},{mom})*VPol(Pol{p.particle.id},{mom})"
         else:
-            ret += f"eps_star(Mu{p.particle.id},Pol{p.particle.id},{mom})"
+            ret += f"eps_star(Mu{p.particle.id},Pol{p.particle.id},{mom})*VPol(Pol{p.particle.id},{mom})"
     if p.spin == 2:
         if not p.particle.is_anti():
             if leg.is_incoming():