From a9580b7fd4be82dcf17300795d404f6e6152f9e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 25 Feb 2018 16:30:42 +0000 Subject: [PATCH] symbolic calculation of translation coefficients in SageMath Former-commit-id: 13aa0ed33d8ecd85bcb4952cbd303a7d2bba4634 --- tests/Amaxcumsum.c | 25 +++++++++++ tests/Bmaxcumsum.c | 34 ++++++++++++++ tests/transcoeff_cruzan.py | 92 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+) create mode 100644 tests/Amaxcumsum.c create mode 100644 tests/Bmaxcumsum.c create mode 100644 tests/transcoeff_cruzan.py diff --git a/tests/Amaxcumsum.c b/tests/Amaxcumsum.c new file mode 100644 index 0000000..63e5467 --- /dev/null +++ b/tests/Amaxcumsum.c @@ -0,0 +1,25 @@ +#include "../qpms/gaunt.h" +#include +#include + +//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; +} diff --git a/tests/Bmaxcumsum.c b/tests/Bmaxcumsum.c new file mode 100644 index 0000000..11ff352 --- /dev/null +++ b/tests/Bmaxcumsum.c @@ -0,0 +1,34 @@ +#include "../qpms/gaunt.h" +#include +#include +#include + +//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; +} diff --git a/tests/transcoeff_cruzan.py b/tests/transcoeff_cruzan.py new file mode 100644 index 0000000..c14a386 --- /dev/null +++ b/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 +