From 155859d8a7f56daae8b71d31ddfdc1ae6862f05b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 29 Jan 2018 17:33:59 +0200 Subject: [PATCH] Fix edge legendre signs again Former-commit-id: 932787e15c101c131bded3ada1034a379da6123c --- misc/test_vswf.gnuplot | 11 ++++++++++ notes/Scattering and Shit.lyx | 38 +++++++++++++++++------------------ qpms/qpms_types.h | 29 +++++++++++++++++++++++--- qpms/test_vswf.c | 8 ++++++-- qpms/vswf.c | 32 ++++++++++++++--------------- 5 files changed, 77 insertions(+), 41 deletions(-) create mode 100644 misc/test_vswf.gnuplot diff --git a/misc/test_vswf.gnuplot b/misc/test_vswf.gnuplot new file mode 100644 index 0000000..7b07259 --- /dev/null +++ b/misc/test_vswf.gnuplot @@ -0,0 +1,11 @@ +f = 'out' +fc = 'outcs' +do for [t=3:137] { +y = ((t-3) % 45)/3 +typ = (t-3) % 3 +n = floor(sqrt(y+1)) +m = y - (n*(n+1)-1) +print 'n = ', n, ', m = ', m, ', typ ', typ +plot f using 1:t w linespoints, fc using 1:t w linespoints +pause -1 +} diff --git a/notes/Scattering and Shit.lyx b/notes/Scattering and Shit.lyx index d332040..1e53c39 100644 --- a/notes/Scattering and Shit.lyx +++ b/notes/Scattering and Shit.lyx @@ -301,8 +301,8 @@ The expressions and \begin_inset Formula \begin{eqnarray*} -\pi_{1\nu}(1-) & = & -\frac{\nu\left(\nu+1\right)\left(\nu+2\right)}{2}\\ -\tau_{1\nu}(1-) & = & \frac{\nu\left(\nu+1\right)\left(\nu+2\right)}{2} +\pi_{1\nu}(+1-) & = & CS\frac{\nu\left(\nu+1\right)}{2}\\ +\tau_{1\nu}(+1-) & = & CS\frac{\nu\left(\nu+1\right)}{2} \end{eqnarray*} \end_inset @@ -314,8 +314,8 @@ and using the parity property \begin_inset Formula \begin{eqnarray*} -\pi_{1\nu}(-1+) & = & \left(-1\right)^{\nu}\frac{\nu\left(\nu+1\right)\left(\nu+2\right)}{2}\\ -\tau_{1\nu}(-1+) & = & \left(-1\right)^{\nu}\frac{\nu\left(\nu+1\right)\left(\nu+2\right)}{2} +\pi_{1\nu}(-1+) & = & -CS\left(-1\right)^{\nu}\frac{\nu\left(\nu+1\right)}{2}\\ +\tau_{1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\frac{\nu\left(\nu+1\right)}{2} \end{eqnarray*} \end_inset @@ -331,10 +331,10 @@ For to get \begin_inset Formula \begin{eqnarray*} -\pi_{-1\nu}(1-) & = & -CS\frac{\nu+2}{2}\\ -\tau_{-1\nu}(1-) & = & CS\frac{\nu+2}{2}\\ -\pi_{-1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\frac{\nu+2}{2}\\ -\tau_{-1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\frac{\nu+2}{2} +\pi_{-1\nu}(+1-) & = & \frac{CS}{2}\\ +\tau_{-1\nu}(+1-) & = & -\frac{CS}{2}\\ +\pi_{-1\nu}(-1+) & = & -\left(-1\right)^{\nu}\frac{CS}{2}\\ +\tau_{-1\nu}(-1+) & = & -\left(-1\right)^{\nu}\frac{CS}{2} \end{eqnarray*} \end_inset @@ -380,10 +380,10 @@ reference "sub:Xu pitau" by the normalisation factor, \begin_inset Formula \begin{eqnarray*} -\tilde{\pi}_{1\nu}(1-) & = & -\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\tau}_{1\nu}(1-) & = & \sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\pi}_{1\nu}(-1+) & = & \left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\tau}_{1\nu}(-1+) & = & \left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2} +\tilde{\pi}_{1\nu}(+1-) & = & CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2}\\ +\tilde{\tau}_{1\nu}(+1-) & = & CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2}\\ +\tilde{\pi}_{1\nu}(-1+) & = & -CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2}\\ +\tilde{\tau}_{1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2} \end{eqnarray*} \end_inset @@ -391,10 +391,10 @@ reference "sub:Xu pitau" \begin_inset Formula \begin{eqnarray*} -\tilde{\pi}_{-1\nu}(1-) & = & -CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\tau}_{-1\nu}(1-) & = & CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\pi}_{-1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ -\tilde{\tau}_{-1\nu}(-1+) & = & CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2} +\tilde{\pi}_{-1\nu}(+1-) & = & CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2}\\ +\tilde{\tau}_{-1\nu}(+1-) & = & -CS\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}}{2}\\ +\tilde{\pi}_{-1\nu}(-1+) & = & -CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2}\\ +\tilde{\tau}_{-1\nu}(-1+) & = & -CS\left(-1\right)^{\nu}\sqrt{\frac{2\nu+1}{4\pi}}\frac{\sqrt{\nu\left(\nu+1\right)}\left(\nu+2\right)}{2} \end{eqnarray*} \end_inset @@ -699,7 +699,7 @@ Definition [T](2.40); \begin_inset Formula \begin{eqnarray*} \widetilde{\vect N}_{mn}^{(j)} & = & \frac{n(n+1)}{kr}\sqrt{\frac{2n+1}{4\pi}\frac{\left(n-m\right)!}{\left(n+m\right)!}}P_{n}^{m}\left(\cos\theta\right)e^{im\phi}z_{n}^{j}\left(kr\right)\hat{\vect r}\\ - & & +\left[\tilde{\tau}_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}+i\tilde{\pi}_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}\frac{1}{kr}\frac{\ud\left(kr\, z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}\\ + & & +\left[\tilde{\tau}_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}+i\tilde{\pi}_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}\frac{1}{kr}\frac{\ud\left(kr\,z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}\\ \widetilde{\vect M}_{mn}^{(j)} & = & \left[i\tilde{\pi}_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}-\tilde{\tau}_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}z_{n}^{j}\left(kr\right) \end{eqnarray*} @@ -720,7 +720,7 @@ are the electric and magnetic waves, respectively: \begin_inset Formula \begin{eqnarray*} \vect N_{mn}^{(j)} & = & \frac{n(n+1)}{kr}P_{n}^{m}\left(\cos\theta\right)e^{im\phi}z_{n}^{j}\left(kr\right)\hat{\vect r}\\ - & & +\left[\tau_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}+i\pi_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}\frac{1}{kr}\frac{\ud\left(kr\, z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}\\ + & & +\left[\tau_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}+i\pi_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}\frac{1}{kr}\frac{\ud\left(kr\,z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}\\ \vect M_{mn}^{(j)} & = & \left[i\pi_{mn}\left(\cos\theta\right)\hat{\vect{\theta}}-\tau_{mn}\left(\cos\theta\right)\hat{\vect{\phi}}\right]e^{im\phi}z_{n}^{j}\left(kr\right) \end{eqnarray*} @@ -760,7 +760,7 @@ Definition [K](2.4.6); to simlify the notation. \begin_inset Formula \begin{eqnarray*} -\left(\vect{u/v/w}\right)_{2lm} & = & \frac{1}{kr}\frac{\ud\left(kr\, z_{l}^{(j)}\left(kr\right)\right)}{\ud\, kr}\vect A_{2lm}\left(\hat{\vect r}\right)+\sqrt{l\left(l+1\right)}\frac{z_{l}^{(j)}(kr)}{kr}\vect A_{3lm}\left(\hat{\vect r}\right)\\ +\left(\vect{u/v/w}\right)_{2lm} & = & \frac{1}{kr}\frac{\ud\left(kr\,z_{l}^{(j)}\left(kr\right)\right)}{\ud\,kr}\vect A_{2lm}\left(\hat{\vect r}\right)+\sqrt{l\left(l+1\right)}\frac{z_{l}^{(j)}(kr)}{kr}\vect A_{3lm}\left(\hat{\vect r}\right)\\ \left(\vect{u/v/w}\right)_{1lm} & = & z_{l}^{(j)}\left(kr\right)\vect A_{1lm}\left(\hat{\vect r}\right) \end{eqnarray*} diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index 4f743d6..bed9fd6 100644 --- a/qpms/qpms_types.h +++ b/qpms/qpms_types.h @@ -17,16 +17,39 @@ typedef enum { } qpms_errno_t; // Normalisations +//const int QPMS_NORMALISATION_T_CSBIT = 128; +#define QPMS_NORMALISATION_T_CSBIT 128 typedef enum { // As in TODO - QPMS_NORMALISATION_XU = 3, // NI! + QPMS_NORMALISATION_XU = 3, + QPMS_NORMALISATION_NONE = QPMS_NORMALISATION_XU, // As in http://www.eit.lth.se/fileadmin/eit/courses/eit080f/Literature/book.pdf, power-normalised - QPMS_NORMALISATION_KRISTENSSON = 2, // NI! - QPMS_NORMALISATION_POWER = QPMS_NORMALISATION_KRISTENSSON, // NI! + QPMS_NORMALISATION_KRISTENSSON = 2, + QPMS_NORMALISATION_POWER = QPMS_NORMALISATION_KRISTENSSON, + // as in TODO QPMS_NORMALISATION_TAYLOR = 1, + QPMS_NORMALISATION_SPHARM = QPMS_NORMALISATION_TAYLOR, + // Variants with Condon-Shortley phase + QPMS_NORMALISATION_XU_CS = QPMS_NORMALISATION_XU | QPMS_NORMALISATION_T_CSBIT, + QPMS_NORMALISATION_NONE_CS = QPMS_NORMALISATION_XU_CS, + QPMS_NORMALISATION_KRISTENSSON_CS = QPMS_NORMALISATION_KRISTENSSON | QPMS_NORMALISATION_T_CSBIT, + QPMS_NORMALISATION_POWER_CS = QPMS_NORMALISATION_KRISTENSSON_CS, + QPMS_NORMALISATION_TAYLOR_CS = QPMS_NORMALISATION_TAYLOR | QPMS_NORMALISATION_T_CSBIT, + QPMS_NORMALISATION_SPHARM_CS = QPMS_NORMALISATION_TAYLOR_CS, QPMS_NORMALISATION_UNDEF = 0 } qpms_normalisation_t; +static inline int qpms_normalisation_t_csphase(qpms_normalisation_t norm) { + return (norm & QPMS_NORMALISATION_T_CSBIT)? -1 : 1; +} + +static inline int qpms_normalisation_t_normonly(qpms_normalisation_t norm) { + return norm & (~QPMS_NORMALISATION_T_CSBIT); +} + + + + typedef enum { QPMS_BESSEL_REGULAR = 1, // regular function j QPMS_BESSEL_SINGULAR = 2, // singular function y diff --git a/qpms/test_vswf.c b/qpms/test_vswf.c index 12d145a..6103851 100644 --- a/qpms/test_vswf.c +++ b/qpms/test_vswf.c @@ -4,12 +4,16 @@ #include const double dtheta = 0.01 * M_PI; const qpms_l_t lMax = 3; - +#ifdef CS +#define CS_OFFSET QPMS_NORMALISATION_T_CSBIT +#else +#define CS_OFFSET 0 +#endif int main() { qpms_y_t nelem = qpms_lMax2nelem(lMax); for (double theta = 0.; theta <= M_PI; theta += dtheta) { printf("%.5e %.5e ", theta, cos(theta)); - for(qpms_normalisation_t norm = 1; norm <= 3; ++norm) {//fujka :D + for(qpms_normalisation_t norm = CS_OFFSET + 1; norm <= CS_OFFSET + 3; ++norm) {//fujka :D qpms_pitau_t pt = qpms_pitau_get(theta, lMax, norm); for (qpms_y_t y = 0; y < nelem; ++y) printf("%.5e %.5e %.5e ", pt.leg[y], pt.pi[y], pt.tau[y]); diff --git a/qpms/vswf.c b/qpms/vswf.c index c58ed19..4837490 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -6,10 +6,6 @@ #include #include -#ifndef CSPHASE -#define CSPHASE (1.) // FIXME this should be later determined by qpms_normalization_t -#endif - // Legendre functions also for negative m, see DLMF 14.9.3 qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase) @@ -72,6 +68,8 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **target, double **dtarget, double qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t norm) { + const double csphase = qpms_normalisation_t_csphase(norm); + norm = qpms_normalisation_t_normonly(norm); qpms_pitau_t res; qpms_y_t nelem = qpms_lMax2nelem(lMax); res.pi = malloc(nelem * sizeof(double)); @@ -88,10 +86,10 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no double p = l*(l+1)/2; const double n = 0.5; int lpar = (l%2)?-1:1; - res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * p; - res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * n * CSPHASE; - res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p; - res.tau[qpms_mn2y(-1, l)] = ((ct>0) ? +1 : lpar) * n * CSPHASE; + res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * p * csphase; + res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * n * csphase; + res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * p * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * n * csphase; } break; case QPMS_NORMALISATION_TAYLOR: @@ -99,10 +97,10 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)*0.25*M_1_PI); int lpar = (l%2)?-1:1; double fl = 0.25 * sqrt((2*l+1)*l*(l+1)*M_1_PI); - res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * fl; - res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; - res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl; - res.tau[qpms_mn2y(-1, l)] = ((ct>0) ? +1 : lpar) * fl * CSPHASE; + res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; } break; case QPMS_NORMALISATION_POWER: @@ -110,10 +108,10 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no res.leg[qpms_mn2y(0, l)] = ((l%2)?ct:1.)*sqrt((2*l+1)/(4*M_PI *l*(l+1))); int lpar = (l%2)?-1:1; double fl = 0.25 * sqrt((2*l+1)*M_1_PI); - res.pi [qpms_mn2y(+1, l)] = ((ct>0) ? -1 : lpar) * fl; - res.pi [qpms_mn2y(-1, l)] = ((ct>0) ? -1 : lpar) * fl * CSPHASE; - res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl; - res.tau[qpms_mn2y(-1, l)] = ((ct>0) ? +1 : lpar) * fl * CSPHASE; + res.pi [qpms_mn2y(+1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + res.pi [qpms_mn2y(-1, l)] = -((ct>0) ? -1 : lpar) * fl * csphase; + res.tau[qpms_mn2y(+1, l)] = ((ct>0) ? +1 : lpar) * fl * csphase; + res.tau[qpms_mn2y(-1, l)] = -((ct>0) ? +1 : lpar) * fl * csphase; } break; @@ -126,7 +124,7 @@ qpms_pitau_t qpms_pitau_get(double theta, qpms_l_t lMax, qpms_normalisation_t no res.leg = malloc(sizeof(double)*qpms_lMax2nelem(lMax)); if (qpms_legendre_deriv_y_fill(res.leg, legder, ct, lMax, norm == QPMS_NORMALISATION_XU ? GSL_SF_LEGENDRE_NONE - : GSL_SF_LEGENDRE_SPHARM, CSPHASE)) + : GSL_SF_LEGENDRE_SPHARM, csphase)) abort(); if (norm == QPMS_NORMALISATION_POWER) /* for Xu (=non-normalized) and Taylor (=sph. harm. normalized)