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 + +#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 #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 + +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