 
   /*********** shared memory distribution **********/
 
-                                                  //  begin,    end,   size
-                                                  // -----------------------
+  //  begin,    end,   size
+  // -----------------------
   int *approx = (int*)&shmem [     0];            //  0*DIM,  2*DIM,  2*DIM
   int *direct = (int*)&approx[ 2*DIM];            //  2*DIM,  3*DIM,  1*DIM
   int *nodes  = (int*)&direct[   DIM];            //  3*DIM, 13*DIM, 10*DIM


}

#ifdef INDSOFT
__device__ float4 get_D04(float ds2, float epsP, float epsQ, int selfGrav = 1) {

//   float eps     = fmaxf(epsP, epsQ);
  float epseff  = epsP + epsQ;
  float ids, ids2, ids3;

  if(ds2 >= (epseff*epseff))
  {
     ids  = rsqrtf(ds2);
     //if(isnan(ids)) ids = 0; not needed if we use non-zero softening
     ids  = ids*selfGrav;       //Prevent selfGravity, instead of using if-statement
     ids3 = ids*ids*ids;
  }
  else
  {
    //these two lines are faster than a real sqrt
    float dist = ds2*selfGrav*rsqrtf(ds2); //Gives NaN is ds is 0
    if(isnan(dist)) dist = 0.0f;
    //float    dist = sqrtf(ds2); //Slower than the two lines above
    //assert(!isnan(dist));

    float epseffi = 1.f/epseff;
    float rhinv   = dist*epseffi;

    if(rhinv < 0.5f)
    {
      ids3  = 4.f/3.f + (rhinv*rhinv)*(4.f*rhinv-4.8f);
      ids   = 1.4f - (rhinv*rhinv)*(8.f/3.f+(rhinv*rhinv)*(3.2f*rhinv-4.8f));
    }
    else
    {
      ids3 = 8.f/3.f-6.f*rhinv+4.8f*(rhinv*rhinv)-4.f/3.f*(rhinv*rhinv*rhinv)-1.f/120.f/(rhinv*rhinv*rhinv);
      ids  = 1.6f-1/(30.f*rhinv)-(rhinv*rhinv)*(16.f/3.f+rhinv*(-8.f+rhinv*(4.8f-rhinv*16.f/15.f)));
    }//end if rhin < 0.5

    ids  = ids*2.f*epseffi;
    ids3 =  ids3*8.f*(epseffi*epseffi*epseffi);

    //Self gravity prevention (and NaN prevention)
    ids  *= selfGrav;
    ids3 *= selfGrav;

  } //end dist >= epseff
  
  ids2 = ids*ids;
  float ids5 = ids3*ids2;
  float ids7 = ids5*ids2;
  return (float4){ids, -ids3, +3.0f*ids5, -15.0f*ids7};
} 

#else

__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

#endif

#ifdef INDSOFT
__device__ float4 add_acc(float4 acc,  float4 pos,
                          float massj, float3 posj, float epsQ,
                          float &ds2, float eps2P, int selfGrav) {

#else

__device__ float4 add_acc(float4 acc,  float4 pos,
			  float massj, float3 posj,
			  float &ds2, float eps2, int selfGrav) {
#endif

  float3 dr = {pos.x - posj.x,
	       pos.y - posj.y,
	       pos.z - posj.z};

  #ifdef INDSOFT
    ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;
    float4 D04 = get_D04(ds2, eps2P, epsQ, selfGrav);
  #else
    ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z + eps2;
    float4 D04 = get_D04(ds2, selfGrav);
  #endif

  float  D0  = D04.x*massj;
  float  D1  = D04.y*massj;

  acc.w -= D0;
  acc.x += D1*dr.x;
  acc.y += D1*dr.y;
  acc.z += D1*dr.z;

  return acc;
}

#ifdef INDSOFT
__device__ float4 add_acc(float4 acc, float4 pos,
                          float mass, float3 com,
                          float3 Q0,  float3 Q1, float epsNode, float eps2) {
#else

__device__ float4 add_acc(float4 acc, float4 pos,
			  float mass, float3 com,
			  float3 Q0,  float3 Q1, float eps2) {
#endif
  float3 dr = {pos.x - com.x,
	       pos.y - com.y,
	       pos.z - com.z};

  #ifdef INDSOFT
    float  ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;
    float4 D04 = get_D04(ds2, epsNode, eps2);
  #else
    float  ds2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z + eps2;
    float4 D04 = get_D04(ds2);
  #endif

  float  D0  = D04.x*mass;
  float  D1  = D04.y*mass;
  float  D2  = D04.z*mass;
  float  D3  = D04.w*mass;

  float oct_q11 = Q0.x;
  float oct_q22 = Q0.y;
  float oct_q33 = Q0.z;
  float oct_q12 = Q1.x;
  float oct_q13 = Q1.y;
  float oct_q23 = Q1.z;

  float Qii = oct_q11 + oct_q22 + oct_q33;
  float QijRiRj =
          (oct_q11*dr.x*dr.x + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z) +
    2.0f*(oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z);
//2.0f was
  //volatile float QijRiRj_1 = oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z;
 //QijRiRj = dr.x*dr.x;// + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z;  //oct_q11;//*

  acc.w        -= D0 + 0.5f*D1*Qii + 0.5f*D2*QijRiRj;
  float C01a    = D1 + 0.5f*D2*Qii + 0.5f*D3*QijRiRj;
  acc.x         += C01a*dr.x + D2*(oct_q11*dr.x + oct_q12*dr.y + oct_q13*dr.z);
  acc.y         += C01a*dr.y + D2*(oct_q12*dr.x + oct_q22*dr.y + oct_q23*dr.z);
  acc.z         += C01a*dr.z + D2*(oct_q13*dr.x + oct_q23*dr.y + oct_q33*dr.z);

  return acc;
}


//Minimum distance opening criteria
#ifdef INDSOFT
  __device__ bool split_node_grav_md(float4 nodeCenter, float4 nodeSize, float4 groupCenter, float4 groupSize,
                                     float group_eps, float node_eps)
#else
  __device__ bool split_node_grav_md(float4 nodeCenter, float4 nodeSize, float4 groupCenter, float4 groupSize)
#endif
{
  //Compute the distance between the group and the cell
  float3 dr = {fabs(groupCenter.x - nodeCenter.x) - (groupSize.x + nodeSize.x),
               fabs(groupCenter.y - nodeCenter.y) - (groupSize.y + nodeSize.y),
               fabs(groupCenter.z - nodeCenter.z) - (groupSize.z + nodeSize.z)};

  dr.x += fabs(dr.x); dr.x *= 0.5f;
  dr.y += fabs(dr.y); dr.y *= 0.5f;
  dr.z += fabs(dr.z); dr.z *= 0.5f;

  //Distance squared, no need to do sqrt since opening criteria has been squared
  float ds2    = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;

  #ifdef INDSOFT
    //Naar idee van Inti nu minder overbodige openingen
    if(ds2 <= ((group_eps + node_eps ) * (group_eps + node_eps) )) return true;
  #endif

  return (ds2 <= fabs(nodeCenter.w));
}

//Improved Barnes Hut criterium
#ifdef INDSOFT
__device__ bool split_node_grav_impbh(float4 nodeCOM, float4 groupCenter, float4 groupSize,
                                      float group_eps, float node_eps)
#else
__device__ bool split_node_grav_impbh(float4 nodeCOM, float4 groupCenter, float4 groupSize)
#endif
{
  //Compute the distance between the group and the cell
  float3 dr = {fabs(groupCenter.x - nodeCOM.x) - (groupSize.x),
               fabs(groupCenter.y - nodeCOM.y) - (groupSize.y),
               fabs(groupCenter.z - nodeCOM.z) - (groupSize.z)};

  dr.x += fabs(dr.x); dr.x *= 0.5f;
  dr.y += fabs(dr.y); dr.y *= 0.5f;
  dr.z += fabs(dr.z); dr.z *= 0.5f;

  //Distance squared, no need to do sqrt since opening criteria has been squared
  float ds2    = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;

  #ifdef INDSOFT
    //Extra test
    if(ds2 <= ((group_eps + node_eps ) * (group_eps + node_eps) )) return true;
  #endif

  return (ds2 <= fabs(nodeCOM.w));
}



    }
  }
  __syncthreads();
#endif

  return acc_i;
}
