Koefficient A a testy. Kompiluje, většinou funguje, někdy nan

Former-commit-id: 3303493088dc8edca027fb93d49d62166e71544e

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
 	

qpms/gaunt.h

@@ -2,12 +2,12 @@
 #define GAUNT_H
 #include <stdlib.h>
 
-#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

qpms/gsltest.c

@@ -0,0 +1,54 @@
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <gsl/gsl_sf_legendre.h>
+
+
+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;
+}
+
+
+

qpms/test_translations.c

@@ -0,0 +1,25 @@
+#include "translations.h"
+#include <stdio.h>
+//#include <math.h>
+#include <complex.h>
+
+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)));
+	}
+}
+
+

qpms/testcases_taylor.REMOVED.git-id

@@ -0,0 +1 @@
+2b61a18b6426306ba80a50ee30339651313fc44f

qpms/translations.c

@@ -1,17 +1,105 @@
 #include <math.h>
 #include "gaunt.h"
+#include "translations.h"
+#include <stdbool.h>
+#include <gsl/gsl_sf_legendre.h>
+#include <gsl/gsl_sf_bessel.h>
+#include <assert.h>
+
 
 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;
+}
+

qpms/translations.h

@@ -0,0 +1,24 @@
+#ifndef QPMS_TRANSLATIONS_H
+#define QPMS_TRANSLATIONS_H
+#include "vectors.h"
+#include <complex.h>
+#include <stdbool.h>
+
+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

qpms/vectors.h

@@ -1,6 +1,7 @@
 #ifndef VECTORS_H
 #define VECTORS_H
 #include <math.h>
+#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);