Polynomial interpolation of ARF values between 12 and 24 hours. See https://tonyladson.wordpress.com/2020/05/23/smooth-interpolation-of-arf-curves/

ARF_poly_interp.R

# Smooth interpolation between long and short duration ARFs using a cubic.




library(tidyverse)
library(pracma)

devtools::source_url("https://gist.githubusercontent.com/TonyLadson/fc870cf7ebfe39ea3d1a812bcc53c8fb/raw/d8112631a92a32be749cabe334a22931c035711e/ARF2019.R?raw=TRUE")


#source(file.path('ARR2019_ARF', "ARF_2019.R"))
#source(file.path('ARR2019_ARF', "ARF_tests.R")) # Check that we pass tests


# Functions and constants ------------------------------------------------------



log_breaks = function(maj, radix=10) {
  function(x) {
    minx         = floor(min(logb(x,radix), na.rm=T)) - 1
    maxx         = ceiling(max(logb(x,radix), na.rm=T)) + 1
    n_major      = maxx - minx + 1
    major_breaks = seq(minx, maxx, by=1)
    if (maj) {
      breaks = major_breaks
    } else {
      steps = logb(1:(radix-1),radix)
      breaks = rep(steps, times=n_major) +
        rep(major_breaks, each=radix-1)
    }
    radix^breaks
  }
  
  
  
  
}
scale_x_log_eng = function(..., radix=10) {
  scale_x_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}
scale_y_log_eng = function(..., radix=10) {
  scale_y_continuous(...,
                     trans=scales::log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}








params <-
  structure(list(`East Coast North` =     c(0.327,  0.241, 0.448, 0.36,  0.00096,   0.48,  -0.21,   0.012,   -0.0013), 
                 `Semi-arid Inland QLD` = c(0.159,  0.283, 0.25,  0.308, 7.3e-07,   1,      0.039,  0,        0), 
                 Tasmania =               c(0.0605, 0.347, 0.2,   0.283, 0.00076,   0.347,  0.0877, 0.012,   -0.00033), 
                 `SW WA` =                c(0.183,  0.259, 0.271, 0.33,  3.845e-06, 0.41,   0.55,   0.00817, -0.00045), 
                 `Central NSW` =          c(0.265,  0.241, 0.505, 0.321, 0.00056,   0.414, -0.021,  0.015,   -0.00033), 
                 `SE Coast` =             c(0.06,   0.361, 0,     0.317, 8.11e-05,  0.651,  0,      0,        0), 
                 `Southern Semi-arid` =   c(0.254,  0.247, 0.403, 0.351, 0.0013,    0.302,  0.058,  0,        0), 
                 `Southern Temperate` =   c(0.158,  0.276, 0.372, 0.315, 0.000141,  0.41,   0.15,   0.01,    -0.0027), 
                 `Northern Coastal` =     c(0.326,  0.223, 0.442, 0.323, 0.0013,    0.58,  -0.374,  0.013,   -0.0015), 
                 `Inland Arid` =          c(0.297,  0.234, 0.449, 0.344, 0.00142,    0.216, 0.129,  0,        0)), 
            class = "data.frame", row.names = c("a", "b", "c", "d", "e", "f", "g", "h", "i"))



# Figure 1 from blog


AEP <- c(0.005)
area <- c(1000)
duration <- 10^seq(0,4,length.out = 500)
duration <- append(duration, c(719, 720, 722, 725, 727, 1440)) %>% sort()


ARF_df <- expand_grid(area = area, AEP = AEP, duration = duration)

ARF_short(area = 30000, duration = 720, aep = 0.005) #[1] 0.5151772


my_line_size <-  0.2
my_text_size <-  1.8
my_point_size <- 1


# log scale
p <- ARF_df %>% 
  rowwise() %>% 
  mutate(ARF_calc = ARF(area, duration, aep = AEP, region = 'Tasmania', params = params))  %>% 
  ggplot(aes(x = duration, y = ARF_calc, colour = factor(area))) +
  geom_line(size = my_line_size) + # size = 0.2
  annotate(geom = 'text', x = 1480, y = 0.25, angle = 90, label = '24 hours', size = my_text_size, hjust = 'left') +
  annotate(geom = 'text', x = 700, y = 0.25, angle = 90, label = '12 hours', size = my_text_size, hjust = 'left') +
  annotate(geom = 'point', x = 720, y = ARF_short(area = 1000, duration = 720, aep = 0.005), colour = scales::hue_pal()(3)[1], size = my_point_size) +
  annotate(geom = 'point', x = 1440, y = ARF_long(area = 1000, duration = 1440, aep = 0.005, params = params,region = 'Tasmania'), colour = scales::hue_pal()(3)[1], size = my_point_size) +
  geom_vline(xintercept = 720, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  geom_vline(xintercept = 1440, linetype = 'dotted', size = 0.2, alpha = 0.5) +
  scale_x_continuous(name = 'Duration (min)', breaks = c(200, 500, 720, 1000, 1440, 2000)) +
  #scale_x_log_eng(name = 'Duration (min)') +
  scale_y_continuous(name = 'Areal Reduction Factor') +
  scale_colour_hue(name = bquote('Area ('~km^2*')'))+
  
  labs(title = 'Tasmania',
       subtitle = bquote('AEP = 0.5%, Area = 1000'~km^2))+
  
  theme_grey(base_size = 6) + #5
  theme(legend.position = c(0.6, 0.3)) + 
  coord_cartesian(xlim = c(200, 2000), ylim = c(0.25,1)) +
  guides(colour = 'none') +
  NULL

p
ggsave(filename = file.path('figure','ARF_smooth_interp_3.png'), plot = p, width = 4, height = 3)








DARF_short <- function(A, D, P){
  
  a <- 0.287	
  b <- 0.265
  c <- 0.439
  d <- 0.36
  e <- 0.00226	
  f <- 0.226
  g <- 0.125
  h <- 0.0141
  i <- -0.021
  j <- 0.213
  
  
  a * d * D^(-d - 1) * (A^b - (c * log(D))/log(10)) + (a * c * D^(-d - 1))/log(10) + 
    e * g * A^f * D^(g - 1) * (log(P)/log(10) + 0.3) + 
    (1/9) * 2^((A * D * i)/1440 - 5) * 5^((A * D * i)/1440 - 1) * A *  h * i * log(10) * (log(P)/log(10) + 0.3)
  
  
  
  
  # (0.10332 * (A^0.265 - 0.190655 * log(D)))/(D^1.36) + 
  #   (0.0002825 *A^0.226 * (log(P)/log(10) + 0.3))/D^0.875 - 
  #   9.46938e-7*10^(-0.0000145833 * (D - 180)^2) * A^0.213 * (D - 180) * (log(P)/log(10) + 0.3) + 
  #   0.0547181/D^1.36
}  

# Plain text
# d/dD(1 - a (A^b - c log(10, D)) D^(-d) + e A^f D^g (0.3 + log(10, P)) + h×10^((j A D)/1440) (0.3 + log(10, P))) = a d D^(-d - 1) (A^b - (c log(D))/log(10)) + (a c D^(-d - 1))/log(10) + e g A^f D^(g - 1) (log(P)/log(10) + 0.3) + 1/9×2^((A D j)/1440 - 5)×5^((A D j)/1440 - 1) A h j log(10) (log(P)/log(10) + 0.3)


i

# Derivative of the long-duration ARF equation

DARF_long <- function(A, D, P, region = 'Tasmania', arf_params = params){ 
  
  
  for(z in 1:9) {
    assign(letters[z], arf_params[ ,region][z])
  }
  
  
  a * d * D^(-d - 1) * (A^b - (c * log(D))/log(10)) + 
    (a * c * D^(-d - 1))/log(10) + 
    e * g * A^f * D^(g - 1) * (log(P)/log(10) + 0.3) + 
    1/9 * 2^((A * D * i)/1440 - 5) * 5^((A * D * i)/1440 - 1) * A * h * i * log(10) * (log(P)/log(10) + 0.3)
  
}





# Values and slopes

ARF_long(area = 1000, duration = 1440, aep = 0.005, region = 'Tasmania', params) # [1] 0.8771158

# numerical derivative
pracma::fderiv(f = ARF_long, 1440, n = 1, method = 'central', area =1000, aep = 0.005, region = 'Tasmania', params = params) #[1] 2.01937e-05

# analytical derivative
DARF_long(A = 1000, D = 1440, P = 0.005) 
#[1] 2.019372e-05



ARF_short(area = 1000, duration = 720, aep = 0.005) # [1] [1] 0.8170708

# numerical derivative
pracma::fderiv(f = ARF_short, 720, n = 1, method = 'central', area =1000, aep = 0.005) #[1] [1] 6.579283e-05

# analytical derivative
DARF_short(A = 1000, D = 720, P = 0.005) 
# [1] 6.554367e-05


# Set up a matrix to solve for the cubic polynomial coefficients

x1 <- 720
x2 <- 1440
A <- 1000
P <-  0.005

B <- matrix(c(x1^3, x1^2, x1, 1,
  x2^3, x2^2, x2, 1,
  3*x1^2, 2*x1, 1, 0,
  3*x2^2, 2*x2, 1, 0), ncol = 4, byrow = TRUE)



F <- matrix(c(
  ARF_short(area = A, duration = 720, aep = P),
  ARF_long(area = A, duration = 1440, aep = P, region = 'Tasmania', params = params),
  DARF_short(A = A, D = 720, P = P),
  DARF_long(A = A, D = 1440, P = P, region = 'Tasmania', arf_params = params)), nrow = 4)

F

# Solve B %*% C = F for C where C is a matrix of coefficients

C <- as.vector(solve(B, F))




# Function describing the polynomial interpolation

Poly_interp <- function(D, coef){

  coef[1]*D^3 + coef[2]*D^2  + coef[3]*D + coef[4]
  
}

# Check the end points
# Value at 720 and 1440 mins

Poly_interp( D = 720, coef = C) # 0.8170708
Poly_interp( D = 1440, coef = C) # 0.8771158


# Add line to the graph


D_seq <- seq(from = 720, to = 1440, length.out = 100)

poly_interp_df <- tibble(duration = D_seq, ARF = Poly_interp(D_seq, coef = C))

# Add to graph

p

p2 <- p +
  geom_line(data = poly_interp_df, aes(x = duration, y = ARF), colour = 'blue') +
  coord_cartesian(xlim = c(200, 1800), ylim = c(0.75,0.9)) +
  annotate(geom = 'text', x = 1480, y = 0.75, angle = 90, label = '24 hours', size = my_text_size, hjust = 'left') +
  annotate(geom = 'text', x = 700, y = 0.75, angle = 90, label = '12 hours', size = my_text_size, hjust = 'left') 

p2


ggsave(filename = file.path('figure','ARF_smooth_interp_poly.png'), plot = p2, width = 4, height = 3)