symbolic calculation of translation coefficients in SageMath

Former-commit-id: 13aa0ed33d8ecd85bcb4952cbd303a7d2bba4634

tests/Amaxcumsum.c

@@ -0,0 +1,25 @@
+#include "../qpms/gaunt.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+//const int lMax = 30;
+
+int main(int argc, char **argv) 
+{
+	int lMax;
+	if (argc < 2)
+		lMax = 30;
+	else lMax = atoi(argv[1]);
+
+	printf("// assuming lMax == %d:\n", lMax);
+	size_t Amaxcumsum = 0;
+	printf("0x%zx,\n", Amaxcumsum);
+	for (int n = 0; n <= lMax; n++)
+		for (int m = -n; m <= n; m++)
+			for (int nu = 0; nu <= lMax; nu++)
+				for (int mu = -nu; mu <= nu; mu++) {
+					Amaxcumsum += gaunt_q_max(-m, n, mu, nu) + 1;
+					printf("0x%zx,\n", Amaxcumsum);
+				}
+	return 0;
+}

tests/Bmaxcumsum.c

@@ -0,0 +1,34 @@
+#include "../qpms/gaunt.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+//const int lMax = 30;
+
+#define MIN(x,y) (((x)>(y))?(y):(x))
+
+static inline int Qmax(int m, int n, int mu, int nu) {
+	return MIN(n, MIN(nu, (n+nu+1-abs(-m+mu))/2));
+}
+
+int main(int argc, char **argv) 
+{
+	int lMax;
+	if (argc < 2)
+		lMax = 30;
+	else lMax = atoi(argv[1]);
+
+	printf("// assuming lMax == %d:\n", lMax);
+	size_t Bmaxcumsum = 0;
+	printf("0x%zx,\n", Bmaxcumsum);
+	for (int n = 0; n <= lMax; n++)
+		for (int m = -n; m <= n; m++)
+			for (int nu = 0; nu <= lMax; nu++)
+				for (int mu = -nu; mu <= nu; mu++) {
+					//FIXME toto je možná úplně blbě
+					assert(gaunt_q_max(-m, n+1, mu, nu) == Qmax(m, n, mu, nu);
+					Bmaxcumsum += gaunt_q_max(-m, n+1, mu, nu) + 1;
+					printf("0x%zx,\n", Bmaxcumsum);
+				}
+	return 0;
+}

tests/transcoeff_cruzan.py

@@ -0,0 +1,92 @@
+# [Xu] = Journal of computational physics 139, 137–165
+
+from __future__ import print_function
+
+def p_q(q, n, nu):
+    return n + nu - 2*q
+
+def qmax(M, n, mu, nu):
+    return floor(min(n, nu, (n+nu-abs(M+mu))/2))
+
+def Qmax(M, n, mu, nu): # [Xu](60)
+    return floor(min(n, nu, (n+nu+1-abs(M+mu))/2))
+
+def gaunta_p(M, n, mu, nu, p): # [Xu](5)
+    #print (M,n,mu,nu,p, file=sys.stderr)
+    return (-1)**(M+mu) * (2*p +1) * sqrt(
+        factorial(n+M) * factorial(nu+mu) * factorial(p-M-mu)
+      / factorial(n-M) / factorial(nu-mu) / factorial(p+M+mu)) * (
+        wigner_3j(n, nu, p, 0, 0, 0) * wigner_3j(n, nu, p, M, mu, -M-mu))
+
+def bCXcoeff(M, n, mu, nu, p): # [Xu](61)
+    return (-1)**(M+mu) * (2*p + 3) * sqrt(
+        factorial(n+M) * factorial(nu+mu) * factorial(p+1-M-mu)
+      / factorial(n-M) / factorial(nu-mu) / factorial(p+1+M+mu)) * (
+        wigner_3j(n, nu, p, 0, 0, 0) * wigner_3j(n, nu, p+1, M, mu, -M-mu))
+
+def ACXcoeff(m, n, mu, nu, q): # [Xu](58)
+    p = p_q(q,n,nu)
+    return ((-1)**m * (2*nu + 1) * factorial(n+m) * factorial(nu-mu) / (
+        2 * n * (n+1) * factorial(n-m) * factorial(nu+mu)) * I**p * 
+        (n*(n+1) + nu*(nu+1) - p*(p+1)) * gaunta_p(-m,n,mu,nu,p))
+
+def BCXcoeff(m, n, mu, nu, q): # [Xu](59)
+    p = p_q(q,n,nu)
+    return ((-1)**(m+1) * (2*nu + 1) * factorial(nu+m) * factorial(nu-mu) / (
+        2 * n * (n+1) * factorial(n-m) * factorial(nu+mu)) * I**(p+1) * 
+        sqrt(((p+1)**2-(n-nu)**2) * ((n+nu+1)**2-(p+1)**2)) 
+        * gaunta_p(-m,n,mu,nu,p))
+
+def printACXcoeffs(lMax, file=sys.stdout):
+    for n in IntegerRange(lMax+1):
+        for nu in IntegerRange(lMax+1):
+            for m in IntegerRange(-n, n+1):
+                for mu in IntegerRange(-nu, nu+1):
+                    for q in IntegerRange(qmax(-m,n,mu,nu)):
+                        #print(m, n, mu, nu, q, p_q(q,n,nu), file=sys.stderr)
+                        coeff= ACXcoeff(m, n, mu, nu, q);
+                        print(N(ACXcoeff(m, n, mu, nu, q), prec=53), 
+                            ", // %d, %d, %d, %d, %d," % (m,n,mu,nu,q),
+                            coeff,
+                            file=file)
+    return
+
+def printBCXcoeffs(lMax, file=sys.stdout):
+    for n in IntegerRange(lMax+1):
+        for nu in IntegerRange(lMax+1):
+            for m in IntegerRange(-n, n+1):
+                for mu in IntegerRange(-nu, nu+1):
+                    for q in IntegerRange(1, Qmax(-m,n,mu,nu)):
+                        #print(m, n, mu, nu, q, p_q(q,n,nu), file=sys.stderr)
+                        coeff= BCXcoeff(m, n, mu, nu, q);
+                        print(N(BCXcoeff(m, n, mu, nu, q), prec=53), 
+                            ", // %d, %d, %d, %d, %d," % (m,n,mu,nu,q),
+                            coeff,
+                            file=file)
+    return
+
+sphericalBessels = (None,
+        spherical_bessel_J,
+        spherical_bessel_Y,
+        spherical_hankel1,
+        spherical_hankel2
+        )
+
+# N.B. sage's gen_legendre_P _does_ include (-1)**m Condon-Shortley phase
+# whereas formulae in [Xu] do not.
+def trcoeff_ACX(m, n, mu, nu, besseltype, kd, th, fi, csphase=1): # [Xu](58)
+    res = 0
+    for q in range(qmax(-m,n,mu,nu)+1):
+        p = p_q(q,n,nu)
+        res += ACXcoeff(m,n,mu,nu,q) * sphericalBessels[besseltype](p,kd) * gen_legendre_P(p, mu-m, cos(th)) * (-csphase)**(mu-m) # compensate for csphase
+    res *= exp(I*(mu-m)*fi)
+    return res
+
+def trcoeff_BCX(m, n, mu, nu, besseltype, kd, th, fi, csphase=1): # [Xu](59) 
+    res = 0
+    for q in range(Qmax(-m,n,mu,nu)+1):
+        p = p_q(q,n,nu)
+        res += BCXcoeff(m,n,mu,nu,q) * sphericalBessels[besseltype](p+1) * gen_legendre_P(p+1, mu-m, cos(th)) * (-csphase)**(mu-m)
+    res *= exp(I*(mu-m)*fi)
+    return res
+