Continue writing translation coeff calculator object, etc.

- Removed unnecessary continue in qpms_trans_single_B_Taylor (it
  was there before previously incorrect calculation of q_max),
- for the same reason, the functions no longer crash for
  certain indices, so the old tests now pass without artificial
  constraints on the indices.


Former-commit-id: 0047fe18596ea006e7a2defb192787572b2a7fc2

qpms/test_translations.c

@@ -16,17 +16,17 @@ testcase_single_trans_t testcases_Taylor[] = {
 
 int main() {
 	for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) {
-		if (tc->n > 12 || tc->nu > 12 || !tc->n || !tc->nu ) continue;
+		//if (tc->n > 40 || tc->nu > 40 ) continue;
 
 		printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu);
 		complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J);
 		complex double B = qpms_trans_single_B_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J);
 		printf("A = %.16f+%.16fj, relerr=%.16f, J=%d\n",
-			creal(A), cimag(A),
+			creal(A), cimag(A),  (0 == cabs(tc->result_A - A)) ? 0 :
 			cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)),
 			tc->J);
 		printf("B = %.16f+%.16fj, relerr=%.16f, J=%d\n",
-			creal(B), cimag(B),
+			creal(B), cimag(B), (0 == cabs(tc->result_B - B)) ? 0 :
 			cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)),
 			tc->J);
 	}

qpms/translations.c

@@ -137,7 +137,8 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd
 	    double a3q_n = a3q[q]/a3q0;
 	    complex double zp_ = bes[p+1];
 	    int Pp_order_ = mu-m;
-	    if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
+	    //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
+	    assert(p+1 >= abs(Pp_order_));
 	    double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))];
 	    if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
 	    complex double summandq = ((2*(n+1)*(nu-mu)*a2q_n
@@ -172,30 +173,82 @@ complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu,
 	return qpms_trans_single_B_Taylor(m,n,mu,nu,kdlj,r_ge_d,J);
 }
 
-#if 0
+void qpms_trans_calculator_free(qpms_trans_calculator *c) {
-typedef struct qpms_trans_calculator {
+	free(c->A_multipliers[0]);
-	int lMax;
+	free(c->A_multipliers);
-	size_t nelem;
+	free(c->B_multipliers[0]);
-	double **A_coeffs;
+	free(c->B_multipliers);
-	double **B_coeffs;
+	free(c);
-	// TODO enum normalization
+}
-} qpms_trans_calculator;
 	
 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__indexcalc(const qpms_trans_calculator *c,
+static inline int qpms_y2n(size_t y) {
+	return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
+}
+
+static inline int qpms_yn2m(size_t y, int n) {
+	return y-qpms_mn2y(0,n);
+}
+
+static inline void qpms_y2mn_p(size_t y, int *m, int *n){
+	*m=qpms_yn2m(y,*n=qpms_y2n(y));
+}
+
+static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c,
 		int m, int n, int mu, int nu){
 	return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu);
 }
 
-struct qpms_trans_calculator
+static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator *c,
-*qpms_trans_calculator_init_Taylor (int lMax) {
+		size_t y, size_t yu) {
+	return c->nelem * y + yu;
+}
+
+#define SQ(x) ((x)*(x))
+
+qpms_trans_calculator
+*qpms_trans_calculator_init (int lMax, qpms_normalization_t normalization) {
 	qpms_trans_calculator *c = malloc(sizeof(qpms_trans_calculator));
 	c->lMax = lMax;
 	c->nelem = lMax * (lMax+2);
-	c->A_coeffs = malloc((1+c->nelem) * sizeof(double *));
+	c->A_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
-	c->B_coeffs = malloc((1+c->nelem) * sizeof(double *));
+	c->B_multipliers = malloc((1+SQ(c->nelem)) * sizeof(complex double *));
+	c->normalization = normalization;
+	size_t *qmaxes = malloc(SQ(c->nelem) * sizeof(size_t));
+	size_t qmaxsum = 0;
+	for(size_t y = 0; y < c->nelem; y++)
+		for(size_t yu = 0; yu < c->nelem; yu++) {
+			int m,n, mu, nu;
+			qpms_y2mn_p(y,&m,&n);
+			qpms_y2mn_p(yu,&mu,&nu);
+			qmaxsum += (qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] = gaunt_q_max(-m,n,mu,nu));
+		}
+	c->A_multipliers[0] = malloc(qmaxsum * sizeof(complex double));
+	for(size_t i = 0; i < SQ(c->nelem); ++i)
+		c->A_multipliers[i+1] = c->A_multipliers[i] + qmaxes[i];
+	// TODO here comes the filling of A_multipliers
+	
+
+
+	qmaxsum = 0;
+	for(size_t y=0; y < c->nelem; y++)
+		for(size_t yu = 0; yu < c->nelem; yu++) {
+			int m, n, mu, nu;
+			qpms_y2mn_p(y,&m,&n);
+			qpms_y2mn_p(yu,&mu,&nu);
+			qmaxsum += (qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] = gaunt_q_max(-m,n+1,mu,nu));
+		}
+	c->B_multipliers[0] = malloc(qmaxsum * sizeof(complex double));
+	for(size_t i = 0; i < SQ(c->nelem); ++i)
+		c->B_multipliers[i+1] = c->B_multipliers[i] + qmaxes[i];
+
+	// TODO here comes the filling of B_multipliers
+
+	free(qmaxes);
+	return c;
+}
+
 
-#endif

qpms/translations.h

@@ -3,10 +3,15 @@
 #include "vectors.h"
 #include <complex.h>
 #include <stdbool.h>
-
+#include <stddef.h>
 double qpms_legendre0(int m, int n); 
 double qpms_legendred0(int m, int n); 
 
+typedef enum {
+	QPMS_NORMALIZATION_TAYLOR = 1,
+	QPMS_NORMALIZATION_UNDEF = 0
+} qpms_normalization_t;
+
 typedef enum {
         QPMS_BESSEL_REGULAR = 1, // regular function j
         QPMS_BESSEL_SINGULAR = 2, // singular function y
@@ -29,4 +34,16 @@ complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, doub
 complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r,
 		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
 
+
+typedef struct qpms_trans_calculator {
+        int lMax;
+        size_t nelem;
+        complex double **A_multipliers;
+        complex double **B_multipliers;
+        qpms_normalization_t normalization;
+} qpms_trans_calculator;
+
+qpms_trans_calculator *qpms_trans_calculator_init(int lMax, qpms_normalization_t nt);
+void qpms_trans_calculator_free(qpms_trans_calculator *);
+
 #endif // QPMS_TRANSLATIONS_H