-// #include "support_kernels.cu"
+  // #include "support_kernels.cu"
 #include <stdio.h>
-#include "../profiling/bonsai_timing.h"
-PROF_MODULE(dev_approximate_gravity);
 
 #include "node_specs.h"
 

        int *prefix0 = &prefix[  0];
        int *prefix1 = &prefix[DIM];

#ifdef OLDPREFIX
        int n_total = calc_prefix<DIM2>(prefix, tid,  nchild);
        prefix[tid] += n_offset - nchild;
        __syncthreads();
#else
        int2 offset2 = inclusive_scan_blockS<DIM2>(prefix_shmem, nchild, tid);        // inclusive scan to compute memory offset of each child
        int  offset = offset2.x;
        int n_total = offset2.y;
        offset += n_offset - nchild;
#endif

        for (int i = n_offset; i < n_offset + n_total; i += DIM)         //nullify part of the array that will be filled with children
          nodes[tid + i] = 0;                                          //but do not touch those parts which has already been filled
        __syncthreads();                                                 //Thread barrier to make sure all warps finished writing data

        bool flag = (split && !leaf) && use_node;                        //Flag = use_node + split + not_a_leaf;Use only non_leaf nodes that are to be split
#if 0
        if (flag) nodes[prefix[tid]] = child;                            //Thread with the node that is about to be split
        __syncthreads();                                                 //writes the first child in the array of nodes

        /*** in the following 8 lines, we calculate indexes of all the children that have to be walked from the index of the first child***/
        if (flag && nodes[prefix[tid] + 1] == 0) nodes[prefix[tid] + 1] = child + 1; __syncthreads();
        if (flag && nodes[prefix[tid] + 2] == 0) nodes[prefix[tid] + 2] = child + 2; __syncthreads();
        if (flag && nodes[prefix[tid] + 3] == 0) nodes[prefix[tid] + 3] = child + 3; __syncthreads();
        if (flag && nodes[prefix[tid] + 4] == 0) nodes[prefix[tid] + 4] = child + 4; __syncthreads();
        if (flag && nodes[prefix[tid] + 5] == 0) nodes[prefix[tid] + 5] = child + 5; __syncthreads();
        if (flag && nodes[prefix[tid] + 6] == 0) nodes[prefix[tid] + 6] = child + 6; __syncthreads();
        if (flag && nodes[prefix[tid] + 7] == 0) nodes[prefix[tid] + 7] = child + 7; __syncthreads();
#else
        /* Thread with the node that is about to be split */
        if (flag                          ) nodes[offset    ] = child;    
        /* in the following 8 lines, we calculate indexes of all the children that have to be walked from the index of the first child */
        if (flag && nodes[offset + 1] == 0) nodes[offset + 1] = child + 1;
        if (flag && nodes[offset + 2] == 0) nodes[offset + 2] = child + 2;
        if (flag && nodes[offset + 3] == 0) nodes[offset + 3] = child + 3;
        if (flag && nodes[offset + 4] == 0) nodes[offset + 4] = child + 4;

texture<float4, 1, cudaReadModeElementType> texMultipole;
texture<float4, 1, cudaReadModeElementType> texBody;

#if 1
template<class T>
 struct ADDOP {
  __device__ static inline T identity()           {return (T)(0);}
