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

Rozbil jsem to, jdu domů. 

Former-commit-id: 77490824d9a8f975ec4ebe98d8cfa104093884e4
This commit is contained in:
Marek Nečada 2018-01-30 01:09:48 +02:00
parent 155859d8a7
commit 2a7015b80b
9 changed files with 209 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
notes/Scattering and Shit.lyx
View File

@ -683,6 +683,19 @@ In the usual cases we encounter, the part described by the magnetic waves
is pretty small. is pretty small.
\end_layout \end_layout
\begin_layout Standard
The expression with Bessel function derivatives appearing below in the electric
waves can be rewritten using (DLMF 10.51.2)
\begin_inset Formula
\[
\frac{1}{kr}\frac{\ud\left(kr\,z_{n}^{j}\left(kr\right)\right)}{\ud(kr)}=\frac{\ud z_{n}^{j}\left(kr\right)}{\ud(kr)}+\frac{z_{n}^{j}\left(kr\right)}{kr}=z_{n-1}^{j}\left(kr\right)-n\frac{z_{n}^{j}\left(kr\right)}{kr}.
\]
\end_inset
\end_layout
\begin_layout Subsection \begin_layout Subsection
Taylor Taylor
\end_layout \end_layout

38
qpms/bessel.c Normal file
View File

@ -0,0 +1,38 @@
#include "translations.h"
#include <stdlib.h>
#include <gsl/gsl_sf_bessel.h>
#include <assert.h>
int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
int retval;
double tmparr[lmax+1];
switch(typ) {
case QPMS_BESSEL_REGULAR:
retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
return retval;
break;
case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one?
retval = gsl_sf_bessel_yl_array(lmax,x,tmparr);
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
return retval;
break;
case QPMS_HANKEL_PLUS:
case QPMS_HANKEL_MINUS:
retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr);
for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l];
if(retval) return retval;
retval = gsl_sf_bessel_yl_array(lmax, x, tmparr);
if (typ==QPMS_HANKEL_PLUS)
for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l];
else
for (int l = 0; l <= lmax; ++l) result_array[l] +=-I * tmparr[l];
return retval;
break;
default:
abort();
//return GSL_EDOM;
}
assert(0);
}

6
qpms/qpms_types.h
View File

@ -63,6 +63,10 @@ typedef struct {
double x, y, z; double x, y, z;
} cart3_t; } cart3_t;
typedef struct {
complex double x, y, z;
} ccart3_t;
typedef struct { typedef struct {
double x, y; double x, y;
} cart2_t; } cart2_t;
@ -80,4 +84,6 @@ typedef struct {
double r, phi; double r, phi;
} pol_t; } pol_t;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#endif // QPMS_TYPES #endif // QPMS_TYPES

0
qpms/test_vswf.c → qpms/test_vswf_pitau.c
View File

57
qpms/translations.c
View File

@ -1,4 +1,5 @@
#include <math.h> #include <math.h>
#include "qpms_types.h"
#include "gaunt.h" #include "gaunt.h"
#include "translations.h" #include "translations.h"
#include "indexing.h" // TODO replace size_t and int with own index types here #include "indexing.h" // TODO replace size_t and int with own index types here
@ -34,7 +35,7 @@ double qpms_legendreD0(int m, int n) {
return -2 * qpms_legendre0(m, n); return -2 * qpms_legendre0(m, n);
} }
int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double *result_array) { int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) {
int retval; int retval;
double tmparr[lmax+1]; double tmparr[lmax+1];
switch(typ) { switch(typ) {
@ -86,7 +87,7 @@ complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj,
double leg[gsl_sf_legendre_array_n(n+nu)]; double leg[gsl_sf_legendre_array_n(n+nu)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
complex double bes[n+nu+1]; complex double bes[n+nu+1];
if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort(); if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
complex double sum = 0; complex double sum = 0;
for(int q = 0; q <= qmax; ++q) { for(int q = 0; q <= qmax; ++q) {
int p = n+nu-2*q; int p = n+nu-2*q;
@ -132,7 +133,7 @@ complex double qpms_trans_single_A_Kristensson(int m, int n, int mu, int nu, sph
double leg[gsl_sf_legendre_array_n(n+nu)]; double leg[gsl_sf_legendre_array_n(n+nu)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
complex double bes[n+nu+1]; complex double bes[n+nu+1];
if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort(); if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
complex double sum = 0; complex double sum = 0;
for(int q = 0; q <= qmax; ++q) { for(int q = 0; q <= qmax; ++q) {
int p = n+nu-2*q; int p = n+nu-2*q;
@ -177,7 +178,7 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
double leg[gsl_sf_legendre_array_n(n+nu)]; double leg[gsl_sf_legendre_array_n(n+nu)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
complex double bes[n+nu+1]; complex double bes[n+nu+1];
if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort(); if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort();
complex double sum = 0; complex double sum = 0;
for(int q = 0; q <= qmax; ++q) { for(int q = 0; q <= qmax; ++q) {
int p = n+nu-2*q; int p = n+nu-2*q;
@ -229,7 +230,7 @@ complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj,
double leg[gsl_sf_legendre_array_n(n+nu+1)]; double leg[gsl_sf_legendre_array_n(n+nu+1)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
complex double bes[n+nu+2]; complex double bes[n+nu+2];
if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort(); if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
complex double sum = 0; complex double sum = 0;
for (int q = 0; q <= Qmax; ++q) { for (int q = 0; q <= Qmax; ++q) {
@ -285,7 +286,7 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd
double leg[gsl_sf_legendre_array_n(n+nu+1)]; double leg[gsl_sf_legendre_array_n(n+nu+1)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
complex double bes[n+nu+2]; complex double bes[n+nu+2];
if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort(); if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort();
complex double sum = 0; complex double sum = 0;
for (int q = 0; q <= Qmax; ++q) { for (int q = 0; q <= Qmax; ++q) {
@ -338,10 +339,6 @@ void qpms_trans_calculator_free(qpms_trans_calculator *c) {
free(c); free(c);
} }
static inline size_t qpms_mn2y(int m, int n) {
return (size_t) n * (n + 1) + m - 1;
}
static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c, static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c,
int m, int n, int mu, int nu){ int m, int n, int mu, int nu){
return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu); return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);
@ -470,9 +467,9 @@ static void qpms_trans_calculator_multipliers_B_Taylor(
} }
} }
int qpms_trans_calculator_multipliers_A(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
switch (norm) { switch (norm) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax); qpms_trans_calculator_multipliers_A_Taylor(dest,m,n,mu,nu,qmax);
return 0; return 0;
break; break;
@ -482,9 +479,9 @@ int qpms_trans_calculator_multipliers_A(qpms_normalization_t norm, complex doubl
assert(0); assert(0);
} }
int qpms_trans_calculator_multipliers_B(qpms_normalization_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) {
switch (norm) { switch (norm) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax); qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax);
return 0; return 0;
break; break;
@ -495,14 +492,14 @@ int qpms_trans_calculator_multipliers_B(qpms_normalization_t norm, complex doubl
} }
qpms_trans_calculator qpms_trans_calculator
*qpms_trans_calculator_init (int lMax, qpms_normalization_t normalization) { *qpms_trans_calculator_init (int lMax, qpms_normalisation_t normalisation) {
assert(lMax > 0); assert(lMax > 0);
qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator)); qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator));
c->lMax = lMax; c->lMax = lMax;
c->nelem = lMax * (lMax+2); c->nelem = lMax * (lMax+2);
c->A_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *)); c->A_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
c->B_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *)); c->B_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
c->normalization = normalization; c->normalisation = normalisation;
size_t *qmaxes = malloc(SQ(c->nelem) * sizeof(size_t)); size_t *qmaxes = malloc(SQ(c->nelem) * sizeof(size_t));
size_t qmaxsum = 0; size_t qmaxsum = 0;
for(size_t y = 0; y < c->nelem; y++) for(size_t y = 0; y < c->nelem; y++)
@ -525,7 +522,7 @@ qpms_trans_calculator
int m, n, mu, nu; int m, n, mu, nu;
qpms_y2mn_p(y, &m, &n); qpms_y2mn_p(y, &m, &n);
qpms_y2mn_p(yu, &mu, &nu); qpms_y2mn_p(yu, &mu, &nu);
qpms_trans_calculator_multipliers_A(normalization, qpms_trans_calculator_multipliers_A(normalisation,
c->A_multipliers[i], m, n, mu, nu, qmaxes[i]); c->A_multipliers[i], m, n, mu, nu, qmaxes[i]);
} }
@ -550,7 +547,7 @@ qpms_trans_calculator
int m, n, mu, nu; int m, n, mu, nu;
qpms_y2mn_p(y, &m, &n); qpms_y2mn_p(y, &m, &n);
qpms_y2mn_p(yu, &mu, &nu); qpms_y2mn_p(yu, &mu, &nu);
qpms_trans_calculator_multipliers_B(normalization, qpms_trans_calculator_multipliers_B(normalisation,
c->B_multipliers[i], m, n, mu, nu, qmaxes[i]); c->B_multipliers[i], m, n, mu, nu, qmaxes[i]);
} }
@ -586,13 +583,13 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c,
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR)
// TODO warn? // TODO warn?
return NAN+I*NAN; return NAN+I*NAN;
switch(c->normalization) { switch(c->normalisation) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
{ {
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu, if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf)) abort();
if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bessel_buf)) abort(); if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort();
return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
} }
@ -631,13 +628,13 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c,
if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR)
// TODO warn? // TODO warn?
return NAN+I*NAN; return NAN+I*NAN;
switch(c->normalization) { switch(c->normalisation) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
{ {
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf)) abort();
if (qpms_sph_bessel_array(J, n+nu+2, kdlj.r, bessel_buf)) abort(); if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort();
return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
} }
@ -660,13 +657,13 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
// TODO warn? different return value? // TODO warn? different return value?
return 0; return 0;
} }
switch(c->normalization) { switch(c->normalisation) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
{ {
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf)) abort();
if (qpms_sph_bessel_array(J, n+nu+2, kdlj.r, bessel_buf)) abort(); if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort();
*Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
*Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
@ -698,13 +695,13 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c,
// TODO warn? different return value? // TODO warn? different return value?
return 0; return 0;
} }
switch(c->normalization) { switch(c->normalisation) {
case QPMS_NORMALIZATION_TAYLOR: case QPMS_NORMALISATION_TAYLOR:
{ {
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf)) abort();
if (qpms_sph_bessel_array(J, 2*c->lMax+2, kdlj.r, bessel_buf)) abort(); if (qpms_sph_bessel_fill(J, 2*c->lMax+2, kdlj.r, bessel_buf)) abort();
size_t desti = 0, srci = 0; size_t desti = 0, srci = 0;
for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) { for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) { for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {

7
qpms/translations.h
View File

@ -9,7 +9,8 @@
double qpms_legendre0(int m, int n); double qpms_legendre0(int m, int n);
// Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2) // Associated Legendre polynomial derivative at zero argument (DLMF 14.5.2)
double qpms_legendred0(int m, int n); double qpms_legendred0(int m, int n);
int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double *result_array); // TODO unify types
int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array);
// TODO replace the xplicit "Taylor" functions with general, // TODO replace the xplicit "Taylor" functions with general,
// taking qpms_bessel_t argument. // taking qpms_bessel_t argument.
@ -30,10 +31,10 @@ typedef struct qpms_trans_calculator {
size_t nelem; size_t nelem;
complex double **A_multipliers; complex double **A_multipliers;
complex double **B_multipliers; complex double **B_multipliers;
qpms_normalization_t normalization; qpms_normalisation_t normalisation;
} qpms_trans_calculator; } qpms_trans_calculator;
qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalization_t nt); qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalisation_t nt);
void qpms_trans_calculator_free(qpms_trans_calculator *); void qpms_trans_calculator_free(qpms_trans_calculator *);
complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c,

29
qpms/vectors.h
View File

@ -6,11 +6,11 @@
//static inline double vectors_h_sq(double x) {return x*x;} //static inline double vectors_h_sq(double x) {return x*x;}
static inline double cart3norm(cart3_t v) { static inline double cart3norm(const cart3_t v) {
return sqrt(v.x*v.x + v.y*v.y + v.z*v.z); return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
} }
static inline sph_t cart2sph(cart3_t cart) { static inline sph_t cart2sph(const cart3_t cart) {
sph_t sph; sph_t sph;
sph.r = cart3norm(cart); sph.r = cart3norm(cart);
sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2; sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
@ -18,7 +18,7 @@ static inline sph_t cart2sph(cart3_t cart) {
return sph; return sph;
} }
static inline cart3_t sph2cart(sph_t sph) { static inline cart3_t sph2cart(const sph_t sph) {
cart3_t cart; cart3_t cart;
double sin_th = sin(sph.theta); double sin_th = sin(sph.theta);
cart.x = sph.r * sin_th * cos(sph.phi); cart.x = sph.r * sin_th * cos(sph.phi);
@ -27,5 +27,28 @@ static inline cart3_t sph2cart(sph_t sph) {
return cart; return cart;
} }
// equivalent to sph_loccart2cart in qpms_p.py
static inline ccart3_t csphvec2cart(const csphvec_t sphvec, const sph_t at) {
ccart3_t res = {0, 0, 0};
const double st = sin(at.theta);
const double ct = cos(at.theta);
const double sf = sin(at.phi);
const double cf = cos(at.phi);
const double rx = st * cf;
const double ry = st * sf;
const double rz = ct;
const double tx = ct * cf;
const double ty = ct * sf;
const double tz = -st;
const double fx = -sf;
const double fy = cf;
const double fz = 0.;
res.x = rx * sphvec.rc + tx * sphvec.thetac + fx * sphvec.phic;
res.y = ry * sphvec.rc + ty * sphvec.thetac + fy * sphvec.phic;
res.z = rz * sphvec.rc + tz * sphvec.thetac + fz * sphvec.phic;
return sphvec;
}
#endif //VECTORS_H #endif //VECTORS_H

92
qpms/vswf.c
View File

@ -1,10 +1,12 @@
#include "vswf.h" #include "vswf.h"
#include "indexing.h" #include "indexing.h"
#include "translations.h" // TODO move qpms_sph_bessel_fill elsewhere
#include <math.h> #include <math.h>
#include <gsl/gsl_math.h> #include <gsl/gsl_math.h>
#include <gsl/gsl_sf_legendre.h> #include <gsl/gsl_sf_legendre.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "assert_cython_workaround.h"
// Legendre functions also for negative m, see DLMF 14.9.3 // 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, qpms_errno_t qpms_legendre_deriv_y_fill(double *target, double *target_deriv, double x, qpms_l_t lMax,
@ -156,3 +158,93 @@ void qpms_pitau_free(qpms_pitau_t x) {
free(x.pi); free(x.pi);
free(x.tau); free(x.tau);
} }
csphvec_t qpms_vswf_single_el(qpms_m_t m, qpms_l_t l, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm) {
lmcheck(l,m);
csphvec_t N;
complex double *bessel = malloc((l+1)*sizeof(complex double));
if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
complex double eimf = cexp(m * kdlj.phi * I);
qpms_y_t y = qpms_mn2y(m,l);
N.rc = l*(l+1) * pt.leg[y] * bessel[l] / kdlj.r * eimf;
complex double besselfac = bessel[l-1] - l * bessel[l] / kdlj.r;
N.thetac = pt.tau[y] * besselfac * eimf;
N.phic = pt.pi[y] * besselfac * I * eimf;
qpms_pitau_free(pt);
free(bessel);
return N;
}
csphvec_t qpms_vswf_single_mg(qpms_m_t m, qpms_l_t l, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm) {
lmcheck(l,m);
csphvec_t M;
complex double *bessel = malloc((l+1)*sizeof(complex double));
if(qpms_sph_bessel_fill(btyp, l, kdlj.r, bessel)) abort();
qpms_pitau_t pt = qpms_pitau_get(kdlj.theta, l, norm);
complex double eimf = cexp(m * kdlj.phi * I);
qpms_y_t y = qpms_mn2y(m,l);
M.rc = 0.;
M.thetac = pt.pi[y] * bessel[l] * I * eimf;
M.phic = -pt.tau[y] * bessel[l] * eimf;
qpms_pitau_free(pt);
free(bessel);
return M;
}
qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm) {
qpms_vswfset_sph_t *res = malloc(sizeof(qpms_vswfset_sph_t));
res->lMax = lMax;
qpms_y_t nelem = qpms_lMax2nelem(lMax);
res->el = malloc(sizeof(csphvec_t)*nelem);
res->mg = malloc(sizeof(csphvec_t)*nelem);
if(QPMS_SUCCESS != qpms_vswf_fill(res->mg, res->el, lMax, kdlj, btyp, norm))
abort(); // or return NULL? or rather assert?
return res;
}
void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *w) {
assert(NULL != w && NULL != w->el && NULL != w->mg);
free(w->el);
free(w->mg);
free(w);
}
qpms_errno_t qpms_vswf_fill(csphvec_t *mgtarget, csphvec_t *eltarget,
qpms_l_t lMax, sph_t kr,
qpms_bessel_t btyp, qpms_normalisation_t norm) {
assert(lMax >= 1);
complex double *bessel = malloc((lMax+1)*sizeof(complex double));
if(qpms_sph_bessel_fill(btyp, lMax, kr.r, bessel)) abort();
qpms_pitau_t pt = qpms_pitau_get(kr.theta, lMax, norm);
complex double const *pbes = bessel + 1; // starting from l = 1
double const *pleg = pt.leg;
double const *ppi = pt.pi;
double const *ptau = pt.tau;
csphvec_t *pmg = mgtarget, *pel = eltarget;
for(qpms_l_t l = 1; l <= lMax; ++l) {
complex double besfac = *pbes / kr.r;
complex double besderfac = *(pbes-1) - l * besfac;
for(qpms_m_t m = -l; m <= l; ++m) {
complex double eimf = cexp(m * kr.phi * I);
pel->rc = l*(l+1) * (*pleg) * besfac * eimf;
pel->thetac = *ptau * besderfac * eimf;
pel->phic = *ppi * besderfac * I * eimf;
pmg->rc = 0.;
pmg->thetac = *ppi * (*pbes) * I * eimf;
pmg->phic = - *ptau * (*pbes) * eimf;
++pleg; ++ppi; ++ptau; ++pel; ++pmg;
}
++pbes;
}
free(bessel);
qpms_pitau_free(pt);
return QPMS_SUCCESS;
}

3
qpms/vswf.h
View File

@ -31,6 +31,9 @@ qpms_errno_t qpms_legendre_deriv_y_get(double **result, double **result_deriv, d
qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x, qpms_errno_t qpms_legendre_deriv_y_fill(double *where, double *where_deriv, double x,
qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase);
qpms_errno_t qpms_vswf_fill(csphvec_t *resultM, csphvec_t *resultN, qpms_l_t lMax, sph_t kdrj,
qpms_bessel_t btyp, qpms_normalisation_t norm);
qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj, qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
qpms_bessel_t btyp, qpms_normalisation_t norm);//NI qpms_bessel_t btyp, qpms_normalisation_t norm);//NI
void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI void qpms_vswfset_sph_pfree(qpms_vswfset_sph_t *);//NI