anchor.diff:         if ( (new_pos(2)<0) .or. (new_pos(2)>L(2)) .or. (new_pos(3)<0) .or. (new_pos(3)>L(3)) ) &
               call yzwall_collision(old_pos, old_vel, new_pos, new_vel, im, L, dt, bc, g)
   
  -       if (new_pos(1)<0) then
  -          new_pos(1) = new_pos(1) + L(1)
  -       else if (new_pos(1)>L(1)) then
  -          new_pos(1) = new_pos(1) - L(1)
  -       end if
  +       im = floor(new_pos/L)
  +       new_pos = new_pos - im*L
          particles%pos(:,i) = new_pos
          particles%vel(:,i) = new_vel
          particles%image(:,i) = particles%image(:,i) + im
  
positive.code:     else if (bc(coll_dim) == SPECULAR_BC) then
         v(coll_dim) = -v(coll_dim)
      else if (bc(coll_dim) == PERIODIC_BC) then
         if (x(coll_dim) <= 0) then
            x(coll_dim) = x(coll_dim) + L(coll_dim)
         else if (x(coll_dim) >= L(coll_dim)) then
            x(coll_dim) = x(coll_dim) - L(coll_dim)
         end if
      end if
  
      wall_loop: do while (.true.)
  
negative.code:        else if (bc(coll_dim) == SPECULAR_BC) then
            v(coll_dim) = -v(coll_dim)
         else if (bc(coll_dim) == PERIODIC_BC) then
            jump = floor(x/L)
            x(coll_dim) = x(coll_dim) - jump(coll_dim)*L(coll_dim)
         end if
      end do wall_loop
  
      t_collision = t_remainder
      x = x + v*t_collision + gvec*t_collision**2 / 2
      v = v + gvec*t_collision
  
    end subroutine yzwall_collision
  
  