library(tidyverse)

rain_design <- 90.9
t_pattern <- c(4.1, 8.0, 11.0, 23.3, 15.3, 8.1, 7.0, 7.0, 5.1, 6.1, 3.1, 1.9)
t_step <-  0.5
t_step_count <- 1:length(t_pattern)
t_hour <- t_step_count * t_step
hyeto <- t_pattern * rain_design/100

hyeto_df <- data_frame(t_step_count, t_hour, rain = hyeto)
hyeto_df


# Chart the hyetograph

hyeto_df %>%
  ggplot(aes( x = t_step_count )) +
  geom_bar(aes(y = rain), stat = 'identity', fill = 'blue', width = 1, colour = 'black') +
  ylab("Rainfall (mm)") +
  scale_x_continuous(name = "Time (hours)",
                     breaks = c(t_step, t_step_count + t_step),
                     labels = c(0, t_hour)) + # some fiddling here to get the x-axis labelling correct
  theme_gray(base_size = 14)


# Function to remove initial loss from a hyetograph
#
# Arguments

# Hyeto = hyetograph (mm of rain for 1 or more time steps)
# IL = Initial Loss (mm)

Remove_IL <- function(Hyeto, IL){

  # IL <- 1
  # Hyeto <- c(30, 7, 4, 3, 7, 20, 15, 3)
  # #
  if(sum(Hyeto) < IL) {
    warning('Initial loss is not satisfied')
    Hyeto[] <- 0
    return(Hyeto)
  }

  rain_cum <- cumsum(Hyeto)

  tzero <- min(which(rain_cum > IL)) - 1

  if(tzero == 0) {
    Hyeto[1] <- Hyeto[1] - IL
    return(Hyeto)
  }

  Hyeto[tzero + 1] <- Hyeto[tzero + 1] - (IL - rain_cum[tzero])
  Hyeto[1:tzero] <- 0
  Hyeto

}

# Function to remove continuing loss from a hyetograph

# Arguments
# Hyetograph - rainfall for each time step
# CL - continuing loss in mm/h
# t_step - length of each time step in hours

Remove_CL <- function(Hyeto, CL, t_step) {

  CL_t_step = CL * t_step # Continuing loss per time step
  Hyeto <- Hyeto - CL_t_step
  Hyeto[Hyeto < 0] <- 0
  Hyeto


}


# Plot rainfall excess hyetograph

hyeto_df %>%
  mutate(rain_IL = Remove_IL(rain, 10),
         rain_ILCL = Remove_CL(rain_IL, 2, 0.5)) %>%
  ggplot(aes( x = t_step_count)) +
  geom_bar(aes(y = rain), stat = 'identity', fill = 'white', width = 1, colour = 'black') +
  geom_bar(aes(y = rain_IL), stat = 'identity', fill = 'yellow', width = 1, colour = 'black') +
  geom_bar(aes(y = rain_ILCL), stat = 'identity', fill = 'blue',  width = 1, colour = 'black') +
  ylab("Rainfall (mm)") +
  scale_x_continuous(name = "Time (hours)",
                     breaks = c(t_step, t_step_count + t_step),
                     labels = c(0, t_hour)) + # some fiddling here to get the x-axis labelling correct
  annotate(geom = 'rect', xmin = 8, xmax = 8.5, ymin = 20, ymax = 21, fill = 'white', colour = 'black') +
  annotate(geom = 'text', x = 8.6, y = 20.5, hjust = 0, label = 'Initial Loss' ) +
  annotate(geom = 'rect', xmin = 8, xmax = 8.5, ymin = 18.5, ymax = 19.5, fill = 'yellow', colour = 'black') +
  annotate(geom = 'text', x = 8.6, y = 19, hjust = 0, label = 'Continuing loss' ) +
  theme_gray(base_size = 14)

# Convert to a hydrograph

sub_area <- 78.7
t_step = 0.5 # time step in hours

# Example calculation

hyeto_df %>%
  mutate(rain_IL = Remove_IL(rain, 10),
         rain_ILCL = Remove_CL(rain_IL, 2, 0.5)) %>%
         select(rain_ILCL) %>%
         unlist(use.names = FALSE) # pull out the rainfall excess

# [1]  0.0000  0.0000  8.9990 20.1797 12.9077  6.3629  5.3630  5.3630  3.6359  4.5449  1.8179  0.7271

# Rainfall in the 3rd time step is 8.9 mm

(8.999 / 0.5) * 1e-3 * 78.7 * 1e+6 / 3600

# Calculate the rainfall excess hydrograph

# Add an extra row to the bottom of the data frame so the flow returns to zero
# and plot hydrograph

hyeto_df %>%
  mutate(rain_excess = Remove_CL(Remove_IL(rain, 10), 2, 0.5)) %>%
  add_row(t_step_count = 13,
          t_hour = 6.5,
          rain = 0,
          rain_excess = 0 ) %>%
  mutate(flow = (rain_excess/t_step) * sub_area * 1e3/3600) %>%
  ggplot(aes(x = t_hour, y = flow)) +
    geom_line(alpha = 0.5) +
    geom_point(shape = 21, color = 'grey', fill = 'blue', size = 3) +
    ylab(bquote('Flow ('*m^-3~s^-1*')')) +
    xlab('Time (hour)') +
    theme_gray(base_size = 14)



hyeto_df %>%
  mutate(rain_excess = Remove_CL(Remove_IL(rain, 10), 2, 0.5)) %>%
  add_row(t_step_count = 13,
          t_hour = 6.5,
          rain = 0,
          rain_excess = 0 ) %>%
  mutate(flow = (rain_excess/t_step) * sub_area * 1e3/3600)

# # A tibble: 13 × 5
# t_step_count t_hour    rain rain_excess      flow
# <dbl>  <dbl>   <dbl>       <dbl>     <dbl>
# 1             1    0.5  3.7269      0.0000   0.00000
# 2             2    1.0  7.2720      0.0000   0.00000
# 3             3    1.5  9.9990      8.9990 393.45628
# 4             4    2.0 21.1797     20.1797 882.30133
# 5             5    2.5 13.9077     12.9077 564.35333
# 6             6    3.0  7.3629      6.3629 278.20013
# 7             7    3.5  6.3630      5.3630 234.48228
# 8             8    4.0  6.3630      5.3630 234.48228
# 9             9    4.5  4.6359      3.6359 158.96963
# 10           10    5.0  5.5449      4.5449 198.71313
# 11           11    5.5  2.8179      1.8179  79.48263
# 12           12    6.0  1.7271      0.7271  31.79043
# 13           13    6.5  0.0000      0.0000   0.00000