anchor.diff:      double precision, dimension(3) :: solvent_v
       integer :: idx, s2, i_r, i_th
   
  +    double precision, allocatable :: this_c(:,:,:)
  +    integer, allocatable :: this_count(:,:,:)
  +    double precision, allocatable :: this_v(:,:,:,:)
  +
       integer :: cell_idx, start, n
  +    integer :: n_species, n_theta, n_r
   
       r_min = this%r_min
       r_min_sq = r_min**2
  
positive.code:     dtheta = this%dtheta
      L = cells%edges
  
      call this%time_polar_update%tic()
      !$omp parallel do &
      !$omp private(cell_idx, idx, start, n, s2, d, r_sq, r, i_r, theta, i_th, one_r, &
      !$omp out_of_plane, one_theta, solvent_v, v_th, v_r)
      do cell_idx = 1, cells%N
  
         ! TODO: check if minimum distance from x to cell corners is below r_max
  
negative.code:              solvent_v = solvent%vel(:,idx) - v
               v_r = dot_product(solvent_v, one_r)
               v_th = dot_product(solvent_v, one_theta)
               this_v(:, s2, i_th, i_r) = this_v(:, s2, i_th, i_r) + [v_r, v_th]
               this_count(s2, i_th, i_r) = this_count(s2, i_th, i_r) + 1
            end if
         end do
      end do
  
      this%v = this%v + this_v
      this%c = this%c + this_c
      this%count = this%count + this_count
  
      deallocate(this_c)
      deallocate(this_v)
      deallocate(this_count)
  
      call this%time_polar_update%tac()
  
    end subroutine update
  