 
         n_approx += n_total;
 
+
         while (n_approx >= DIM) 
         {
           n_approx -= DIM;
           const int address      = (approx[n_approx + tid] << 1) + approx[n_approx + tid];
 #ifndef TEXTURES
           const float4 monopole  = multipole_data[address    ];
-#if 0
-          float4 octopole0 = multipole_data[address + 1];
-          float4 octopole1 = multipole_data[address + 2];
-#endif
 #else
           const float4 monopole  = tex1Dfetch(texMultipole, address);
-#if 0
-          float4 octopole0 = tex1Dfetch(texMultipole, address + 1);
-          float4 octopole1 = tex1Dfetch(texMultipole, address + 2);
-#endif
 #endif
 
           node_mon0[tid] = monopole.w;
           node_mon1[tid] = make_float3(monopole.x,  monopole.y,  monopole.z);
           __syncthreads();
-
-#if 0
-          const float f_dm   = 0.0f;
-          const float f_star = 1.0f
-            const float darkMatterMass = f_dm   * octopole1.w;
-          /* eg: we need to be careful with the line below to avoid truncation error due to 
-             subtraction of two large numbers, monopole.w and darkMatterMass both could be
-             very large.
-             Instead, we can use octopole1.w to be stellar mass, and DM mass to be 
-             monopole.w, then we add the two together to get total mass, but this will
-             require more changes to the kernel */
-          const float    stellarMass = f_star * (monopole.w - darkMatterMass);
-          const float hinv = 1.0f/hi;   /* eg: this can be precomputing to avoid division */
-          density += interact(
-              make_float3(pos_i.x, pos_i.y, pos_i.z), h, hinv,
-              make_float3(monopole.x, monople.y, monopole.z), darkMatterMass + stellarMass);
-          /* eg: the interact function still calls sqrtf(f), which invloves 1 div and 1 rsqrtf,
-             so ideally we would like to take advantage of rsqrtf in add_acc, and then we only
-             do 1 div */
-#endif
-
 
 #if 1
 #pragma unroll 16

        body_list[tid] = direct[tid];                                            //copy list of bodies from previous pass to body_list
        sh_body  [tid] = jbody;                                                  //store the leafs first body id into shared memory

        // step 1
#if 0
        int v0 = inclusive_scan_block<DIM2>(prefix_shmem, (int)flag, tid);       // inclusive scan on flags to construct array
        prefix0[tid] = v0;
        __syncthreads();

        if (flag) prefix1[prefix0[tid] - 1] = tid;                             //with tidś whose leaves have to be opened
        __syncthreads();                                                      //thread barrier, make sure all warps completed the job
#else
        offset = inclusive_scan_block<DIM2>(prefix_shmem, (int)flag, tid);       // inclusive scan on flags to construct array
        if (flag) prefix1[offset - 1] = tid;                             //with tidś whose leaves have to be opened
#endif

        // step 2
#if 0
        int v0 = inclusive_scan_block<DIM2>(prefix_shmem, nbody, tid);        // inclusive scan to compute memory offset for each body
        prefix0[tid] = v0;
        __syncthreads();
        int n_bodies = prefix0[blockDim.x - 1];                            //Total number of bides extract from the leaves
        __syncthreads();                                                   // thread barrier to make sure that warps completed their jobs

        direct [tid]  = prefix0[tid];                                       //Store a copy of inclusive scan in direct
        prefix0[tid] -= nbody;                                              //convert inclusive int oexclusive scan
        prefix0[tid] += 1;                                                  //add unity, since later prefix0[tid] == 0 used to check barrier
#else
        offset2 = inclusive_scan_blockS<DIM2>(prefix_shmem, nbody, tid);        // inclusive scan to compute memory offset for each body
        int offset1  = offset2.x;
        int n_bodies = offset2.y;

        direct[tid] = offset1;
        offset1    -= nbody;
        offset1    += 1;
        prefix0[tid] = offset1;
        __syncthreads();
#endif

        int nl_pre = 0;                                                     //Number of leaves that have already been processed



}

#endif

#if 1

#define WARP_SIZE2 5
#define WARP_SIZE  (1 << WARP_SIZE2)

#if 0

__device__ __forceinline__ uint shfl_scan_add_step(uint partial, uint up_offset)
{
