def interpolate(xa, ya, queryPoints):
    """
    A cubic spline interpolation on a given set of points (x,y)
    Recalculates everything on every call which is far from efficient but does the job for now
    should eventually be replaced by an external helper class
    """
    n = len(xa)
    u, y2 = [0] * n, [0] * n
    for i in range(1, n - 1):
        wx = xa[i + 1] - xa[i - 1]
        sig = (xa[i] - xa[i - 1]) / wx
        p = sig * y2[i - 1] + 2.0
        y2[i] = (sig - 1.0) / p
        ddydx = (ya[i + 1] - ya[i]) / (xa[i + 1] - xa[i]) - (ya[i] - ya[i - 1]
            ) / (xa[i] - xa[i - 1])
        u[i] = (6.0 * ddydx / wx - sig * u[i - 1]) / p
    y2[n - 1] = 0
    for i in range(n - 2, -1, -1):
        y2[i] = y2[i] * y2[i + 1] + u[i]
    results = [0] * n
    for i in range(len(queryPoints)):
        klo = 0
        khi = n - 1
        while khi - klo > 1:
            k = khi + klo >> 1
            if xa[k] > queryPoints[i]:
                khi = k
            else:
                klo = k
        h = xa[khi] - xa[klo]
        a = (xa[khi] - queryPoints[i]) / h
        b = (queryPoints[i] - xa[klo]) / h
        results[i] = a * ya[klo] + b * ya[khi] + ((a * a * a - a) * y2[klo] +
            (b * b * b - b) * y2[khi]) * (h * h) / 6.0
    return results
