 
   //  begin,    end,   size
   // -----------------------
-#ifdef _ORIG_SHMEM_
-
-  int *approxS = (int*)&shmem  [     0];            //  0*DIM,  2*DIM,  2*DIM
-  int *directS = (int*)&approxS[ 2*DIM];            //  2*DIM,  3*DIM,  1*DIM
-  int *nodesS = (int*)&directS [   DIM];            //  3*DIM, 12*DIM,  9*DIM
-  int *prefix = (int*)&nodesS  [9 *DIM];            // 12*DIM, 14*DIM,  2*DIM
-  int *sh_body = &approxS[DIM];
-
-  int *prefix0 = &prefix[  0];
-  int *prefix1 = &prefix[DIM];
-
-  const int NJMAX = DIM*2;
-  int    *body_list = (int*   )&nodesS   [  DIM]; //  4*DIM,  6*DIM,  2*DIM
-  float  *sh_mass   = (float* )&body_list[NJMAX]; //  6*DIM,  7*DIM,  1*DIM
-  float3 *sh_pos    = (float3*)&sh_mass  [  DIM]; //  7*DIM, 10*DIM   3*DIM
-
-  int *approxM = approxS;
-  int *directM = directS;
-  int * nodesM =  nodesS;
-
-#else   /* !_ORIG_SHMEM_ */
-
   const int stack_sz = (LMEM_STACK_SIZE << SHIFT) << DIM2;  /* stack allocated per thread-block */
   int *approxL = lmem + stack_sz; 
 
   int *directS = shmem;                              //  0*DIM,  1*DIM,  1*DIM
   int *nodesS  = directS + DIM;                      //  1*DIM, 10*DIM,  9*DIM
-  int *prefixS = nodesS  + DIM*9;                    // 10*DIM, 12*DIM,  2*DIM
-
-  int *prefix  = prefixS;
-  int *prefix0 = prefixS;
-  int *prefix1 = prefixS + DIM;
+  int *prefix  = nodesS  + DIM*9;                    // 10*DIM, 11*DIM,  1*DIM
 
   const int NJMAX = DIM*3;
   int    *body_list = (int*   )&nodesS   [  DIM]; //  2*DIM,   5*DIM,  2*DIM
   float  *sh_mass   = (float* )&body_list[NJMAX]; //  5*DIM,   6*DIM,  1*DIM
   float3 *sh_pos    = (float3*)&sh_mass  [  DIM]; //  6*DIM,   9*DIM   3*DIM
-  int    *sh_body   = nodesS + DIM*8;             //  9*DIM,  10*DIM,  1*DIM
 
   int *approxM = approxL;
   int *directM = directS;
   int * nodesM =  nodesS;
-
-#endif /* _ORIG_SHMEM_ */
 
 
   float  *node_mon0 = sh_mass;

          //the amount of allocated shared memory

          // step 0                                                      //nullify part of the body_list that will be filled with bodies
          for (int i = n_direct; i < n_direct + nb; i += DIM){           //from the leaves that are being processed
            body_list[i + tid] = 0;
          }
          __syncthreads();

          //step 1:
          if (flag && (directM[tid] <= nb) && (prefix0[tid] > 0))        //make sure that the thread indeed carries a leaf
            body_list[n_direct + prefix0[tid] - 1] = 1;                 //whose bodies will be extracted
          __syncthreads();

          //step 2:
          int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);   // inclusive scan to compute number of leaves to process
          // to make sure that there is enough shared memory for bodies
          nb = directM[prefix1[nl_pre + nl - 1]];                       // number of bodies stored in these leaves

          // step 3:
          for (int i = n_direct; i < n_direct + nb; i += DIM) {          //segmented fill of the body_list
            int j = prefix1[nl_pre + body_list[i + tid] - 1];            // compute the first body in shared j-body array
            body_list[i + tid] = (i + tid - n_direct) -                 //add to the index of the first j-body in a child
              (prefix0[j] - 1) + sh_body[j];         //the index of the first child in body_list array
          }
          __syncthreads();


          /**************************************************
           *  example of what is accomplished in steps 0-4   *
           *       ---------------------------               *
           * step 0: body_list = 000000000000000000000       *
           * step 1: body_list = 100010001000000100100       *
           * step 2: body_list = 111122223333333444555       *
           * step 3: body_list = 012301230123456012012       *
           *         assuming that sh_body[j] = 0            *
           ***************************************************/

          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
          directM[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) 

PROF_MODULE(dev_approximate_gravity);

#include "node_specs.h"

#ifdef WIN32
#define M_PI        3.14159265358979323846264338328
#endif

#define WARP_SIZE2 5
#define WARP_SIZE  32

#if 1
#define laneId (threadIdx.x & (WARP_SIZE - 1))
#define warpId (threadIdx.x >> WARP_SIZE2)
#endif

__forceinline__ __device__ float Wkernel(const float q)
{
