From 6474ceff0b7d22b699528149323cbacb00f007ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 7 Mar 2018 18:44:32 +0200 Subject: [PATCH] =?UTF-8?q?Testing=20translation=20coefficients=20?= =?UTF-8?q?=E2=80=93=20found=20some=20irregularities.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: 1f4791f3c6447de48229ccc4eb5a59c6cd7a2968 --- notes/tftests.lyx | 93 +++++++++++++++++++++++++ qpms/test_translations_xu_multipliers.c | 46 ++++++++++++ tests/cruzan_A.m | 26 +++++++ tests/cruzan_B.m | 26 +++++++ tests/transcoeff_cruzan.py | 2 +- 5 files changed, 192 insertions(+), 1 deletion(-) create mode 100644 qpms/test_translations_xu_multipliers.c create mode 100644 tests/cruzan_A.m create mode 100644 tests/cruzan_B.m diff --git a/notes/tftests.lyx b/notes/tftests.lyx index e9f2cdd..359d623 100644 --- a/notes/tftests.lyx +++ b/notes/tftests.lyx @@ -312,6 +312,99 @@ Vytvoř náhodně vektor přesunu . +\end_layout + +\begin_layout Enumerate +TODO +\end_layout + +\begin_layout Subsection +Nalezené nesrovnalosti +\end_layout + +\begin_layout Standard +Xuovy vzorce ve starší práci +\begin_inset CommandInset citation +LatexCommand cite +after "(77–80)" +key "xu_calculation_1996" + +\end_inset + + a v novější práci +\begin_inset CommandInset citation +LatexCommand cite +after "(63–65, ...)" +key "xu_efficient_1998" + +\end_inset + + pro koefficienty +\begin_inset Formula $B_{mn\mu\nu}$ +\end_inset + + se liší v několika ohledech: +\end_layout + +\begin_layout Enumerate +Ve starší práci suma začíná na +\begin_inset Formula $q=0$ +\end_inset + +, kdežto v novější práci až na +\begin_inset Formula $q=1$ +\end_inset + +. + Ovšem členy s +\begin_inset Formula $q=0$ +\end_inset + + jsou identicky nulové, takže je zbytečné začínat na nule. + (Ověřeno numericky – i tam jsou to přesně nuly.) +\end_layout + +\begin_layout Enumerate +Ve starší práci je poslední člen sumy +\begin_inset Formula $q=\min\left(n+1,\nu,\frac{n+\nu+1-\left|\mu-m\right|}{2}\right)$ +\end_inset + +, zatímco v novější je to +\begin_inset Formula $q=\min\left(n,\nu,\frac{n+\nu+1-\left|\mu-m\right|}{2}\right)$ +\end_inset + +. + Tyto hodnoty se pochopitelně mohou lišit, například pro +\begin_inset Formula $\left(m,n,\mu,\nu\right)=\left(-1,1,-1,3\right)$ +\end_inset + +. + Numericky ověřeno, že „přebytečné“ členy ze starší práce jsou nulové (avšak + vypočtené hodnoty nejsou přesně nuly, něco zbude kvůli zaokrouhlovacích + chyb). +\end_layout + +\begin_layout Enumerate +!!! Některé hodnoty nesedějí, například pro +\begin_inset Formula $\left(m,n,\mu,\nu\right)=\left(0,1,-1,1\right)$ +\end_inset + +!!! +\end_layout + +\begin_layout Enumerate +A nakonec samotné vzorce pro sčítance mají poněkud jiný tvar. +\end_layout + +\begin_layout Standard +\begin_inset CommandInset bibtex +LatexCommand bibtex +bibfiles "/l/necadam1/repo/qpms/Electrodynamics" +options "plain" + +\end_inset + + \end_layout \end_body diff --git a/qpms/test_translations_xu_multipliers.c b/qpms/test_translations_xu_multipliers.c new file mode 100644 index 0000000..43e678b --- /dev/null +++ b/qpms/test_translations_xu_multipliers.c @@ -0,0 +1,46 @@ +#include "translations.h" +#include "gaunt.h" +#include +//#include +#include + + +void qpms_trans_calculator_multipliers_B_general( + qpms_normalisation_t norm, + complex double *dest, int m, int n, int mu, int nu, int Qmax) ; + +int lMax=13; + +#define MIN(x,y) ((x)<(y)?(x):(y)) + +// Python test: Qmax(M, n, mu, nu) = floor(min(n,nu,(n+nu+1-abs(M+mu))/2)) +// q in IntegerRange(1, Qmax(-m,n,mu,nu)) + +int main() { + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_XU); + complex double dest[lMax + 2]; + + for(int n = 1; n <= lMax; ++n) + for(int nu = 1; nu <= lMax; ++nu) + for(int m = -n; m <= n; ++m) + for(int mu = -nu; mu <= nu; ++mu){ + int Qmax = gaunt_q_max(-m, n+1, mu, nu); + int Qmax_alt = MIN(n,MIN(nu,(n+nu+1-abs(mu-m)))); + qpms_trans_calculator_multipliers_B_general(QPMS_NORMALISATION_XU_CS, + dest, m, n, mu, nu, Qmax); + for(int q = 0; q <= Qmax; ++q) { + // int p = n + nu - 2*q; + int tubig = cabs(dest[q]) > 1e-8; + printf("%.16g + %.16g*I, // %d, %d, %d, %d, %d,%s\n", + creal(dest[q]), cimag(dest[q]), m, n, mu, nu, q, + q > Qmax_alt ? (tubig?" //tubig":" //tu") : ""); + } + fflush(stdout); + } + + + + qpms_trans_calculator_free(c); +} + + diff --git a/tests/cruzan_A.m b/tests/cruzan_A.m new file mode 100644 index 0000000..59c41be --- /dev/null +++ b/tests/cruzan_A.m @@ -0,0 +1,26 @@ +(* Vector translation coefficients as in Journal of Computational Physics 139, 137--165 (1998), eqs. (58), (59), (61) *) +gaunt[m_, n_, mu_, nu_, + p_] := (-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]] ThreeJSymbol[{n, + 0}, {nu, 0}, {p, 0}] ThreeJSymbol[{n, m}, {nu, + mu}, {p, -m - mu}]; +bCXcoeff[m_,n_,mu_,nu_,p_]:=(-1)^(mu+m)(2p+3)Sqrt[(n+m)!(nu+mu)!(p-m-mu+1)!/(n-m)!/(nu-mu)!/(p+m+mu+1)!]ThreeJSymbol[{n,m},{nu,mu},{p+1,-m-mu}]ThreeJSymbol[{n,0},{nu,0},{p,0}] +p[q_,n_,nu_]:=n+nu-2q; +ACXcoeff[m_,n_,mu_,nu_,q_]:=(-1)^m (2nu+1)(nu+m)!(nu-mu)!/2/n/(nu+1)/(nu-m)!/(nu+m)!I^p[q,n,nu](n(n+1)+nu(nu+1)-p[q,n,nu](p[q,n,nu]+1))gaunt[-m,n,mu,nu,p[q,n,nu]] +BCXcoeff[m_,n_,mu_,nu_,q_]:=(-1)^(m+1)(2nu+1)(n+m)!(nu-m)!/2/n/(n+1)/(n-m)!(nu+mu)!I^(p[q,n,nu]+1)Sqrt[((p[q,n,nu]+1)^2-(n-nu)^2)((n+nu+1)^2-(p[q,n,nu]+1)^2)]bCXcoeff[-m,n,mu,nu,p[q,n,nu]]` + +lMax := 5 +For[n = 0, n <= lMax, n++, + For[nu = 0, nu <= lMax, nu++, + For[m = -n, m <= n, m++, + For[mu = -nu, mu <= nu, mu++, + For[q = 0, q <= Min[n, nu, (n + nu - Abs[m + mu])/2],q++, + Print[CForm[N[ACXcoeff[m,n,mu,nu,q],16]]] + ] + ] + ] + ] + ] + diff --git a/tests/cruzan_B.m b/tests/cruzan_B.m new file mode 100644 index 0000000..481bfd3 --- /dev/null +++ b/tests/cruzan_B.m @@ -0,0 +1,26 @@ +(* Vector translation coefficients as in Journal of Computational Physics 139, 137--165 (1998), eqs. (58), (59), (61) *) +gaunt[m_, n_, mu_, nu_, + p_] := (-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]] ThreeJSymbol[{n, + 0}, {nu, 0}, {p, 0}] ThreeJSymbol[{n, m}, {nu, + mu}, {p, -m - mu}]; +bCXcoeff[m_,n_,mu_,nu_,p_]:=(-1)^(mu+m)(2p+3)Sqrt[(n+m)!(nu+mu)!(p-m-mu+1)!/(n-m)!/(nu-mu)!/(p+m+mu+1)!]ThreeJSymbol[{n,m},{nu,mu},{p+1,-m-mu}]ThreeJSymbol[{n,0},{nu,0},{p,0}] +p[q_,n_,nu_]:=n+nu-2q; +ACXcoeff[m_,n_,mu_,nu_,q_]:=(-1)^m (2nu+1)(nu+m)!(nu-mu)!/2/n/(nu+1)/(nu-m)!/(nu+m)!I^p[q,n,nu](n(n+1)+nu(nu+1)-p[q,n,nu](p[q,n,nu]+1))gaunt[-m,n,mu,nu,p[q,n,nu]] +BCXcoeff[m_,n_,mu_,nu_,q_]:=(-1)^(m+1)(2nu+1)(n+m)!(nu-m)!/2/n/(n+1)/(n-m)!(nu+mu)!I^(p[q,n,nu]+1)Sqrt[((p[q,n,nu]+1)^2-(n-nu)^2)((n+nu+1)^2-(p[q,n,nu]+1)^2)]bCXcoeff[-m,n,mu,nu,p[q,n,nu]] + +lMax := 18 +For[n = 0, n <= lMax, n++, + For[nu = 0, nu <= lMax, nu++, + For[m = -n, m <= n, m++, + For[mu = -nu, mu <= nu, mu++, + For[q = 1, q <= Min[n, nu, (n + nu + 1 - Abs[-m + mu])/2],q++, + Print[CForm[N[BCXcoeff[m,n,mu,nu,q],16]]] + ] + ] + ] + ] + ] + diff --git a/tests/transcoeff_cruzan.py b/tests/transcoeff_cruzan.py index f4ae08a..912c497 100644 --- a/tests/transcoeff_cruzan.py +++ b/tests/transcoeff_cruzan.py @@ -57,7 +57,7 @@ def printBCXcoeffs(lMax, file=sys.stdout): 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)): + for q in IntegerRange(1, Qmax(-m,n,mu,nu) +1 ): #print(m, n, mu, nu, q, p_q(q,n,nu), file=sys.stderr) coeff= BCXcoeff(m, n, mu, nu, q); print(N(coeff, prec=53),