examples/qft.py

   1 #!/usr/bin/env python
   2 import iam_sympy_example
   3
   4 from sympy import Basic,exp,Symbol,sin,Rational,I,Mul, Matrix, \
   5     ones, sqrt, pprint, simplify, trim, Eq, sympify
   6
   7 from sympy.physics import msigma, mgamma
   8
   9 #gamma^mu
  10 gamma0=mgamma(0)
  11 gamma1=mgamma(1)
  12 gamma2=mgamma(2)
  13 gamma3=mgamma(3)
  14 gamma5=mgamma(5)
  15
  16 #sigma_i
  17 sigma1=msigma(1)
  18 sigma2=msigma(2)
  19 sigma3=msigma(3)
  20
  21 E = Symbol("E", real=True)
  22 m = Symbol("m", real=True)
  23
  24 def u(p,r):
  25     """ p = (p1, p2, p3); r = 0,1 """
  26     assert r in [1,2]
  27     p1,p2,p3 = p
  28     if r == 1:
  29         ksi = Matrix([ [1],[0] ])
  30     else:
  31         ksi = Matrix([ [0],[1] ])
  32     a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E+m) * ksi
  33     if a ==0:
  34         a = zeros((2, 1))
  35     return sqrt(E+m) * Matrix([ [ksi[0,0]], [ksi[1,0]], [a[0,0]], [a[1,0]] ])
  36
  37 def v(p,r):
  38     """ p = (p1, p2, p3); r = 0,1 """
  39     assert r in [1,2]
  40     p1,p2,p3 = p
  41     if r == 1:
  42         ksi = Matrix([ [1],[0] ])
  43     else:
  44         ksi = -Matrix([ [0],[1] ])
  45     a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E+m) * ksi
  46     if a ==0:
  47         a = zeros((2,1))
  48     return sqrt(E+m) * Matrix([ [a[0,0]], [a[1,0]], [ksi[0,0]], [ksi[1,0]] ])
  49
  50 def pslash(p):
  51     p1,p2,p3 = p
  52     p0 = sqrt(m**2+p1**2+p2**2+p3**2)
  53     return gamma0*p0-gamma1*p1-gamma2*p2-gamma3*p3
  54
  55 def Tr(M):
  56     assert M.lines == M.cols
  57     t = 0
  58     for i in range(M.lines):
  59         t+=M[i,i]
  60     return t
  61
  62 a=Symbol("a", real=True)
  63 b=Symbol("b", real=True)
  64 c=Symbol("c", real=True)
  65
  66 p = (a,b,c)
  67
  68 assert u(p, 1).D * u(p, 2) == Matrix(1, 1, [0])
  69 assert u(p, 2).D * u(p, 1) == Matrix(1, 1, [0])
  70
  71 p1,p2,p3 =[Symbol(x, real=True) for x in ["p1","p2","p3"]]
  72 pp1,pp2,pp3 =[Symbol(x, real=True) for x in ["pp1","pp2","pp3"]]
  73 k1,k2,k3 =[Symbol(x, real=True) for x in ["k1","k2","k3"]]
  74 kp1,kp2,kp3 =[Symbol(x, real=True) for x in ["kp1","kp2","kp3"]]
  75
  76 p = (p1,p2,p3)
  77 pp = (pp1,pp2,pp3)
  78
  79 k = (k1,k2,k3)
  80 kp = (kp1,kp2,kp3)
  81
  82 mu = Symbol("mu")
  83
  84 e = (pslash(p)+m*ones(4))*(pslash(k)-m*ones(4))
  85 f = pslash(p)+m*ones(4)
  86 g = pslash(p)-m*ones(4)
  87
  88
  89 def xprint(lhs, rhs):
  90     pprint( Eq(sympify(lhs), rhs ) )
  91
  92 #pprint(e)
  93 xprint( 'Tr(f*g)', Tr(f*g) )
  94 #print Tr(pslash(p) * pslash(k)).expand()
  95
  96 M0 = [ ( v(pp, 1).D * mgamma(mu) * u(p, 1) ) * ( u(k, 1).D * mgamma(mu,True) * \
  97         v(kp, 1) ) for mu in range(4)]
  98 M = M0[0]+M0[1]+M0[2]+M0[3]
  99 M = M[0]
 100 assert isinstance(M, Basic)
 101 #print M
 102 #print trim(M)
 103
 104 d=Symbol("d", real=True) #d=E+m
 105
 106 xprint('M', M)
 107 print "-"*40
 108 M = ((M.subs(E,d-m)).expand() * d**2 ).expand()
 109 xprint('M2', 1/(E+m)**2 * M)
 110 print "-"*40
 111 x,y= M.as_real_imag()
 112 xprint('Re(M)', x)
 113 xprint('Im(M)', y)
 114 e = x**2+y**2
 115 xprint('abs(M)**2', e)
 116 print "-"*40
 117 xprint('Expand(abs(M)**2)', e.expand())
 118
 119 #print Pauli(1)*Pauli(1)
 120 #print Pauli(1)**2
 121 #print Pauli(1)*2*Pauli(1)