C indexing functions to separate file, new normalizations in progress
Former-commit-id: 30b4ab302fe69fedc8c12bdb7b3b35f532cfc35d
This commit is contained in:
Marek Nečada
parent
55d618143b
commit
bb289a46b9
|
|
@ -0,0 +1,27 @@
|
|||
#ifndef QPMS_INDEXING_H
|
||||
#define QPMS_INDEXING_H
|
||||
|
||||
#include "qpms_types.h"
|
||||
|
||||
static inline qpms_y_t qpms_mn2y(qpms_m_t m, qpms_l_t n) {
|
||||
return (qpms_y_t) n * (n + 1) + m - 1;
|
||||
}
|
||||
|
||||
static inline qpms_lm_t qpms_y2n(qpms_y_t y) {
|
||||
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
|
||||
return sqrt(y+1);
|
||||
}
|
||||
|
||||
static inline qpms_m_t qpms_yn2m(qpms_y_t y, qpms_l_t n) {
|
||||
return y-qpms_mn2y(0,n);
|
||||
}
|
||||
|
||||
static inline void qpms_y2mn_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
|
||||
*m=qpms_yn2m(y,*n=qpms_y2n(y));
|
||||
}
|
||||
|
||||
static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
|
||||
return lmax * ((qpms_y_t)lmax + 2);
|
||||
}
|
||||
|
||||
#endif //QPMS_INDEXING_H
|
||||
|
|
@ -7,6 +7,8 @@
|
|||
|
||||
// integer index types
|
||||
typedef int qpms_lm_t;
|
||||
typedef unsigned int qpms_l_t;
|
||||
typedef qpms_lm_t qpms_m_t;
|
||||
typedef size_t qpms_y_t;
|
||||
|
||||
// Normalisations
|
||||
|
|
|
|||
|
|
@ -1,11 +1,17 @@
|
|||
#include <math.h>
|
||||
#include "gaunt.h"
|
||||
#include "translations.h"
|
||||
#include "indexing.h" // TODO replace size_t and int with own index types here
|
||||
#include <stdbool.h>
|
||||
#include <gsl/gsl_sf_legendre.h>
|
||||
#include <gsl/gsl_sf_bessel.h>
|
||||
#include "assert_cython_workaround.h"
|
||||
#include <stdlib.h> //abort()
|
||||
|
||||
/*
|
||||
* References:
|
||||
* [1] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996)
|
||||
*/
|
||||
|
||||
static const double sqrtpi = 1.7724538509055160272981674833411451827975494561223871;
|
||||
//static const double ln2 = 0.693147180559945309417232121458176568075500134360255254120;
|
||||
|
|
@ -61,9 +67,10 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double
|
|||
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
|
||||
// [1], eq. (82)
|
||||
complex double qpms_trans_single_A_Xu(int m, int n, int mu, int nu, sph_t kdlj,
|
||||
bool r_ge_d, qpms_bessel_t J) {
|
||||
abort(); // FIXME, THIS IS STILL TAYLOR
|
||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
||||
|
||||
double costheta = cos(kdlj.theta);
|
||||
|
|
@ -105,6 +112,155 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
|
|||
return (presum / prenormratio) * sum;
|
||||
}
|
||||
|
||||
|
||||
|
||||
complex double qpms_trans_single_A_Kristensson(int m, int n, int mu, int nu, sph_t kdlj,
|
||||
bool r_ge_d, qpms_bessel_t J) {
|
||||
abort();// FIXME, THIS IS STILL TAYLOR
|
||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
||||
|
||||
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[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;
|
||||
int Pp_order = mu-m;
|
||||
//if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze
|
||||
assert(p >= abs(Pp_order));
|
||||
double a1q_n = a1q[q] / a1q0;
|
||||
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;
|
||||
}
|
||||
|
||||
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));
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
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) {
|
||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
||||
|
||||
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[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;
|
||||
int Pp_order = mu-m;
|
||||
//if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze
|
||||
assert(p >= abs(Pp_order));
|
||||
double a1q_n = a1q[q] / a1q0;
|
||||
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;
|
||||
}
|
||||
|
||||
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));
|
||||
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;
|
||||
}
|
||||
|
||||
// [1], eq. (83)
|
||||
complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj,
|
||||
bool r_ge_d, qpms_bessel_t J) {
|
||||
abort(); // FIXME, this is still Taylor
|
||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
||||
double costheta = cos(kdlj.theta);
|
||||
|
||||
int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
|
||||
int Qmax = gaunt_q_max(-m,n+1,mu,nu);
|
||||
double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0;
|
||||
int err;
|
||||
if (mu == nu) {
|
||||
for (int q = 0; q <= q2max; ++q)
|
||||
a2q[q] = 0;
|
||||
a2q0 = 1;
|
||||
}
|
||||
else {
|
||||
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
|
||||
a2q0 = a2q[0];
|
||||
}
|
||||
gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort();
|
||||
a3q0 = a3q[0];
|
||||
|
||||
double leg[gsl_sf_legendre_array_n(n+nu+1)];
|
||||
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
|
||||
complex double bes[n+nu+2];
|
||||
if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort();
|
||||
|
||||
complex double sum = 0;
|
||||
for (int q = 0; q <= Qmax; ++q) {
|
||||
int p = n+nu-2*q;
|
||||
double a2q_n = a2q[q]/a2q0;
|
||||
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
|
||||
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
|
||||
-(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
|
||||
*min1pow(q) * zp_ * Pp_);
|
||||
sum += summandq;
|
||||
}
|
||||
|
||||
double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
|
||||
+lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
|
||||
+lgamma(n+nu+2) - lgamma(2*(n+nu)+3));
|
||||
complex double presum = exp(exponent);
|
||||
presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / (
|
||||
(4*n)*(n+1)*(n+m+1));
|
||||
|
||||
// Taylor normalisation v2, proven to be equivalent
|
||||
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)));
|
||||
|
||||
|
||||
|
||||
complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
|
||||
bool r_ge_d, qpms_bessel_t J) {
|
||||
if(r_ge_d) J = QPMS_BESSEL_REGULAR;
|
||||
|
|
@ -186,19 +342,6 @@ static inline size_t qpms_mn2y(int m, int n) {
|
|||
return (size_t) n * (n + 1) + m - 1;
|
||||
}
|
||||
|
||||
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
|
||||
return sqrt(y+1);
|
||||
}
|
||||
|
||||
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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue