//compute some cell lists real quick
    for (int pi = 0; pi < _px.size(); pi++)
    {
        int i = (_px[pi]-_offset[0])/_h[0];
        int j = (_py[pi]-_offset[1])/_h[1];
        int k = (_pz[pi]-_offset[2])/_h[2];

        _cells[i*_ny*_nz+j*_nz+k].push_back(pi);
    }

    for (int i = 0; i < _nx; i++) {
        for (int j = 0; j < _ny; j++) {
            for (int k = 0; k < _nz; k++) {

                //bounds
                int low_i = std::max(0,i-1);
                int low_j = std::max(0,j-1);
                int low_k = std::max(0,k-1);

                int high_i = std::min(i+1,_nx-1);
                int high_j = std::min(j+1,_ny-1);
                int high_k = std::min(k+1,_nz-1);

                //+1/-1 support
                for (int ii = low_i; ii <= high_i; ii++)
                    for (int jj = low_j; jj <= high_j; jj++)
                        for (int kk = low_k; kk <= high_k; kk++)
                        {
                            //target cell center
                            float posii = ii * _h[0] + _offset[0] + 0.5f*_h[0];
                            float posjj = jj * _h[1] + _offset[1] + 0.5f*_h[1];
                            float poskk = kk * _h[2] + _offset[2] + 0.5f*_h[2];

                            for (int it = 0; it < _cells[ii*_ny*_nz+jj*_nz+kk].size(); it++)
                            {
                                float px = _px[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];
                                float py = _py[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];
                                float pz = _pz[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];

                                //tent kernel interpol
                                float m_mass = 1.f*tent(px-posii,_h[0])*tent(py-posjj,_h[1])*tent(pz-poskk,_h[2]);

                                _weights[it] += m_mass;
                            }
                        }
            }
        }
    }

    for (int i = 0; i < _nx; i++) {
        for (int j = 0; j < _ny; j++) {
            for (int k = 0; k < _nz; k++) {

                //bounds
                int low_i = std::max(0,i-1);
                int low_j = std::max(0,j-1);
                int low_k = std::max(0,k-1);

                int high_i = std::min(i+1,_nx-1);
                int high_j = std::min(j+1,_ny-1);
                int high_k = std::min(k+1,_nz-1);

                //+1/-1 support
                for (int ii = low_i; ii <= high_i; ii++)
                    for (int jj = low_j; jj <= high_j; jj++)
                        for (int kk = low_k; kk <= high_k; kk++)
                        {
                            //target cell center
                            float posii = ii * _h[0] + _offset[0] + 0.5f*_h[0];
                            float posjj = jj * _h[1] + _offset[1] + 0.5f*_h[1];
                            float poskk = kk * _h[2] + _offset[2] + 0.5f*_h[2];

                            for (int it = 0; it < _cells[ii*_ny*_nz+jj*_nz+kk].size(); it++)
                            {
                                if (_weights[it] == 0.f) continue;

                                float px = _px[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];
                                float py = _py[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];
                                float pz = _pz[_cells[ii*_ny*_nz+jj*_nz+kk].at(it)];

                                //tent kernel interpol
                                float p_mass = _h_part_density[it] * tent(px-posii,_h[0])*tent(py-posjj,_h[1])*tent(pz-poskk,_h[2]) / _weights[it];

                                _h_density->at(ii,jj,kk) += p_mass;
                            }
                        }
            }
        }
    }


    //sanity check
    float max_dens = 0.f;

    for (int i = 0; i < _nx; i++) {
        for (int j = 0; j < _ny; j++) {
            for (int k = 0; k < _nz; k++) {
                max_dens = std::max(max_dens,_h_density->at(i,j,k));
            }
        }
    }