/* -*- mode: c++; c-basic-offset: 4 -*- */

#define NO_IMPORT_ARRAY

#include 

// utilities for irregular grids
void _bin_indices_middle(
    unsigned int *irows, int nrows, const float *ys1, unsigned long ny, float dy, float y_min)
{
    int i, j, j_last;
    unsigned int *rowstart = irows;
    const float *ys2 = ys1 + 1;
    const float *yl = ys1 + ny;
    float yo = y_min + dy / 2.0;
    float ym = 0.5f * (*ys1 + *ys2);
    // y/rows
    j = 0;
    j_last = j;
    for (i = 0; i < nrows; i++, yo += dy, rowstart++) {
        while (ys2 != yl && yo > ym) {
            ys1 = ys2;
            ys2 = ys1 + 1;
            ym = 0.5f * (*ys1 + *ys2);
            j++;
        }
        *rowstart = j - j_last;
        j_last = j;
    }
}

void _bin_indices_middle_linear(float *arows,
                                unsigned int *irows,
                                int nrows,
                                const float *y,
                                unsigned long ny,
                                float dy,
                                float y_min)
{
    int i;
    int ii = 0;
    int iilast = (int)ny - 1;
    float sc = 1 / dy;
    int iy0 = (int)floor(sc * (y[ii] - y_min));
    int iy1 = (int)floor(sc * (y[ii + 1] - y_min));
    float invgap = 1.0f / (iy1 - iy0);
    for (i = 0; i < nrows && i  iy1 && ii < iilast) {
            ii++;
            iy0 = iy1;
            iy1 = (int)floor(sc * (y[ii + 1] - y_min));
            invgap = 1.0f / (iy1 - iy0);
        }
        if (i >= iy0 && i  0) {
        int ii = 0;
        int iilast = (int)ny - 1;
        int iy0 = (int)floor(sc * (y[ii] - offs));
        int iy1 = (int)floor(sc * (y[ii + 1] - offs));
        for (i = 0; i < nrows && i < iy0; i++) {
            irows[i] = -1;
        }
        for (; i < nrows; i++) {
            while (i > iy1 && ii < iilast) {
                ii++;
                iy0 = iy1;
                iy1 = (int)floor(sc * (y[ii + 1] - offs));
            }
            if (i >= iy0 && i  iy1 && ii > 1) {
                ii--;
                iy0 = iy1;
                iy1 = (int)floor(sc * (y[ii - 1] - offs));
            }
            if (i >= iy0 && i  0) {
        int ii = 0;
        int iilast = (int)ny - 1;
        int iy0 = (int)floor(sc * (y[ii] - offs));
        int iy1 = (int)floor(sc * (y[ii + 1] - offs));
        float invgap = 1.0 / (iy1 - iy0);
        for (i = 0; i < nrows && i < iy0; i++) {
            irows[i] = -1;
        }
        for (; i < nrows; i++) {
            while (i > iy1 && ii < iilast) {
                ii++;
                iy0 = iy1;
                iy1 = (int)floor(sc * (y[ii + 1] - offs));
                invgap = 1.0 / (iy1 - iy0);
            }
            if (i >= iy0 && i  iy1 && ii > 1) {
                ii--;
                iy0 = iy1;
                iy1 = (int)floor(sc * (y[ii - 1] - offs));
                invgap = 1.0 / (iy1 - iy0);
            }
            if (i >= iy0 && i