anchor.diff:    double precision :: rho(r_l1:r_h1,r_l2:r_h2,r_l3:r_h3)
   
     integer          :: i, j, k
  -  double precision :: x, y, z
     double precision :: dV
  -  double precision :: r_P, r_S
  +  double precision :: r(3), r_P, r_S, grav_force_P, grav_force_S
   
     dV  = dx(1)*dx(2)*dx(3)
     volp = ZERO
     vols = ZERO
   
     do k = lo(3), hi(3)
  -     z = problo(3) + (dble(k) + HALF) * dx(3)
  +     r(3) = problo(3) + (dble(k) + HALF) * dx(3)
        do j = lo(2), hi(2)
  -        y = problo(2) + (dble(j) + HALF) * dx(2)
  +        r(2) = problo(2) + (dble(j) + HALF) * dx(2)
           do i = lo(1), hi(1)
  -           x = problo(1) + (dble(i) + HALF) * dx(1)
  +           r(1) = problo(1) + (dble(i) + HALF) * dx(1)
   
              if (rho(i,j,k) > rho_cutoff) then
   
  -              r_P = ( (x - com_p(1))**2 + (y - com_p(2))**2 + (z - com_p(3))**2 )**0.5
  -              r_S = ( (x - com_s(1))**2 + (y - com_s(2))**2 + (z - com_s(3))**2 )**0.5
  +              r_P = sqrt(sum((r - com_p)**2))
  +              r_S = sqrt(sum((r - com_s)**2))
  +
  +              grav_force_P = mass_P / r_P
  +              if (single_star .or. mass_S < 1.d-12 * mass_P) then
  +                 grav_force_S = ZERO
  +              else
  +                 grav_force_S = mass_S / r_S
  +              endif
                 
  -              if (r_P < roche_rad_P .or. single_star) then
  +              if (grav_force_P > grav_force_S) then
                    volp = volp + dV
  -              elseif (r_S < roche_rad_S) then
  +              else
                    vols = vols + dV
                 endif
                
  
positive.code: ! and given a density cutoff, computes the total volume of all zones
  ! whose density is greater or equal to that density cutoff.
  ! We also impose a distance requirement so that we only look 
  ! at zones that are within the effective Roche radius of the white dwarf.
  
  subroutine ca_volumeindensityboundary(rho,r_l1,r_l2,r_l3,r_h1,r_h2,r_h3,lo,hi,dx,&
                                        volp,vols,rho_cutoff)
  
    use bl_constants_module
    use probdata_module, only: mass_P, mass_S, com_P, com_S, roche_rad_P, roche_rad_S, single_star
    use prob_params_module, only: problo
  
    implicit none
  
negative.code: 
    use bl_constants_module, only: HALF, ZERO, ONE
    use prob_params_module, only: problo
    use probdata_module, only: mass_P, mass_S, com_P, com_S, single_star
  
    implicit none
  
  