anchor.diff:            i1 = links(1,i_link)
             i2 = links(2,i_link)
             d = distances(i_link)
  -          mass1 = p%mass(p%species(i1))
  -          mass2 = p%mass(p%species(i2))
  -          inv_mass = 1/mass1 + 1/mass2
  +          mass1 = i_mass(p%species(i1))
  +          mass2 = i_mass(p%species(i2))
  +          inv_mass = mass1 + mass2
   
             s = rel_pos(p% pos(:,i1), p% pos(:,i2), edges)
   
  
positive.code:        if (error < precision) exit rattle_loop
      end do rattle_loop
      call p%time_rattle_vel%tac()
  
      if (rattle_i==rattle_max) write(*,*) 'rattle_max reached in rattle_body_vel'
  
  
negative.code:           i1 = links(1,i_link)
            i2 = links(2,i_link)
            d = distances(i_link)
            mass1 = i_mass(p%species(i1))
            mass2 = i_mass(p%species(i2))
            inv_mass = mass1 + mass2
            s = rel_pos(p% pos(:,i1), p% pos(:,i2), edges)
            k = abs(dot_product(p%vel(:,i1)-p%vel(:,i2), s) / (d**2*inv_mass))
            if (k>error) error = k
  