library(ggplot2)
library(deSolve)

# Define the function that calculates the derivative of the water volume with respect to time
tank <- function(time, volume, params) {
  # Inputs:
  # time: current time
  # volume: current water volume in the tank
  # params: list of parameters (here, the flow rate and the outflow rate series)
  
  # Extract the flow rate and outflow rate series from the params list
  flow_rate <- params$flow_rate
  outflow_rate_series <- params$outflow_rate_series
  
  # Calculate the current outflow rate based on the time
  time_step <- which(time == params$time_interval)  # Get the current time step
  outflow_rate <- outflow_rate_series[time_step]  # Get the current outflow rate
  
  # Calculate the rate of change of the water volume based on the flow and outflow rates
  d_volume_dt <- -flow_rate + outflow_rate
  
  # Return the rate of change of the water volume as a list
  list(d_volume_dt)
}

# Set the initial water volume and flow rate
initial_volume <- 100  # litres
flow_rate <- 0.0  # litres per second
outflow_rate_series <- c(rep(c(rep(0.5, 12), rep(0, 12)), 5), rep(c(rep(0.5, 12), rep(0, 11)), 5))  # litres per second
length(outflow_rate_series)

# Set the time interval for the simulation
time_interval <- seq(from = 0, to = 123, by = 1)  # seconds

# Set the parameters for the simulation
params <- list(flow_rate = flow_rate, outflow_rate_series = outflow_rate_series, time_interval = time_interval)

# Run the simulation using the ode function
solution <- ode(y = initial_volume, times = time_interval, func = tank, parms = params)

head(solution)
# Create a data frame with the solution
solution_df <- data.frame(time = solution[, 1], volume = solution[, 2])

# Plot the solution using ggplot
p <- ggplot(data = solution_df, aes(x = time, y = outflow_rate)) +
  geom_line() +
  labs(x = "Time (s)", y = "Water volume (L)", title = "Water flow out of a tank")

# Save the plot as a PDF file
ggsave("water_flow_plot.pdf", plot = p)