 }
 
 
-
-/*static __device__ float4 get_D04(float ds2, int selfGrav = 1) {
-#if 1
-  float ids  = rsqrtf(ds2);  //Does not work with zero-softening
-  //   if(isnan(ids)) ids = 0;               //This does work with zero-softening, few percent performance drop
-  //float ids  = (1.0f / sqrtf(ds2)) * selfGrav; Slower in Pre CUDA4.1
-  ids *= selfGrav;
-#else
-  const float ids = (ds2 > 0.0f) ? rsqrtf(ds2) : 0.0f;
-#endif
-  const float ids2 = ids*ids;
-  float ids3 = ids *ids2;
-  float ids5 = ids3*ids2;
-  float ids7 = ids5*ids2;
-  return make_float4(ids, -ids3, +3.0f*ids5, -15.0f*ids7);
-}  // 9 flops*/
 #ifdef _QUADRUPOLE_
-
 
 static __device__ __forceinline__ float4 add_acc(
     float4 acc, 

    }

    bid   = shmemv[0];
    if (bid >= n_active_groups) return;

    int grpOffset = 0;

    /*********** set necessary thread constants **********/
#ifdef DO_BLOCK_TIMESTEP
    real4 curGroupSize    = groupSizeInfo[active_groups[bid + grpOffset]];
#else
    real4 curGroupSize    = groupSizeInfo[bid + grpOffset];
#endif
    const int   groupData       = __float_as_int(curGroupSize.w);
    const uint body_addr        =   groupData & CRITMASK;
    const uint nb_i             = ((groupData & INVCMASK) >> CRITBIT) + 1;

#ifdef DO_BLOCK_TIMESTEP
    real4 group_pos       = groupCenterInfo[active_groups[bid + grpOffset]];
#else
    real4 group_pos       = groupCenterInfo[bid + grpOffset];
#endif

    uint body_i[2];
    int ni = nb_i <= WARP_SIZE ? 1 : 2;
    body_i[0] = body_addr + laneId%nb_i;
    body_i[1] = body_addr + WARP_SIZE + laneId%(nb_i - WARP_SIZE);

    float4 pos_i[2];
    float4 acc_i[2];

    pos_i[0] = group_body_pos[body_i[0]];
    if(ni > 1) //Only read if we actually have ni == 2
      pos_i[1] = group_body_pos[body_i[1]];

    acc_i[0] = acc_i[1] = make_float4(0.0f, 0.0f, 0.0f, 0.0f);



    int apprCount = 0;
    int direCount = 0;

    const float group_eps  = 0;
    int ngb_i;
#if 1
    if (ni == 1)
      approximate_gravity<0, blockDim2, 1>(
          pos_i, group_pos,
          eps2, node_begend,
          multipole_data, body_pos,
          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
          group_eps, 
          acc_i);
    else
      approximate_gravity<0, blockDim2, 2>(
          pos_i, group_pos,
          eps2, node_begend,
          multipole_data, body_pos,
          shmem, lmem, ngb_i, apprCount, direCount, boxSizeInfo, curGroupSize, boxCenterInfo,
          group_eps, 
          acc_i);
#endif

    assert(apprCount >= 0);
#if 0 /* this increase lmem spill count */
    if(apprCount < 0)
    {

      //Try to get access to the big stack, only one block per time is allowed
      if(laneId == 0)
      {
        int res = atomicExch(&active_inout[n_bodies+1], 1); //If the old value (res) is 0 we can go otherwise sleep


#define TEXTURES

template<int SHIFT, int BLOCKDIM2, int NI, bool INTCOUNT>
static __device__ 
uint2 approximate_gravity(
    float4 acc_i[NI],
    const float4 _pos_i[NI],
    const float4 groupPos,
