QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Minor stuff. 

Former-commit-id: c0de35f6a52ce4cc4f34cf4b417364d92e8dac4d
This commit is contained in:
Marek Nečada 2019-10-05 20:38:04 +03:00
parent aa86fcfb08
commit 2dd9fe2f34
3 changed files with 17 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
qpms/ewaldsf.c
View File

@ -37,6 +37,8 @@ void IgnoreUnderflowsGSLErrorHandler (const char * reason,
// DLMF 8.7.3 (latter expression) for complex second argument // DLMF 8.7.3 (latter expression) for complex second argument
// BTW if a is large negative, it might take a while to evaluate. // BTW if a is large negative, it might take a while to evaluate.
// This can't be used for non-positive integer a due to
// Г(a) in the formula.
int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) { int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) {
if (a <= 0 && a == (int) a) { if (a <= 0 && a == (int) a) {
result->val = NAN + NAN*I; result->val = NAN + NAN*I;
@ -160,7 +162,7 @@ int cx_gamma_inc_CF_e(const double a, const complex double z, qpms_csf_result *r
} }
// incomplete gamma for complex second argument /// Incomplete gamma function for complex second argument.
int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) { int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) {
if (creal(x) >= 0 && if (creal(x) >= 0 &&
(0 == fabs(cimag(x)) || // x is real positive; just use the real fun (0 == fabs(cimag(x)) || // x is real positive; just use the real fun
@ -172,9 +174,14 @@ int complex_gamma_inc_e(double a, complex double x, qpms_csf_result *result) {
return retval; return retval;
} else if (creal(x) >= 0 && cabs(x) > 0.5) } else if (creal(x) >= 0 && cabs(x) > 0.5)
return cx_gamma_inc_CF_e(a, x, result); return cx_gamma_inc_CF_e(a, x, result);
else else if (QPMS_LIKELY(a > 0 || (a % 1.0)))
return cx_gamma_inc_series_e(a, x, result); return cx_gamma_inc_series_e(a, x, result);
else
/* FIXME cx_gamma_inc_series_e() probably fails for non-positive integer a.
* This does not matter for 2D lattices in 3D space,
* but it might cause problems in the other cases.
*/
QPMS_NOT_IMPLEMENTED("Incomplete Gamma function with non-positive integer a.");
} }
// Exponential integral for complex argument; !UNTESTED! and probably not needed, as I expressed everything in terms of inc. gammas anyways. // Exponential integral for complex argument; !UNTESTED! and probably not needed, as I expressed everything in terms of inc. gammas anyways.

6
qpms/qpms_c.pyx
View File

@ -592,6 +592,12 @@ def pitau(double theta, qpms_l_t lMax, double csphase = -1):
qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase) qpms_pitau_fill(&leg[0], &pi[0], &tau[0], theta, lMax, csphase)
return (lega, pia, taua) return (lega, pia, taua)
def linton_gamma(cdouble x):
return clilgamma(x)
def linton_gamma_real(double x):
return lilgamma(x)
def gamma_inc(double a, cdouble x): def gamma_inc(double a, cdouble x):
cdef qpms_csf_result res cdef qpms_csf_result res
with pgsl_ignore_error(15): #15 is underflow with pgsl_ignore_error(15): #15 is underflow

2
qpms/qpms_p.py
View File

@ -1,5 +1,5 @@
import numpy as np import numpy as np
from qpms_c import * from .qpms_c import *
ň = np.newaxis ň = np.newaxis
import scipy import scipy
from scipy.constants import epsilon_0 as ε_0, c, pi as π, e, hbar as ℏ, mu_0 as μ_0 from scipy.constants import epsilon_0 as ε_0, c, pi as π, e, hbar as ℏ, mu_0 as μ_0