Commit bddd29f1 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

implement correction factors

parent 4893a233
......@@ -148,6 +148,63 @@ inline double fmodulo (double v1, double v2)
// return (v1>=0) ? ((v1<v2) ? v1 : fmod(v1,v2)) : (fmod(v1,v2)+v2);
}
static inline double one_minus_x2 (double x)
{ return (fabs(x)>0.1) ? (1.+x)*(1.-x) : 1.-x*x; }
void legendre_prep(int n, vector<double> &x, vector<double> &w)
{
constexpr double pi = 3.141592653589793238462643383279502884197;
constexpr double eps = 3e-14;
int m = (n+1)>>1;
x.resize(m);
w.resize(m);
double t0 = 1 - (1-1./n) / (8.*n*n);
double t1 = 1./(4.*n+2.);
#pragma omp parallel
{
int i;
#pragma omp for schedule(dynamic,100)
for (i=1; i<=m; ++i)
{
double x0 = cos(pi * ((i<<2)-1) * t1) * t0;
int dobreak=0;
int j=0;
double dpdx;
while(1)
{
double P_1 = 1.0;
double P0 = x0;
double dx, x1;
for (int k=2; k<=n; k++)
{
double P_2 = P_1;
P_1 = P0;
// P0 = ((2*k-1)*x0*P_1-(k-1)*P_2)/k;
P0 = x0*P_1 + (k-1.)/k * (x0*P_1-P_2);
}
dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
/* Newton step */
x1 = x0 - P0/dpdx;
dx = x0-x1;
x0 = x1;
if (dobreak) break;
if (abs(dx)<=eps) dobreak=1;
if (++j>=100) throw runtime_error("convergence problem");
}
x[m-i] = x0;
w[m-i] = 2. / (one_minus_x2(x0) * dpdx * dpdx);
}
} // end of parallel region
}
using a_i_c = py::array_t<int, py::array::c_style | py::array::forcecast>;
using a_d_c = py::array_t<double, py::array::c_style | py::array::forcecast>;
using a_c_c = py::array_t<complex<double>,
......@@ -450,6 +507,31 @@ a_c_c from_grid_pre (const a_d_c &grid_)
return res;
}
a_d_c correction_factors (size_t n, size_t nval, double epsilon)
{
auto w = int(-log10(epsilon)+1.9999);
auto beta = 2.3*w;
auto p = int(1.5*w+2);
double h = 2*pi/n;
double alpha = pi*w/n;
vector<double> x, wgt;
legendre_prep(2*p,x,wgt);
auto psi = x;
for (auto &v:psi)
v = exp(beta*(sqrt(1-v*v)-1.));
int odim[] = {nval};
a_d_c res(odim);
auto val = res.mutable_data();
for (size_t k=0; k<nval; ++k)
{
double tmp=0;
for (int i=0; i<p; ++i)
tmp += wgt[i]*psi[i]*cos(alpha*k*x[i]);
val[k] = 1./(w*tmp);
}
return res;
}
} // unnamed namespace
PYBIND11_MODULE(nifty_gridder, m)
......@@ -461,4 +543,5 @@ PYBIND11_MODULE(nifty_gridder, m)
m.def("to_grid_post",&to_grid_post);
m.def("from_grid",&from_grid);
m.def("from_grid_pre",&from_grid_pre);
m.def("correction_factors",&correction_factors);
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment