makeacg.py

import numpy


def acg(N):
    """ approximated confined gaussian window of support N;

        https://en.wikipedia.org/wiki/Window_function#Confined_Gaussian_window
    """
    # sigma = 1.0
    s = 1.0
    # edge
    A = (N - 1) / 2.0
    x = numpy.linspace(0, N * 0.5, 8192, endpoint=True)
    y = x + A

    def G(y):
        return numpy.exp(-0.25 * ((y - A)/ s) ** 2)

    phi = G(y) - G(-0.5) * (G(y + N) + G(y - N)) / (G(-0.5 + N) + G(-0.5 - N))

    sum = 2 * numpy.trapz(phi, x)
    phi /= sum
    print(N, phi[-1])
    return phi, x

def genacg(n):
    phi, x = acg(n)
    name = 'acg%d' % n
    support = n

    vnumbers = ["%.8f, %.8f, %.8f, %.8f" % tuple(a) for a in phi.reshape(-1, 4)]
    step = numpy.diff(x).mean()

    template = """
    static double _%(funcname)s_vtable[] = %(vtable)s;
    static double _%(funcname)s_nativesupport = %(support)g;
    static double _%(funcname)s_kernel(double x)
    {
        x = fabs(x);
        double f = x / %(step)e;
        int i = f;
        if (i < 0) return 0;
        if (i >= %(tablesize)d - 1) return 0;
        f -= i;
        return _%(funcname)s_vtable[i] * (1 - f)
             + _%(funcname)s_vtable[i+1] * f;
    }
    static double _%(funcname)s_diff(double x)
    {
        double factor;
        if(x >= 0) {
            factor = 1;
        } else {
            factor = -1;
            x = -x;
        }

        int i = x / %(step)e;
        if (i < 0) return 0;
        if (i >= %(tablesize)d - 1) return 0;
        double f = _%(funcname)s_vtable[i+1] - _%(funcname)s_vtable[i];
        return factor * f / %(step)e;
    }
    """
    return template % {
            'vtable' : "{\n" + ",\n".join(vnumbers) + "}",
            'hsupport' : support * 0.5,
            'support' : support,
            'funcname' : name,
            'step' : step,
            'tablesize' : len(phi),
    }

with open('pmesh/_window_acg.h', 'wt') as f:
    f.write("""
/*
 * do not modify this file
 * generated by makeacg.py
 *
 */
""")
    f.write(genacg(2))
    f.write(genacg(3))
    f.write(genacg(4))
    f.write(genacg(5))
    f.write(genacg(6))