  	  node_oct0[tid] = make_float3(octopole0.x, octopole0.y, octopole0.z);
  	  node_oct1[tid] = make_float3(octopole1.x, octopole1.y, octopole1.z);
 
-          //float darkMatterMass = octopole1.w;
-          //float stellarMass = monopole.w - darkMatterMass;
-
-          #ifdef INDSOFT
-            float temp          = node_eps[tid]; //Backup value in the shmem into register
-            node_eps[tid]       = octopole0.w;
-          #endif
-
-	  __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
+
+#ifdef INDSOFT
+    float temp          = node_eps[tid]; //Backup value in the shmem into register
+    node_eps[tid]       = octopole0.w;
+#endif
+
+    __syncthreads();
 #pragma unroll
-	  for (int i = 0; i < DIMx; i++)
-          {
-            apprCount++;
+    for (int i = 0; i < DIMx; i++)
+    {
+      apprCount++;
 #ifdef DOGRAV
-            #ifdef INDSOFT
-              acc_i = add_acc(acc_i, pos_i,
-                              node_mon0[offs + i], node_mon1[offs + i],
-                              node_oct0[offs + i], node_oct1[offs + i], node_eps[offs + i], eps2);
-
-            #else
-              acc_i = add_acc(acc_i, pos_i,
-                              node_mon0[offs + i], node_mon1[offs + i],
-                              node_oct0[offs + i], node_oct1[offs + i], eps2);
-            #endif
-
-#endif
+#ifdef INDSOFT
+      acc_i = add_acc(acc_i, pos_i,
+          node_mon0[offs + i], node_mon1[offs + i],
+          node_oct0[offs + i], node_oct1[offs + i], node_eps[offs + i], eps2);
+
+#else
+      acc_i = add_acc(acc_i, pos_i,
+          node_mon0[offs + i], node_mon1[offs + i],
+          node_oct0[offs + i], node_oct1[offs + i], eps2);
+#endif
+
+#endif
+    }
+    __syncthreads();
+#ifdef INDSOFT
+    node_eps[tid] = temp; //Restore original value in shmem
+    __syncthreads();
+#endif
+  }
+  __syncthreads();
+#endif
+
+#if 1
+  /***********************************/
+  /******       DIRECT          ******/
+  /***********************************/
+
+  int *sh_body = &approx[DIM];
+
+  flag         = split && leaf && use_node;                                //flag = split + leaf + use_node
+  int  jbody   = node_data & BODYMASK;                                     //the first body in the leaf
+  int  nbody   = (((node_data & INVBMASK) >> LEAFBIT)+1) & BTEST(flag);    //number of bodies in the leaf masked with the flag
+
+  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
+#ifdef OLDPREFIX
+  calc_prefix<DIM2>(prefix0, tid, flag);
+#else
+  inclusive_scan_block<ADDOP<int>, int>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
+#endif
+
+  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
+
+  // step 2
+#ifdef OLDPREFIX
+  int n_bodies  = calc_prefix<DIM2>(prefix0, tid, nbody);
+#else
+  inclusive_scan_block<ADDOP<int>, int>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
+  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
+#endif
+
+  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
+
+  int nl_pre = 0;                                                     //Number of leaves that have already been processed
+
+#define NJMAX (DIM*4)
+  while (n_bodies > 0) {
+    int nb    = min(n_bodies, NJMAX - n_direct);                    //Make sure number of bides to be extracted does not exceed
+    //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 && (direct[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:
+#ifdef OLDPREFIX
+    int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);
+#else
+    int nl = inclusive_scan_array<ADDOP<int>, int>              // inclusive scan to compute number of leaves to process
+      (&body_list[n_direct], nb, tid);            // to make sure that there is enough shared memory for bodies
+#endif
+    nb = direct[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
+    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) {
+      n_direct -= DIM;
+
+
+      float4 posj  = body_pos[body_list[n_direct + tid]];
+      //             float4 posj  = tex1Dfetch(texBody, body_list[n_direct + tid]);
+      sh_mass[tid] = posj.w;
+      sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
+      sh_jid [tid] = body_list[n_direct + tid];
+
+#ifdef INDSOFT
+      float   temp = sh_eps [tid];        //Store the value from Shmem into a register
+      sh_eps [tid] = body_vel[body_list[n_direct + tid]].w;  //Load the softening
+#endif
+      __syncthreads();
+#pragma unroll
+      for (int j = 0; j < DIMx; j++)
+      {
+        int selfGrav = (body_i != sh_jid[offs + j]);
+        // 	      if (body_i != sh_jid[offs + j]) //If statement replaced by multiplication
+        {
+          float ds2 =  1.0e10f;
+          direCount++;
+#ifdef DOGRAV
+#ifdef INDSOFT
+          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], sh_eps[offs + j], ds2, eps2, selfGrav);
+#else
+          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], ds2, eps2, selfGrav);
+#endif
+#endif
+          ds2 += selfGrav*1.0e10f;
+          if (ds2 < ds2_min) {
+            ngb     = sh_jid[offs + j];
+            ds2_min = ds2;
           }
- 	  __syncthreads();
-          #ifdef INDSOFT
-             node_eps[tid] = temp; //Restore original value in shmem
-             __syncthreads();
-          #endif
- 	}
-	__syncthreads();
-#endif
-
-#if 1
-	/***********************************/
-	/******       DIRECT          ******/
-	/***********************************/
-
-        int *sh_body = &approx[DIM];
-
-        flag         = split && leaf && use_node;                                //flag = split + leaf + use_node
-        int  jbody   = node_data & BODYMASK;                                     //the first body in the leaf
-        int  nbody   = (((node_data & INVBMASK) >> LEAFBIT)+1) & BTEST(flag);    //number of bodies in the leaf masked with the flag
-
-        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
-        #ifdef OLDPREFIX
-          calc_prefix<DIM2>(prefix0, tid, flag);
-        #else
-          inclusive_scan_block<ADDOP<int>, int>(prefix0, (int)flag, tid);       // inclusive scan on flags to construct array
-        #endif
-
-        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
-
-        // step 2
-        #ifdef OLDPREFIX
-          int n_bodies  = calc_prefix<DIM2>(prefix0, tid, nbody);
-        #else
-          inclusive_scan_block<ADDOP<int>, int>(prefix0, nbody, tid);        // inclusive scan to compute memory offset for each body
-          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
-        #endif
-
-	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
-
-	int nl_pre = 0;                                                     //Number of leaves that have already been processed
-
-        #define NJMAX (DIM*4)
-	while (n_bodies > 0) {
-	  int nb    = min(n_bodies, NJMAX - n_direct);                    //Make sure number of bides to be extracted does not exceed
-                                                                          //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 && (direct[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:
-         #ifdef OLDPREFIX
-            int nl = calc_prefix<DIM2>(nb, &body_list[n_direct], tid);
-          #else
-            int nl = inclusive_scan_array<ADDOP<int>, int>              // inclusive scan to compute number of leaves to process
-                            (&body_list[n_direct], nb, tid);            // to make sure that there is enough shared memory for bodies
-          #endif
-	  nb = direct[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
-	  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) {
-	    n_direct -= DIM;
-
-
-	    float4 posj  = body_pos[body_list[n_direct + tid]];
-//             float4 posj  = tex1Dfetch(texBody, body_list[n_direct + tid]);
-	    sh_mass[tid] = posj.w;
-	    sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
-	    sh_jid [tid] = body_list[n_direct + tid];
-
-            #ifdef INDSOFT
-              float   temp = sh_eps [tid];        //Store the value from Shmem into a register
-              sh_eps [tid] = body_vel[body_list[n_direct + tid]].w;  //Load the softening
-            #endif
-	    __syncthreads();
-#pragma unroll
-	    for (int j = 0; j < DIMx; j++)
-            {
-              int selfGrav = (body_i != sh_jid[offs + j]);
-// 	      if (body_i != sh_jid[offs + j]) //If statement replaced by multiplication
-              {
-		float ds2 =  1.0e10f;
-                direCount++;
-#ifdef DOGRAV
-                #ifdef INDSOFT
-                  acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], sh_eps[offs + j], ds2, eps2, selfGrav);
-                #else
-                  acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], ds2, eps2, selfGrav);
-                #endif
-#endif
-                ds2 += selfGrav*1.0e10f;
-		if (ds2 < ds2_min) {
-		  ngb     = sh_jid[offs + j];
-		  ds2_min = ds2;
-		}
-	      }
-	    }
-            #ifdef INDSOFT
-              //Restore the shmem value after all threads have used it
-              __syncthreads();
-              sh_eps[tid] = temp;
-            #endif
-	    __syncthreads();
-	  }
-
-	}
-	direct[tid] = body_list[tid];
-	__syncthreads();
+        }
+      }
+#ifdef INDSOFT
+      //Restore the shmem value after all threads have used it
+      __syncthreads();
+      sh_eps[tid] = temp;
+#endif
+      __syncthreads();
+    }
+
+  }
+  direct[tid] = body_list[tid];
+  __syncthreads();
 #endif
       } //end lvl
 


      float4 posj = body_pos[direct[tid]];

  //     float4 posj  = tex1Dfetch(texBody, direct[tid]);
      sh_mass[tid] = posj.w;
      sh_pos [tid] = make_float3(posj.x, posj.y, posj.z);
      sh_jid [tid] = direct[tid];

      #ifdef INDSOFT
        sh_eps [tid] = body_vel[direct[tid]].w;
      #endif
    } else {
      sh_mass[tid] = 0.0f;
      sh_pos [tid] = make_float3(1.0e10f, 1.0e10f, 1.0e10f);
      sh_jid [tid] = -1;

      #ifdef INDSOFT
        sh_eps [tid] = 0.01f;
      #endif
    }

    __syncthreads();
  #pragma unroll
    for (int j = 0; j < DIMx; j++) {
      if ((sh_jid[offs + j] >= 0)) {
        int selfGrav = (body_i != sh_jid[offs + j]);
        float ds2 =  1.0e10f;
        direCount++;
  #ifdef DOGRAV
        #ifdef INDSOFT
          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], sh_eps[offs + j], ds2, eps2, selfGrav);
        #else
          acc_i = add_acc(acc_i, pos_i, sh_mass[offs + j], sh_pos[offs + j], ds2, eps2, selfGrav);
        #endif
  #endif
        ds2 += selfGrav*1.0e10f;
        if (ds2 < ds2_min) {
          ngb     = sh_jid[offs + j];

        if((n_stack1 - c_stack0) >= (LMEM_STACK_SIZE << SHIFT))
        {
          //We overwrote our current stack
          apprCount = -1; return acc_i;	 
        }
      }
      __syncthreads();


      /***
       **** --> copy nodes1 to nodes0: done by reassigning the pointers
       ***/
      n_nodes0    = n_nodes1;

      n_stack_pre = n_stack0;
