From b8eff3eb77426ae1b34f15a4e737816b960e999e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek.necada@aalto.fi>
Date: Mon, 10 Apr 2017 13:48:06 +0300
Subject: [PATCH] Dudom

Former-commit-id: cf2eb302d096612df04c204ca0dcee71a0776672
---
 qpms/gaunt.c        | 14 ++++++++++++--
 qpms/gaunt.h        | 11 ++++++++++-
 qpms/translations.c | 15 +++++++++++++++
 qpms/vectors.h      | 46 +++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 83 insertions(+), 3 deletions(-)
 create mode 100644 qpms/vectors.h

diff --git a/qpms/gaunt.c b/qpms/gaunt.c
index 146fa57..4652b49 100644
--- a/qpms/gaunt.c
+++ b/qpms/gaunt.c
@@ -1167,13 +1167,23 @@ void gaunt_xu(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error)
 	} // qmax_case
 } // gaunt_xu
 
-#ifdef GAUNTTEST
+	
 #define MIN(x,y) (((x) > (y)) ? (y) : (x))
-		
 static inline int q_max(int m, int n, int mu, int nu) {
 	return MIN(n, MIN(nu,(n+nu-abs(m+mu))/2));
 }
 
+int gaunt(int m, int n, int mu, int nu, double *v) {
+	int err = 0;
+	int qmax = q_max(m,n,mu,nu);
+	if (!v) v = calloc(qmax+1, sizeof(double));
+	if (!v) return -1;
+	gaunt_xu(m, n, mu, nu, qmax, v, &err);
+	return err;
+}
+
+#ifdef GAUNTTEST
+	
 void __vec_trans_MOD_gaunt_xu(const double *m, const double *n, const double *mu, const double *nu, const int *qmax, double *v_aq, int *err);
 int main()
 {
diff --git a/qpms/gaunt.h b/qpms/gaunt.h
index bc6b99c..8536075 100644
--- a/qpms/gaunt.h
+++ b/qpms/gaunt.h
@@ -1,4 +1,13 @@
 #ifndef GAUNT_H
 #define GAUNT_H
-void gaunt_xu(int m, int n, int mu, int nu, int qmax, double v_aq, int *err);
+#include <stdlib.h>
+
+#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);
 #endif //GAUNT_H
diff --git a/qpms/translations.c b/qpms/translations.c
index 18018e0..bb170de 100644
--- a/qpms/translations.c
+++ b/qpms/translations.c
@@ -1,2 +1,17 @@
 #include <math.h>
 #include "gaunt.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);
+}
+
+// Derivative of associated Legendre polynomial at zero argument (DLMF 14.5.2)
+double legendreD0(int m, int n) {
+	return -2 * legendre0(m, n);
+}
+
+
diff --git a/qpms/vectors.h b/qpms/vectors.h
new file mode 100644
index 0000000..eafbfaa
--- /dev/null
+++ b/qpms/vectors.h
@@ -0,0 +1,46 @@
+#ifndef VECTORS_H
+#define VECTORS_H
+#include <math.h>
+
+typedef struct {
+	double x, y, z;
+} cart3_t;
+
+typedef struct {
+	double x, y;
+} cart2_t;
+
+typedef struct {
+	double r, theta, phi;
+} sph_t;
+
+typedef struct {
+	double r, phi;
+} pol_t;
+
+
+//static inline double vectors_h_sq(double x) {return x*x;}
+
+static inline cart3norm(cart3_t v) {
+	return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+}
+
+static inline sph_t cart2sph(cart3_t cart) {
+	sph_t sph;
+	sph.r = cart3norm(cart);
+	sph.theta = sph.r ? acos(cart.z / sph.r) : M_PI_2;
+	sph.phi = atan2(cart.y, cart.x);
+	return sph;
+}
+
+static inline cart_t sph2cart(sph_t sph) {
+	cart_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);
+	cart.z = sph.r * cos(sph.theta);
+	return cart;
+}
+
+
+#endif //VECTORS_H