"Dirac equation 1"

P = (E / c, px, py, pz)
X = (c t, x, y, z)

gmunu = ((1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,-1))

psi1 = exp(-i dot(P,gmunu,X) / hbar) w1
psi2 = exp(-i dot(P,gmunu,X) / hbar) w2

psi3 = exp(i dot(P,gmunu,X) / hbar) w3
psi4 = exp(i dot(P,gmunu,X) / hbar) w4

w1 = sqrt((E + m c^2) / (2 m c^2)) *
     (1, 0, pz c / (E + m c^2), (px + i py) c / (E + m c^2))

w2 = sqrt((E + m c^2) / (2 m c^2)) *
     (0, 1, (px - i py) c / (E + m c^2), -pz c / (E + m c^2))

w3 = sqrt((E + m c^2) / (2 m c^2)) *
     (pz c / (E + m c^2), (px + i py) c / (E + m c^2), 1, 0)

w4 = sqrt((E + m c^2) / (2 m c^2)) *
     ((px - i py) c / (E + m c^2), -pz c / (E + m c^2), 0, 1)

E = sqrt(px^2 c^2 + py^2 c^2 + pz^2 c^2 + m^2 c^4)

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))

D(psi) = 1/c dot(gamma0,d(psi,t)) +
         dot(gamma1,d(psi,x)) +
         dot(gamma2,d(psi,y)) +
         dot(gamma3,d(psi,z))

"Verify wavefunctions"

check(i hbar D(psi1) == m c psi1)
check(i hbar D(psi2) == m c psi2)
check(i hbar D(psi3) == m c psi3)
check(i hbar D(psi4) == m c psi4)

"ok"

"Verify normalization"

check(dot(psi1,conj(psi1)) == E / (m c^2))
check(dot(psi2,conj(psi2)) == E / (m c^2))
check(dot(psi3,conj(psi3)) == E / (m c^2))
check(dot(psi4,conj(psi4)) == E / (m c^2))

"ok"