 
         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

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

#include "node_specs.h"


          n_bodies     -= nb;                                   //subtract from n_bodies number of bodies that have been extracted
          nl_pre       += nl;                                   //increase the number of leaves that where processed
          direct [tid] -= nb;                                   //subtract the number of extracted bodies in this pass
          prefix0[tid] = max(prefix0[tid] - nb, 0);             //same here, but do not let the number be negative (GT200 bug!?)
          n_direct     += nb;                                  //increase the number of bodies to be procssed

          while(n_direct >= DIM) 
