heat_index.R

library(tidyverse)

CELSIUS = FALSE

heat_index16term = function(T, RH) {
  retval = 16.923 + 1.85212 * 1e-1 * T + 5.37941 * RH - 1.00254 * 1e-1 * T *
    RH +
    9.41695 * 1e-3 * T ^ 2 + 7.28898 * 1e-3 * RH ^ 2 + 3.45372 * 1e-4 *
    T ^ 2 * RH - 8.14971 * 1e-4 * T * RH ^ 2 +
    1.02102 * 1e-5 * T ^ 2 * RH ^ 2 - 3.8646 * 1e-5 * T ^ 3 + 2.91583 *
    1e-5 * RH ^ 3 + 1.42721 * 1e-6 * T ^ 3 * RH +
    1.97483 * 1e-7 * T * RH ^ 3 - 2.18429 * 1e-8 * T ^ 3 * RH ^ 2 + 8.43296 *
    1e-10 * T ^ 2 * RH ^ 3 - 4.81975 * 1e-11 * T ^ 3 * RH ^ 3
}

heat_index8term = function(T, RH) {
  retval = -42.379 + 2.04901523 * T + 10.14333127 * RH - 0.22475541 * T *
    RH - 6.83783 * 1e-3 * T ^ 2 -
    5.481717 * 1e-2 * RH ^ 2 + 1.22874 * 1e-3 * T ^ 2 * RH + 8.5282 * 1e-4 *
    T * RH ^ 2 - 1.99 * 1e-6 * T ^ 2 * RH ^ 2
}

heuristic = function(T, RH) {
  if (RH == 90) {
    retval = 5 + 3 * (T - 80) + T
  } else if (RH == 70) {
    retval = 2 * (T - 80) + T
  } else if (RH == 50) {
    retval = 1 * (T - 80) + T
  } else
    (stop("no match in heuristic"))
  retval
}

vheuristic = Vectorize(heuristic)

toC = function(F) {
  (F - 32) * 5 / 9
}

df = expand.grid(T = seq(80, 100, by = 2.5), RH = c(50, 70, 90)) %>%
  mutate(
    Humidity = factor(RH, levels = c(90, 70, 50)),
    hi16 = round(heat_index16term(T, RH), 0),
    hi8 = round(heat_index8term(T, RH), 0),
    blend = (hi16 + hi8) / 2,
    heur = vheuristic(T, RH),
    ad = abs(blend - heur),
    line = -55 / 18 * T + 3430 / 9
  ) %>%
  filter(RH <= line)

fm = lm(data = df, blend - T ~ (T + RH) * Humidity)

df = df %>%
  mutate(linfit = predict(fm),
         ad2 = abs(T + linfit - blend))

if (CELSIUS) {
  df = df %>% mutate(
    T = toC(T),
    hi16 = toC(hi16),
    hi8 = toC(hi8),
    blend = toC(blend),
    heur = toC(heur),
    linfit = linfit * 5 / 9
  )
  PRE = "C"
} else {
  PRE = "F"
}

p = ggplot(df, aes(
  x = T,
  y = hi16 - T,
  group = Humidity,
  color = Humidity
))
p = p + geom_line(size = 1, linetype = "dashed") +
  theme_bw() +
  labs(x = "Temperature", y = "Amount to add to temperature\nto get heat index") +
  theme(legend.position = "bottom")
p
ggsave(
  plot = p,
  file = paste(PRE, "hi16.png"),
  height = 4,
  width = 6
)

p = p + geom_line(aes(y = hi8 - T), size = 1, linetype = "dotdash")
p
ggsave(
  plot = p,
  file = paste(PRE, "hi8hi16.png"),
  height = 4,
  width = 6
)

p = p + geom_line(aes(y = blend - T), size = 2)
p
ggsave(
  plot = p,
  file = paste(PRE, "blendhi8hi16.png"),
  height = 4,
  width = 6
)

p = p + geom_line(data = df, aes(y = linfit))
p
ggsave(
  plot = p,
  file = paste(PRE, "linblendhi8hi16.png"),
  height = 4,
  width = 6
)

p = p + geom_point(aes(y = heur - T), size = 2)
p
ggsave(
  plot = p,
  file = paste(PRE, "heurlinblendhi8hi16.png"),
  height = 4,
  width = 6
)

p = ggplot(df, aes(
  x = T,
  y = blend - T,
  group = Humidity,
  color = Humidity
))
p = p + geom_line(aes(y = linfit), size = 2) +
  geom_point(aes(y = heur - T), size = 2) +
  theme_bw() +
  labs(x = "Temperature", y = "Amount to add to temperature to get heat index") +
  theme(legend.position = "bottom")
p
ggsave(
  file = paste(PRE, "heurlin.png"),
  height = 4,
  width = 6
)