wpoely86 · March 2, 2015 11:10
diff --git a/Clebsch-Gordan.py b/Clebsch-Gordan.py
 #!/usr/bin/env python3
 """
 Generates all Clebsh Gordan coefficients starting from a given j1 and j2 and
 outputs a html file. How to use:
 ./Clebsch-Gordan.py 0.5 0.5 > output.html
 Ward Poelmans <[email protected]> (2014)
 """

 import sys
 import numpy as np
 import scipy.misc as sc
 import math
 import fractions


 def factorial(n):
    if n < 0:
        return -100000000000
    else:
        return math.factorial(n)
 #        return sc.factorial(n, exact=False)


 def clebsh(j1,j2,j,m):
    z = 0
    a = 0.0
    b = 0.0
    c = 0.0
    l = ""
    for m1 in np.arange(-j1-1,(j1+1)):
        for m2 in np.arange(-j2-1,(j2+1)):
            #These are just some conditions which need to give zero.
            if (m1 + m2) != m:
                d = 0
            elif (j1 - j2 - m)%1 != 0.0:
                d = 0
            elif (j1 - j2 + m)%1 != 0.0:
                d = 0
            else:
                c = 0.0
                z = 0
                #This is the algorithm itself, Not very pretty, but that is what it is.
                at = ((factorial(j1+j2-j)*factorial(j1-j2+j)*factorial(-j1+j2+j)))
                an = (factorial(j1+j2+j+1.0)) #**.5
                a = math.sqrt(at/an)
                b = (factorial(j1+m1)*factorial(j1-m1)*factorial(j2+m2)*factorial(j2-m2)*factorial(j-m)*factorial(j+m)) #**.5
                while z<(j1-m1+3):
                    ct = ((-1.0)**(z+j1-j2+m))
                    cn = (factorial(z)*factorial(j1+j2-j-z)*factorial(j1-m1-z)*factorial(j2+m2-z)*factorial(j-j2+m1+z)*factorial(j-j1-m2+z))
                    c += ct/cn;
                    z += 1
                pre = ((-1.0)**(j1-j2+m))
                presqrt = (2.0*j+1.0)
                d = pre*math.sqrt(presqrt)*math.sqrt(at/an)*math.sqrt(b)*c
                if abs(d)> .0001:
                    breuk = fractions.Fraction(int(at*b*presqrt*c*c),int(an)).limit_denominator(10000)
                    pre_str = "-" if pre < 0 else ""
                    l += "\n<tr><td><math><mn>%s</mn></math></td><td><math><mn>%s</mn></math></td><td><math><mn>%s</mn> <msqrt><mfrac><mn>%s</mn><mn>%s</mn></mfrac></msqrt></math></td></tr>" % (m1,m2,pre_str,breuk.numerator,breuk.denominator)
    return l


 if len(sys.argv) != 3:
    print("Usage: %s j1 j2" % sys.argv[0])
    sys.exit(0)

 j1 = float(sys.argv[1])
 j2 = float(sys.argv[2])

 print("<html>\n<title>Clebsch Gordan</title>\n<body>\n<div>\n")

 for j in np.arange(abs(j1-j2),j1+j2+1):
    for m in np.arange(0,j+1):
        print("<math><msub><mi>j</mi><mn>1</mn></msub>=<mn>%s</mn></math> <math><msub><mi>j</mi><mn>2</mn></msub>=<mn>%s</mn></math>  <math>j=<mn>%s</mn></math>  <math>m=<mn>%s</mn></math><br/>\n" % (fractions.Fraction(j1), fractions.Fraction(j2), fractions.Fraction(j), fractions.Fraction(m)))
        print("<table border=1>")
        print("<thead><tr> <th><math><msub><mi>m</mi><mn>1</mn></msub></math></th> <th><math><msub><mi>m</mi><mn>2</mn></msub></math></th> </tr></thead>\n")
        print("%s" % clebsh(j1,j2,j,m))
        print("</table><br/><br/>\n")

 print("</div>\n</body>\n</html>")
	#!/usr/bin/env python3
	"""
	Generates all Clebsh Gordan coefficients starting from a given j1 and j2 and
	outputs a html file. How to use:
	./Clebsch-Gordan.py 0.5 0.5 > output.html
	Ward Poelmans <[email protected]> (2014)
	"""

	import sys
	import numpy as np
	import scipy.misc as sc
	import math
	import fractions


	def factorial(n):
	if n < 0:
	return -100000000000
	else:
	return math.factorial(n)
	# return sc.factorial(n, exact=False)


	def clebsh(j1,j2,j,m):
	z = 0
	a = 0.0
	b = 0.0
	c = 0.0
	l = ""
	for m1 in np.arange(-j1-1,(j1+1)):
	for m2 in np.arange(-j2-1,(j2+1)):
	#These are just some conditions which need to give zero.
	if (m1 + m2) != m:
	d = 0
	elif (j1 - j2 - m)%1 != 0.0:
	d = 0
	elif (j1 - j2 + m)%1 != 0.0:
	d = 0
	else:
	c = 0.0
	z = 0
	#This is the algorithm itself, Not very pretty, but that is what it is.
	at = ((factorial(j1+j2-j)factorial(j1-j2+j)factorial(-j1+j2+j)))
	an = (factorial(j1+j2+j+1.0)) #**.5
	a = math.sqrt(at/an)
	b = (factorial(j1+m1)factorial(j1-m1)factorial(j2+m2)factorial(j2-m2)factorial(j-m)factorial(j+m)) #*.5
	while z<(j1-m1+3):
	ct = ((-1.0)**(z+j1-j2+m))
	cn = (factorial(z)factorial(j1+j2-j-z)factorial(j1-m1-z)factorial(j2+m2-z)factorial(j-j2+m1+z)*factorial(j-j1-m2+z))
	c += ct/cn;
	z += 1
	pre = ((-1.0)**(j1-j2+m))
	presqrt = (2.0*j+1.0)
	d = premath.sqrt(presqrt)math.sqrt(at/an)math.sqrt(b)c
	if abs(d)> .0001:
	breuk = fractions.Fraction(int(atbpresqrtcc),int(an)).limit_denominator(10000)
	pre_str = "-" if pre < 0 else ""
	l += "\n<tr><td><math><mn>%s</mn></math></td><td><math><mn>%s</mn></math></td><td><math><mn>%s</mn> <msqrt><mfrac><mn>%s</mn><mn>%s</mn></mfrac></msqrt></math></td></tr>" % (m1,m2,pre_str,breuk.numerator,breuk.denominator)
	return l


	if len(sys.argv) != 3:
	print("Usage: %s j1 j2" % sys.argv[0])
	sys.exit(0)

	j1 = float(sys.argv[1])
	j2 = float(sys.argv[2])

	print("<html>\n<title>Clebsch Gordan</title>\n<body>\n<div>\n")

	for j in np.arange(abs(j1-j2),j1+j2+1):
	for m in np.arange(0,j+1):
	print("<math><msub><mi>j</mi><mn>1</mn></msub>=<mn>%s</mn></math> <math><msub><mi>j</mi><mn>2</mn></msub>=<mn>%s</mn></math> <math>j=<mn>%s</mn></math> <math>m=<mn>%s</mn></math><br/>\n" % (fractions.Fraction(j1), fractions.Fraction(j2), fractions.Fraction(j), fractions.Fraction(m)))
	print("<table border=1>")
	print("<thead><tr> <th><math><msub><mi>m</mi><mn>1</mn></msub></math></th> <th><math><msub><mi>m</mi><mn>2</mn></msub></math></th> </tr></thead>\n")
	print("%s" % clebsh(j1,j2,j,m))
	print("</table><br/><br/>\n")

	print("</div>\n</body>\n</html>")