Plane wave decomposition compiles but yields wrong results.

Former-commit-id: 0cfe82f1efe0e68d41a2f4bf315dfac4618cf435
This commit is contained in:
Marek Nečada 2018-02-02 14:23:15 +00:00
parent 7988b4e942
commit f458e94995
4 changed files with 20 additions and 4 deletions

View File

@ -1,7 +1,7 @@
#include <assert.h>
#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;

View File

@ -7,6 +7,9 @@
#ifndef M_PI_2
#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487)
#endif
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
// integer index types
typedef int qpms_lm_t;

View File

@ -34,6 +34,16 @@ static inline cart3_t cart3_add(const cart3_t a, const cart3_t b) {
return res;
}
static inline cart3_t cart3_scale(const double c, const cart3_t v) {
cart3_t res = {c * v.x, c * v.y, c * v.z};
return res;
}
static inline ccart3_t ccart3_scale(const complex double c, const ccart3_t v) {
ccart3_t res = {c * v.x, c * v.y, c * v.z};
return res;
}
static inline csphvec_t csphvec_add(const csphvec_t a, const csphvec_t b) {
csphvec_t res = {a.rc + b.rc, a.thetac + b.thetac, a.phic + b.phic};
return res;
@ -44,13 +54,17 @@ static inline csphvec_t csphvec_scale(complex double c, const csphvec_t v) {
return res;
}
static inline complex double csphvec_dotnc(csphvec_t a, csphvec_t b) {
static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b) {
//N.B. no complex conjugation done here
return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
}
static inline double cart3_dot(const cart3_t a, const cart3_t b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
// equivalent to sph_loccart2cart in qpms_p.py
static inline ccart3_t csphvec2cart(const csphvec_t sphvec, const sph_t at) {
static inline ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) {
const double st = sin(at.theta);
const double ct = cos(at.theta);
const double sf = sin(at.phi);

View File

@ -400,7 +400,6 @@ static inline complex double ipowl(qpms_l_t l) {
qpms_errno_t qpms_planewave2vswf_fill_sph(sph_t wavedir, csphvec_t amplitude,
complex double *target_longcoeff, complex double *target_mgcoeff,
complex double *target_elcoeff, qpms_l_t lMax, qpms_normalisation_t norm) {
abort(); //NI
qpms_y_t nelem = qpms_lMax2nelem(lMax);
csphvec_t * const dual_A1 = malloc(3*nelem*sizeof(csphvec_t)), *const dual_A2 = dual_A1 + nelem,
* const dual_A3 = dual_A2 + nelem;