cross(a,b) = (
dot(a[2],b[3]) - dot(a[3],b[2]),
dot(a[3],b[1]) - dot(a[1],b[3]),
dot(a[1],b[2]) - dot(a[2],b[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 = (gamma1,gamma2,gamma3)

alpha1 = dot(gamma0,gamma1)
alpha2 = dot(gamma0,gamma2)
alpha3 = dot(gamma0,gamma3)

alpha = (alpha1,alpha2,alpha3)

beta = gamma0

phat(f) = -i (d(f,x), d(f,y), d(f,z))

Psi = psi(x,y,z)

H(f) = test(rank(f) == 0, H1(f), H2(f))

-- f is a scalar function

H1(f) = -i (dot(alpha1,d(f,x)) + dot(alpha2,d(f,y)) + dot(alpha3,d(f,z)) + m dot(beta,f))

-- f is a vector

H2(f) = (H1(f[1]),H1(f[2]),H1(f[3]))

"Verify commutator [H,L]"

L(f) = cross((x,y,z),phat(f))

check(H(L(Psi)) - L(H(Psi)) == i cross(phat(Psi),alpha))

"ok"

"Verify commutator [H,S]"

check(cross(gamma,gamma) == -cross(alpha,alpha))

S = -i/4 cross(alpha,alpha)

HS = H(S Psi)

SH = (dot(S[1],H(Psi)), dot(S[2],H(Psi)), dot(S[3],H(Psi)))

check(HS - SH == -i cross(phat(Psi),alpha))

"ok"