Rozdělané Legendrovy polynomy por záporná m (jdu spát)

Former-commit-id: c0a0eba0bdc993de7199e57cc9f8a7c9964667a7

qpms/qpms_types.h

@@ -11,6 +11,11 @@ typedef unsigned int qpms_l_t;
 typedef qpms_lm_t qpms_m_t;
 typedef size_t qpms_y_t;
 
+typedef enum {
+	QPMS_SUCCESS = 0;
+	QPMS_ERROR = 1;
+} qpms_errno_t;
+
 // Normalisations
 typedef enum {
 	QPMS_NORMALIZATION_XU = 3, // NI!

qpms/vswf.c

@@ -0,0 +1,52 @@
+#include "vswf.h"
+
+// 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,
+		gsl_sf_legendre_t lnorm, double csphase)
+{
+	size_t n = gsl_sf_legenre_array_n(lMax);
+	double *legendre_tmp = malloc(n * sizeof(double));
+	double *legendre_deriv_tmp = malloc(n * sizeof(double));
+	int gsl_errno = gsl_sf_legendre_deriv_array_e(
+			lnorm, (size_t)lMax, x, csphase, legendre_tmp,legendre_tmp_deriv);
+	for (qpms_l_t l = 0; l <= lMax; ++l)
+			for (qpms_m_t m = 0; m <= l; ++m) {
+				qpms_y_t y = gpms_mn2y(m,l);
+				size_t i = gsl_sf_legenre_array_index(l,m);
+				target[y] = legendre_tmp[i];
+				target_deriv[y] = legendre_deriv_tmp[i];
+			}
+	switch(lnorm) {
+		case GSL_SF_LEGEDRE_NONE:
+			for (qpms_l_t l = 0; l <= lMax; ++l)
+				for (qpms_m_t m = 1; m <= l; ++m) {
+				qpms_y_t y = gpms_mn2y(-m,l);
+				size_t i = gsl_sf_legenre_array_index(l,m);
+				// viz DLMF 14.9.3, čert ví, jak je to s cs fasí.
+				double factor = exp(lgamma(l-m+1)-lgamma(n+m+1))*((m%2)?-1:1);
+				target[y] = factor * legendre_tmp[i];
+				target_deriv[y] = factor * legendre_deriv_tmp;
+				}
+			break;
+		case GSL_SF_LEGENDRE_SCHMIDT:
+		case GSL_SF_LEGENDRE_SPHARM:
+		case GSL_SF_LEGENDRE_FULL:
+		default:
+			abort(); //NI
+			break;
+	}
+	free(legendre_tmp);
+	free(legendre_deriv_tmp);
+	return QPMS_SUCCESS;
+}
+
+// FIXME ZAPOMNĚL JSEM NA POLE DERIVACÍ (TÉŽ HLAVIČKOVÝ SOUBOR)
+double *qpms_legendre_deriv_y_get(double x, qpms_l_t lMax, gsle_sf_legendre_t lnorm,
+		double csphase)
+{
+	double *ar = malloc(sizeof(double)*qpms_lMaxnelem(lMax));
+	if (qpms_legendre_deriv_y_fill(ar, x, lMax, lnorm, csphase))
+		abort();
+}
+
+

qpms/vswf.h

@@ -1,6 +1,7 @@
 #ifndef QPMS_VSWF_H
 #define QPMS_VSWF_H
 #include "qpms_types.h"
+#include <gsl/gsl_sf_legendre.h>
 
 // Electric wave N; NI
 complex double qpms_vswf_single_el(int m, int n, sph_t kdlj,
@@ -18,10 +19,21 @@ typedef struct {
 	csphvec_t *el, *mg;
 } qpms_vswfset_sph_t;
 
+
+double *qpms_legenre_deriv_y_get(double x, qpms_l_t lMax, 
+		gsle_sf_legendre_t lnorm, double csphase);
+qpms_errno_t qpms_legenre_deriv_y_fill(double *where, double x, 
+		qpms_l_t lMax, gsl_sf_legendre_t lnorm, double csphase); //NI
+
 qpms_vswfset_sph_t *qpms_vswfset_make(qpms_l_t lMax, sph_t kdlj,
 	qpms_bessel_t btyp, qpms_normalization_t norm);//NI
 void qpms_vswfst_sph_pfree(qpms_vswfset_t *);//NI
 
+double *qpms_legendre_y_get(double x, qpms_l_t lMax, qpms_normalisation_t norm);//NI
+double *qpms_legendre0d_y_get(qpms_l_t lMax, qpms_normalization_t norm); //NI
+double *qpms_legendre_plus1d_y_get(qpms_l_t lMax, qpms_normalization_t norm); //NI
+double *qpms_legendre_minus1d_y_get(qpms_l_t lMax, qpms_normalization_t norm); //NI
+
 
 
 // array of pi, tau auxillary function (see [1,(37)])

qpms/vswf_xu.c

@@ -1,4 +1,6 @@
 #include "translations.h"
+#include "qpms_types.h"
+#include "indexing.h"
 #include <math.h>
 #include <stdlib.h> //abort();
 
@@ -10,6 +12,7 @@
 /* 
  * The following value delimits the interval for which the computation
  * of the pi and tau functions [1,(37)] actually takes place.
+ * see also DLMF §14.8
  * For x in [-1, -1 + PITAU_0THRESHOLD] the x->-1 limit is taken,
  * for x in [-1 + PITAU_0THRESHOLD, 1 - PITAU_0THRESHOLD] value for x is calculated,
  * for x in [1 - PITAU_THRESHOLD, 1] the x->1 limit is taken.
@@ -18,28 +21,61 @@
  */
 #define PITAU_0THRESHOLD 1e-7
 
-int taumncos_Xu_zerolim_fill(int lmax, double *where);
-double * taumncos_Xu_zerolim_get(int lmax);
-int taumncos_Xu_pilim_fill(int lmax, double *where);
-double * taumncos_Xu_pilim_get(int lmax);
+// Legendre functions for negative m; see DLMF 14.9.3
+// 
 
-int pimncos_Xu_zerolim_fill(int lmax, double *where);
-double * pimncos_Xu_zerolim_get(int lmax);
-int pimncos_Xu_pilim_fill(int lmax, double *where);
-double * pimncos_Xu_pilim_get(int lmax);
+
+// tau, pi at singularities and zero; see also DLMF 14.8
+int taumncos_Xu_zerolim_fill(qpms_l_t lmax, double *where)
+{	
+
+	return 0;
+}
+
+double * taumncos_Xu_zerolim_get(qpms_l_t lmax) {
+	double *ar = malloc(qpms_lMax2nelem(lmax) * sizeof(double));
+	taumncos_Xu_zerolim_fill(lmax,ar);
+	return ar;
+}
+
+int taumncos_Xu_pilim_fill(qpms_l_t lmax, double *where);
+
+double * taumncos_Xu_pilim_get(qpms_l_t lmax) {
+	double *ar = malloc(qpms_lMax2nelem(lmax) * sizeof(double));
+	taumncos_Xu_pilim_fill(lmax,ar);
+	return ar;
+}
+
+int taumncos_Xu_pihalflim_fill(qpms_l_t lmax, double *where);
+
+double * taumncos_Xu_pihalflim_get(qpms_l_t lmax) {
+	double *ar = malloc(qpms_lMax2nelem(lmax) * sizeof(double));
+	taumncos_Xu_pihalflim_fill(lmax,ar);
+	return ar;
+}
+
+
+int pimncos_Xu_zerolim_fill(qpms_l_t lmax, double *where);
+double * pimncos_Xu_zerolim_get(qpms_l_t lmax){
+	double *ar = malloc(qpms_lMax2nelem(lmax) * sizeof(double));
+	pimncos_Xu_pihalflim_fill(lmax,ar);
+	return ar;
+}
+int pimncos_Xu_pilim_fill(qpms_l_t lmax, double *where);
+double * pimncos_Xu_pilim_get(qpms_l_t lmax);
 
 
 // [1] (37)
 static inline double pimncos_Xu(double theta) {
-	abort();
+	abort();//NI
 }
 static inline double taumncos_Xu(double theta) {
-	abort();
+	abort();//NI
 }
 
 // [1] (36)
-complex double qpms_vswf_single_mg_Xu(int m, int n, sph_t kdlj,
+complex double qpms_vswf_single_mg_Xu(qpms_m_t m, qpms_l_t n, sph_t kdlj,
 		qpms_bessel_t btyp)
 {
-	
+	abort();//NI	
 }