QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fix edge legendre signs again 

Former-commit-id: 932787e15c101c131bded3ada1034a379da6123c
This commit is contained in:
Marek Nečada 2018-01-29 17:33:59 +02:00
parent bd5aacd7d9
commit 155859d8a7
5 changed files with 77 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
misc/test_vswf.gnuplot Normal file
View File

@ -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
}

38
notes/Scattering and Shit.lyx
View File

@ -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*}

29
qpms/qpms_types.h
View File

@ -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

8
qpms/test_vswf.c
View File

@ -4,12 +4,16 @@
#include <gsl/gsl_math.h>
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]);

32
qpms/vswf.c
View File

@ -6,10 +6,6 @@
#include <stdlib.h>
#include <string.h>
#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)