def METHOD_NAME(Nin, xin, yin, Np, xp, yp, xo, dx, _in):
    # Walk the theory spline and the resolution spline together, computing
    # the integral of the pairs of line segments.  Since the spline knots
    # do not necessarily correspond, we need to treat each individual segment
    # in both curves piece-wise, integrating from knot to knot in the union
    # of the two knot sets.
    # Need an undefined value in case the algorithm is broken and we don't
    # initialize the initial lines.  This would ideal be NaN, but that isn't
    # available in MSVC, so set it to something very large.  In order to avoid
    # normalizing that very large value to 1., only set the theory line to
    # undefined, and set the resolution line to 0.
    # const double undefined = 0./0.;
    undefined = 1e308
    m1 = undefined
    b1 = undefined
    m2 = 0.
    b2 = 0.
    sum = 0.
    norm = 0.
    # size_t p;
    # double delta,delta2,delta3;
    # double x, next_x, next_xin, next_xp;
    # Set the target start of the integral to xin[_in].  This may or may
    # not be beyond the end of the resolution function, depending on whether
    # we are at the left edge of the data, or somewhere inside.
    next_xin = xin[_in]
    # Set p to the point just before next_xin
    for p in range(1, Np):
        if (xo + dx*xp[p] > next_xin):
            break
    p -= 1
    next_xp = xo + dx*xp[p]
    # Choose the larger of next_xp and next_xin as the starting point of the
    # integral.  This will force us to step both xin and xp in the first
    # iteration of the loop, computing new slope/intercepts for both lines.
    # If the theory line extends beyond the resolution function, then we will
    # be called with xin[_in] before the first resolution point, and next_xp
    # will be bigger than next_xin.  If the resolution extends beyond the
    # end of the theory, then the next_xp will be set before next_xin, and
    # next_xin will be the start of the integral.  The integral ends when
    # either the resolution or the theory runs out.
    #
    # We are tracking the area under the resolution as well as the area
    # under the product of theory and resolution so that we can properly
    # normalize the data when less than the full resolution is included.
    # This means that we do not need to provide a resolution function
    # normalized to a total area of 1 as our input.
    x = next_xp if next_xp > next_xin else next_xin
    # printf("  point xo:%g dx:%g x:%g, p:%ld, in:%ld, xp:%g, xin:%g\n",
    #       xo,dx,x,p,in,next_xp,next_xin);
    while True:
        # Step xin if we are at the next theory point
        if (next_xin <= x):
            _in += 1
            if (_in >= Nin):
                break  # At the right edge of the data
            next_xin = xin[_in]
            m1 = (yin[_in] - yin[_in-1]) / (xin[_in] - xin[_in-1])
            b1 = yin[_in] - m1*xin[_in]
        # Step xp if we are at the next resolution point
        if (next_xp <= x):
            p += 1
            if (p >= Np):
                break  # At the right edge of the resolution
            next_xp = xo + dx*xp[p]
            m2 = (yp[p] - yp[p-1]) / (xp[p] - xp[p-1]) / dx
            b2 = yp[p] - m2*next_xp
        # Find the next node
        next_x = next_xin if next_xin < next_xp else next_xp
        # Compute the convolution and norm between the current and next node
        delta = next_x - x
        delta2 = next_x*next_x - x*x
        delta3 = next_x*next_x*next_x - x*x*x
        norm += 0.5*m2*delta2 + b2*delta
        sum += m1*m2/3.0*delta3 + 0.5*(m1*b2+m2*b1)*delta2 + b1*b2*delta
        # printf("  delta:%g delta2:%g delta3:%g norm:%g sum:%g m1:%g b1:%g m2:%g b2:%g x:%g nx:%g ni:%g np:%g\n",
        #  delta, delta2, delta3, norm, sum, m1, b1, m2, b2, x, next_x, next_xin, next_xp);
        # Move to the next node
        x = next_x
    return sum/norm