"Verify Gordon decomposition"
E1 = sqrt(p1x^2 + p1y^2 + p1z^2 + m^2)
E2 = sqrt(p2x^2 + p2y^2 + p2z^2 + m^2)
p1 = (E1, p1x, p1y, p1z)
p2 = (E2, p2x, p2y, p2z)
u11 = (E1 + m, 0, p1z, p1x + i p1y)
u12 = (0, E1 + m, p1x - i p1y, -p1z)
u21 = (E2 + m, 0, p2z, p2x + i p2y)
u22 = (0, E2 + m, p2x - i p2y, -p2z)
u1 = (u11,u12)
u2 = (u21,u22)
I = ((1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1))
gmunu = ((1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,-1))
gamma0 = ((1,0,0,0),(0,1,0,0),(0,0,-1,0),(0,0,0,-1))
gamma1 = ((0,0,0,1),(0,0,1,0),(0,-1,0,0),(-1,0,0,0))
gamma2 = ((0,0,0,-i),(0,0,i,0),(0,i,0,0),(-i,0,0,0))
gamma3 = ((0,0,1,0),(0,0,0,-1),(-1,0,0,0),(0,1,0,0))
gamma = (gamma0,gamma1,gamma2,gamma3)
u2bar = dot(conj(u2),gamma0) -- adjoint of u2
T = transpose(dot(gamma,transpose(gamma)),2,3)
sigmamunu = i/2 (T - transpose(T))
S = transpose(sigmamunu,2,3,3,4)
J = (outer(p2+p1,I) + i dot(S,gmunu,p2-p1)) / (2 m)
for(s1,1,2, -- for each spin state
for(s2,1,2,
for(mu,1,4, -- for each gamma
 check( dot(u2bar[s2],gamma[mu],u1[s1]) == dot(u2bar[s2],J[mu],u1[s1]) )
)))
"ok"