From 34940fc7e92d8d9a03968226f29a82fb6f6714f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 12 Apr 2017 14:17:07 +0300 Subject: [PATCH] =?UTF-8?q?Koefficient=20A=20a=20testy.=20Kompiluje,=20v?= =?UTF-8?q?=C4=9Bt=C5=A1inou=20funguje,=20n=C4=9Bkdy=20nan?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: 3303493088dc8edca027fb93d49d62166e71544e --- qpms/gaunt.c | 2 + qpms/gaunt.h | 4 +- qpms/gsltest.c | 54 ++++++++++++++++ qpms/test_translations.c | 25 ++++++++ qpms/testcases_taylor.REMOVED.git-id | 1 + qpms/translations.c | 96 ++++++++++++++++++++++++++-- qpms/translations.h | 24 +++++++ qpms/vectors.h | 7 +- 8 files changed, 204 insertions(+), 9 deletions(-) create mode 100644 qpms/gsltest.c create mode 100644 qpms/test_translations.c create mode 100644 qpms/testcases_taylor.REMOVED.git-id create mode 100644 qpms/translations.h diff --git a/qpms/gaunt.c b/qpms/gaunt.c index 4652b49..013f39b 100644 --- a/qpms/gaunt.c +++ b/qpms/gaunt.c @@ -1173,6 +1173,7 @@ static inline int q_max(int m, int n, int mu, int nu) { return MIN(n, MIN(nu,(n+nu-abs(m+mu))/2)); } +/* THIS THING IS TOTALLY WRONG int gaunt(int m, int n, int mu, int nu, double *v) { int err = 0; int qmax = q_max(m,n,mu,nu); @@ -1181,6 +1182,7 @@ int gaunt(int m, int n, int mu, int nu, double *v) { gaunt_xu(m, n, mu, nu, qmax, v, &err); return err; } +*/ #ifdef GAUNTTEST diff --git a/qpms/gaunt.h b/qpms/gaunt.h index 8536075..a28dafa 100644 --- a/qpms/gaunt.h +++ b/qpms/gaunt.h @@ -2,12 +2,12 @@ #define GAUNT_H #include -#define _GAUNT_H_MIN(x,y) (((x) > (y) ? (y) : (x)) +#define _GAUNT_H_MIN(x,y) (((x) > (y)) ? (y) : (x)) static inline int gaunt_q_max(int m, int n, int mu, int nu) { return _GAUNT_H_MIN(n, _GAUNT_H_MIN(nu, n+nu-abs(m+mu))); } #undef _GAUNT_H_MIN void gaunt_xu(int m, int n, int mu, int nu, int qmax, double *v_aq, int *err); -int gaunt(int m, int n, int mu, int nu, double *v_aq); +//int gaunt(int m, int n, int mu, int nu, double *v_aq); #endif //GAUNT_H diff --git a/qpms/gsltest.c b/qpms/gsltest.c new file mode 100644 index 0000000..c6764e8 --- /dev/null +++ b/qpms/gsltest.c @@ -0,0 +1,54 @@ +#include +#include +#include +#include + + +int main(int argc, char **argv) { + int lmax; + int scanned; + double x; + while (EOF != (scanned = scanf("%d %lf", &lmax, &x))) { + if (scanned != 2) continue; + size_t as = gsl_sf_legendre_array_n(lmax); + double *r_none = calloc(as, sizeof(double)); + double *d_none = calloc(as, sizeof(double)); + double *r_schmidt = calloc(as, sizeof(double)); + double *d_schmidt = calloc(as, sizeof(double)); + double *r_spharm = calloc(as, sizeof(double)); + double *d_spharm = calloc(as, sizeof(double)); + double *r_full = calloc(as, sizeof(double)); + double *d_full = calloc(as, sizeof(double)); + if( 0 + || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_NONE,lmax,x,-1,r_none,d_none) + || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_SCHMIDT,lmax,x,-1,r_schmidt,d_schmidt) + || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_SPHARM,lmax,x,-1,r_spharm,d_spharm) + || gsl_sf_legendre_deriv_alt_array_e(GSL_SF_LEGENDRE_FULL,lmax,x,-1,r_full,d_full) + ) fprintf(stderr, "Something gone wrong for lmax = %d, x = %.15e!\n", lmax, x); + for (int l = 0; l <= lmax; ++l) + for (int m = 0; m <= l; ++m){ + size_t i = gsl_sf_legendre_array_index(l,m); + printf("P(%d,%d)\t%.16e\t%.16e\t%.16e\t%.16e\n", l, m, + r_none[i], r_schmidt[i], r_spharm[i], r_full[i]); + printf("dP(%d,%d)\t%.16e\t%.16e\t%.16e\t%.16e\n", l, m, + d_none[i], d_schmidt[i], d_spharm[i], d_full[i]); + + } + + + + + free(r_none); + free(d_none); + free(r_schmidt); + free(d_schmidt); + free(r_spharm); + free(d_spharm); + free(r_full); + free(d_full); + } + return 0; +} + + + diff --git a/qpms/test_translations.c b/qpms/test_translations.c new file mode 100644 index 0000000..a6210c3 --- /dev/null +++ b/qpms/test_translations.c @@ -0,0 +1,25 @@ +#include "translations.h" +#include +//#include +#include + +typedef struct { + int m, n, mu, nu; + sph_t kdlj; + qpms_bessel_t J; + complex double result_A, result_B; +} testcase_single_trans_t; + +testcase_single_trans_t testcases_Taylor[] = { +#include "testcases_taylor" +}; + +int main() { + for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) { + complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J); + printf("m=%d, n=%d, mu=%d, nu=%d, relerr=%.16f\n", tc->m,tc->n,tc->mu,tc->nu, + cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A))); + } +} + + diff --git a/qpms/testcases_taylor.REMOVED.git-id b/qpms/testcases_taylor.REMOVED.git-id new file mode 100644 index 0000000..ba0ba20 --- /dev/null +++ b/qpms/testcases_taylor.REMOVED.git-id @@ -0,0 +1 @@ +2b61a18b6426306ba80a50ee30339651313fc44f \ No newline at end of file diff --git a/qpms/translations.c b/qpms/translations.c index bb170de..ec06610 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1,17 +1,105 @@ #include #include "gaunt.h" +#include "translations.h" +#include +#include +#include +#include + static const double sqrtpi = 1.7724538509055160272981674833411451827975494561223871; //static const double ln2 = 0.693147180559945309417232121458176568075500134360255254120; // Associated Legendre polynomial at zero argument (DLMF 14.5.1) -double legendre0(int m, int n) { - return pow(2,m) * sqrtpi / gamma(.5*n - .5*m + .5) / gamma(.5*n-.5*m); +double qpms_legendre0(int m, int n) { + return pow(2,m) * sqrtpi / tgamma(.5*n - .5*m + .5) / tgamma(.5*n-.5*m); +} + +static inline int min1pow(int x) { + return (x % 2) ? -1 : 1; +} + +static inline complex double ipow(int x) { + return cpow(I, x); } // Derivative of associated Legendre polynomial at zero argument (DLMF 14.5.2) -double legendreD0(int m, int n) { - return -2 * legendre0(m, n); +double qpms_legendreD0(int m, int 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 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); } +complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J) { // TODO make J enum + if(r_ge_d) J = QPMS_BESSEL_REGULAR; + double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) + +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) + +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)); + double costheta = cos(kdlj.theta); + int qmax = gaunt_q_max(-m,n,mu,nu); // nemá tu být +m? + // N.B. -m !!!!!! + double a1q[qmax+1]; + int err; + gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); + double a1q0 = a1q[0]; + if (err) abort(); + //double *leg = malloc(sizeof(double)*gsl_sf_legendre_array_n(n+nu)); + //if (!leg) abort(); + 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(); + complex double bes[n+nu+1]; + if (qpms_sph_bessel_array(J, n+nu, kdlj.r, bes)) abort(); + complex double sum = 0; + for(int q = 0; q <= qmax; ++q) { + int p = n+nu-2*q; + double a1q_n = a1q[q] / a1q0; + int Pp_order = mu-m; + double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))]; + if (Pp_order < 0) Pp *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); + complex double zp = bes[p]; + complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp; + sum += summandq; + } + //free(leg); + //free(bes); + complex double presum = exp(exponent); + presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n); + + complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( + lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); + return (presum / prenormratio) * sum; +} + diff --git a/qpms/translations.h b/qpms/translations.h new file mode 100644 index 0000000..bc673ea --- /dev/null +++ b/qpms/translations.h @@ -0,0 +1,24 @@ +#ifndef QPMS_TRANSLATIONS_H +#define QPMS_TRANSLATIONS_H +#include "vectors.h" +#include +#include + +double qpms_legendre0(int m, int n); +double qpms_legendred0(int m, int n); + +typedef enum { + QPMS_BESSEL_REGULAR = 1, // regular function j + QPMS_BESSEL_SINGULAR = 2, // singular function y + QPMS_HANKEL_PLUS = 3, // hankel function h1 = j + I*y + QPMS_HANKEL_MINUS = 4, // hankel function h2 = j - I*y + QPMS_BESSEL_UNDEF = 0 +} qpms_bessel_t; + +int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double *result_array); + +complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj, + bool r_ge_d, qpms_bessel_t J); + + +#endif // QPMS_TRANSLATIONS_H diff --git a/qpms/vectors.h b/qpms/vectors.h index eafbfaa..3230223 100644 --- a/qpms/vectors.h +++ b/qpms/vectors.h @@ -1,6 +1,7 @@ #ifndef VECTORS_H #define VECTORS_H #include +#define M_PI_2 (1.570796326794896619231321691639751442098584699687552910487) typedef struct { double x, y, z; @@ -21,7 +22,7 @@ typedef struct { //static inline double vectors_h_sq(double x) {return x*x;} -static inline cart3norm(cart3_t v) { +static inline double cart3norm(cart3_t v) { return sqrt(v.x*v.x + v.y*v.y + v.z*v.z); } @@ -33,8 +34,8 @@ static inline sph_t cart2sph(cart3_t cart) { return sph; } -static inline cart_t sph2cart(sph_t sph) { - cart_t cart; +static inline cart3_t sph2cart(sph_t sph) { + cart3_t cart; double sin_th = sin(sph.theta); cart.x = sph.r * sin_th * cos(sph.phi); cart.y = sph.r * sin_th * sin(sph.phi);