How to increase from a single particle multiple in a simulation

Here is what I have so far for the single particle:

library(tidyverse)
library(magick)
library(gganimate)
library(stringi)
library(ggtext)

particle_mass <- 4.65 #(x 10^-26) 
## A constant that dictates the number of decimal places that the data will be rounded to. Having this set as a constant avoids having to manually adjust a number which appears multiples times throughout the code if you wanted to change the number of decimal places for any reason.
rnd <- 4 
## initial particle coordinates 
x_i <- round(runif(n = 1, min = 0,max = 1), rnd) 
y_i <- round(runif(n = 1, min = 0,max = 1), rnd)

x_velocity <- 400 ## A constant for the x velocity in m/s
y_velocity <- 200 ##A  constant for the y velocity in m/s
time_step <- 0.001 ## A constant for the duration of a single time step in seconds
frame_total <- 200 ## A constant for the number of frames that the animation will be split into (two additional frames are added on at the end)
rate_factor <- 35 ## A factor that controls the playback rate of the animation. Increase this number and the animation will be slowed 

## The x and y displacement components per time step with a factor to slow down the animation to a more easily perceptible speed.
x_jump <- (x_velocity/rate_factor)/frame_total 
y_jump <- (y_velocity/rate_factor)/frame_total 

## Vectors of equal length for x position, y position, and corresponding time steps are initialized.
x_pos <- round(seq(x_i, (x_i + (frame_total*x_jump))-x_jump, x_jump), rnd)
y_pos <- round(seq(y_i, (y_i + (frame_total*y_jump))-y_jump, y_jump), rnd) 
times <- seq(0, (frame_total*time_step)-time_step, time_step)

## A loop that adjusts the direction of the x velocity vector when x = 1 or x = 0 to simulate wall contact
for (i in 2:length(x_pos)) {
  if (x_pos[i-1] > 1 | x_pos[i-1] < 0 ){x_jump = x_jump*(-1)}
  x_pos[[i]] <- round(x_pos[[i-1]] + x_jump, rnd)
}

##  A loop that adjusts the direction of the y velocity vector when y = 1 or y = 0 to simulate wall contact
for (i in 2:length(y_pos)) {
  if (y_pos[i-1] > 1 | y_pos[i-1] < 0 ){y_jump = y_jump*(-1)}
  y_pos[[i]] <- round(y_pos[[i-1]] + y_jump, rnd)
}


bounce <- data.frame(times, x_pos, y_pos) ##time and coordinate position information
bounce$x_dif <- c(0, diff(bounce$x_pos))## change in x position between time steps
bounce$x_turn <- c(0, diff(sign(bounce$x_dif))) ## helps detect changes in direction

## impulse magnitude experienced by the right wall at each time step
bounce <- bounce %>% mutate(impulse_right = case_when((bounce$x_turn) == -2 ~ 2*(x_velocity*particle_mass),
                                                      (bounce$x_turn) != -2 ~ 0))

bounce$cume_impulse_right <- cumsum(bounce$impulse_right) #cumulative impulse on the right wall 
## static plots of the moving particle
bounce_plot <- ggplot(bounce) +
  geom_point(aes(x_pos, y_pos), color = "azure4", size = 1) +
  labs(title = "Single gas particle model") + 
  ylab("") + xlab("") +
  theme_classic() +
  ## axis lines and tick marks
  theme(axis.line = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(), axis.ticks.y = element_blank(),
        axis.text.y = element_blank()) +
  ## plot title aesthetics
  theme(plot.title = element_text(family = "Arial", face = "bold", hjust = 0.5, 
                                  size = 9, color = "gray38")) +
  ## The vertical and horizontal lines create a border around the 1 x 1 box that contains the particle
  geom_vline(xintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=0, color="snow3", size=2) + 
  geom_hline(yintercept=1, color="snow3", size=2) +
  geom_vline(xintercept=1, color="#CC6666", size=1) 

bounce_anim <- bounce_plot + transition_time(times)

## The desired number of frames and duration of the animation are set within the 'animate' function.
bounce_gif <- animate(bounce_anim, width = 3, height = 3, units = "in", res = 200,
                      nframes = frame_total, renderer = magick_renderer())

My plan would be to just create a new data frame for each new particle, and then add a new geom_point line under the ggplot function using that new data. I'm concerned that if I was to add 5, 10 or more particles that the code would end up being very long. Is there a way to add multiple particles to the simulation without having to write so much code? Any help would be greatly appreciated.

Solution

Here's an approach based on making a longer dataframe to hold the journey of multiple id's:

particle_mass <- 4.65
ids = 5
frames = 200
n = nrow(df)
df <- expand.grid(frame = 1:frames, id = 1:ids)
df$x_pos = runif(n)
df$y_pos = runif(n)
df$x_vel = runif(n, min = -0.1, max = 0.1) 
df$y_vel = runif(n, min = -0.1, max = 0.1) 
df$impulse_right = 0

for(i in 1:ids) {
  for(f in 2:frames) {
    r = f + (i - 1) * frames
    df$x_pos[r] = df$x_pos[r-1] + df$x_vel[r-1]
    df$x_vel[r] = df$x_vel[r-1] * ifelse(df$x_pos[r] < 0 | df$x_pos[r] > 1, -1, 1)
    df$impulse_right[r] = ifelse(df$x_vel[r] < df$x_vel[r-1], 
                                 df$x_vel[r-1] * particle_mass, 0)
    df$y_pos[r] = df$y_pos[r-1] + df$y_vel[r-1]
    df$y_vel[r] = df$y_vel[r-1] * ifelse(df$y_pos[r] < 0 | df$y_pos[r] > 1, -1, 1)
  }
}

ggplot(df, aes(x_pos, y_pos)) +
  geom_point(aes(color = id)) +
  transition_states(frame)