読者です 読者をやめる 読者になる 読者になる

fortran66のブログ

fortran について書きます。

Python で有理数 6j symbol

メモ帳 とりあえず定義式通りに。計算量が少ない変形もあるんですが・・・

import sympy
import numpy

def fa(n):
    if n < 0:
        return 0
    else:
        return math.factorial(n)

def dicsort(d):
    return sorted(d.items(), key=lambda x: x[0])

def pr(i):
    a = sympy.factorint((i * i).p)
    b = sympy.factorint((i * i).q)
    return a, b #dicsort(a), dicsort(b)

def triangle(a, b, c):
    if a + b < c: return 0
    if b + c < a: return 0
    if c + a < b: return 0 
    a = sympy.Rational(a)
    b = sympy.Rational(b)
    c = sympy.Rational(c)
    s1 =  a + b - c
    s2 =  a - b + c
    s3 = -a + b + c
    s4 =  a + b + c + 1
    t = sympy.Rational( fa(s1) * fa(s2) * fa(s3), fa(s4) )
    return sympy.sqrt(t)

def sixJ(j1, j2, j3, l1, l2, l3):
    j1 = sympy.Rational(j1)
    j2 = sympy.Rational(j2)
    j3 = sympy.Rational(j3)
    l1 = sympy.Rational(l1)
    l2 = sympy.Rational(l2)
    l3 = sympy.Rational(l3)
    factor = triangle(j1, j2, j3) * triangle(j1, l2, l3) * triangle(l1, j2, l3) * triangle(l1, l2, j3)
    J = j1 + j2 + j3
    jl1 = j1 + l2 + l3
    jl2 = l1 + j2 + l3
    jl3 = l1 + l2 + j3
    jl12 = j1 + j2 + l1 + l2
    jl23 = j2 + j3 + l2 + l3
    jl31 = j3 + j1 + l3 + l1
    kmin = max(J, jl1, jl2, jl3)
    kmax = min(jl12, jl23, jl31)
    s = 0
    for k in arange(kmin, kmax + 1):
        s = s + sympy.Rational( (-1)**k * fa(k + 1),
                fa(k - J) * fa(k - jl1) * fa(k - jl2) * fa(k - jl3) *
                fa(jl12 -k) * fa(jl23 - k) * fa(jl31 - k) )
    return factor * s

def chk6(j1, j2, j3, l1, l2, l3):
    x = sixJ(j1, j2, j3, l1, l2, l3)
    print x.evalf(), x, pr(x)