#!/usr/bin/env python

# Plot Humidex graph with Matplotlib

from pylab import *

figure(figsize = (10, 9))
tt = range(20, 51)

def humidex(t, r):
    "Calculate Humidex value for T (°C) and RH (%)"

    # https://en.wikipedia.org/wiki/Dew_point
    b = 17.67
    c = 243.5
    gamma = log(r / 100) + (b * t) / (c + t)
    t_dew = c * gamma / (b - gamma)

    # https://en.wikipedia.org/wiki/Humidex
    hh = 5417.753 * (1 / 273.16 - 1 / (273.15 + t_dew))
    h = t + .5555 * (6.11 * exp(hh) - 10)
    return round(h)

for r in range(100, 0, -10):
    plot(tt, [humidex(t, r) for t in tt], label = "RH = %u%%" % r,
        lw = 1.5, ls = "dashed", marker = "o")

legend()
fill_between(tt, 19.5, 29.5, color = "#0f0", alpha = .3)
fill_between(tt, 29.5, 39.5, color = "#ff0", alpha = .3)
fill_between(tt, 39.5, 45.5, color = "#f80", alpha = .3)
fill_between(tt, 45.5, 200,  color = "#f00", alpha = .3)

xticks(range(100))
yticks(range(150))
axis([20, 50, 20, 60])
xlabel("T [°C]")
ylabel("Humidex")
title("Humidex based on temperature and relative humidity")
grid()
savefig("humidex_plot_20-50.svg")

show()