forked from sympy/sympy
-
Notifications
You must be signed in to change notification settings - Fork 1
/
qft.py
executable file
·140 lines (104 loc) · 3.32 KB
/
qft.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/usr/bin/env python
"""Quantum field theory example
* http://en.wikipedia.org/wiki/Quantum_field_theory
This particular example is a work in progress. Currently it calculates the
scattering amplitude of the process:
electron + positron -> photon -> electron + positron
in QED (http://en.wikipedia.org/wiki/Quantum_electrodynamics). The aim
is to be able to do any kind of calculations in QED or standard model in
SymPy, but that's a long journey.
"""
from sympy import Basic, exp, Symbol, sin, Rational, I, Mul, Matrix, \
ones, sqrt, pprint, simplify, Eq, sympify
from sympy.physics import msigma, mgamma
#gamma^mu
gamma0 = mgamma(0)
gamma1 = mgamma(1)
gamma2 = mgamma(2)
gamma3 = mgamma(3)
gamma5 = mgamma(5)
#sigma_i
sigma1 = msigma(1)
sigma2 = msigma(2)
sigma3 = msigma(3)
E = Symbol("E", real=True)
m = Symbol("m", real=True)
def u(p, r):
""" p = (p1, p2, p3); r = 0,1 """
assert r in [1, 2]
p1, p2, p3 = p
if r == 1:
ksi = Matrix([[1], [0]])
else:
ksi = Matrix([[0], [1]])
a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E + m) * ksi
if a == 0:
a = zeros(2, 1)
return sqrt(E + m) * Matrix([[ksi[0, 0]], [ksi[1, 0]], [a[0, 0]], [a[1, 0]]])
def v(p, r):
""" p = (p1, p2, p3); r = 0,1 """
assert r in [1, 2]
p1, p2, p3 = p
if r == 1:
ksi = Matrix([[1], [0]])
else:
ksi = -Matrix([[0], [1]])
a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E + m) * ksi
if a == 0:
a = zeros(2, 1)
return sqrt(E + m) * Matrix([[a[0, 0]], [a[1, 0]], [ksi[0, 0]], [ksi[1, 0]]])
def pslash(p):
p1, p2, p3 = p
p0 = sqrt(m**2 + p1**2 + p2**2 + p3**2)
return gamma0*p0 - gamma1*p1 - gamma2*p2 - gamma3*p3
def Tr(M):
return M.trace()
def xprint(lhs, rhs):
pprint( Eq(sympify(lhs), rhs ) )
def main():
a = Symbol("a", real=True)
b = Symbol("b", real=True)
c = Symbol("c", real=True)
p = (a, b, c)
assert u(p, 1).D * u(p, 2) == Matrix(1, 1, [0])
assert u(p, 2).D * u(p, 1) == Matrix(1, 1, [0])
p1, p2, p3 = [Symbol(x, real=True) for x in ["p1", "p2", "p3"]]
pp1, pp2, pp3 = [Symbol(x, real=True) for x in ["pp1", "pp2", "pp3"]]
k1, k2, k3 = [Symbol(x, real=True) for x in ["k1", "k2", "k3"]]
kp1, kp2, kp3 = [Symbol(x, real=True) for x in ["kp1", "kp2", "kp3"]]
p = (p1, p2, p3)
pp = (pp1, pp2, pp3)
k = (k1, k2, k3)
kp = (kp1, kp2, kp3)
mu = Symbol("mu")
e = (pslash(p) + m*ones(4))*(pslash(k) - m*ones(4))
f = pslash(p) + m*ones(4)
g = pslash(p) - m*ones(4)
#pprint(e)
xprint( 'Tr(f*g)', Tr(f*g) )
#print Tr(pslash(p) * pslash(k)).expand()
M0 = [ ( v(pp, 1).D * mgamma(mu) * u(p, 1) ) * ( u(k, 1).D * mgamma(mu, True) *
v(kp, 1) ) for mu in range(4)]
M = M0[0] + M0[1] + M0[2] + M0[3]
M = M[0]
assert isinstance(M, Basic)
#print M
#print simplify(M)
d = Symbol("d", real=True) # d=E+m
xprint('M', M)
print "-"*40
M = ((M.subs(E, d - m)).expand() * d**2 ).expand()
xprint('M2', 1/(E + m)**2 * M)
print "-"*40
x, y = M.as_real_imag()
xprint('Re(M)', x)
xprint('Im(M)', y)
e = x**2 + y**2
xprint('abs(M)**2', e)
print "-"*40
xprint('Expand(abs(M)**2)', e.expand())
#print Pauli(1)*Pauli(1)
#print Pauli(1)**2
#print Pauli(1)*2*Pauli(1)
if __name__ == "__main__":
main()