# How much heat does a warm drink bring, compared with the body's own heat flows?
C_WATER = 4.186            # kJ kg⁻¹ K⁻¹ (a thin decoction is close to water)
C_BODY = 3.47              # kJ kg⁻¹ K⁻¹, mean specific heat of the human body
mass, core = 70.0, 37.0    # kg, °C

for ml, temp in ((200, 45), (200, 55), (300, 50)):
    q = ml / 1000 * C_WATER * (temp - core)                 # kJ given up as the drink reaches 37 °C
    rise = q / (mass * C_BODY)                              # if all of it stayed in the body
    minutes = q / (80 * 60 / 1000)                          # resting heat production ≈ 80 W
    print(f"{ml} mL at {temp} °C: {q:4.1f} kJ, mean body temperature +{rise:.3f} K at most, "
          f"= {minutes:3.1f} min of resting heat production")

# For comparison: a 40 °C bath that raises core temperature by 0.5 K stores about
print(f"\nbath raising the core by 0.5 K: ≈ {0.5 * mass * C_BODY * 0.7:.0f} kJ "
      f"(taking the warmed mass as 70 % of the body)")
# 200 mL at 45 °C:  6.7 kJ, mean body temperature +0.028 K at most, = 1.4 min of resting heat production
# 200 mL at 55 °C: 15.1 kJ, mean body temperature +0.062 K at most, = 3.1 min of resting heat production
# 300 mL at 50 °C: 16.3 kJ, mean body temperature +0.067 K at most, = 3.4 min of resting heat production
#
# bath raising the core by 0.5 K: ≈ 85 kJ (taking the warmed mass as 70 % of the body)

import math

def setpoint(h):
    """Circadian set point of core temperature: cosine, mesor 36.8 °C, amplitude 0.4 K, lowest at 05:00."""
    return 36.8 - 0.4 * math.cos(2 * math.pi * (h - 5) / 24)

def simulate(heat_kj=0.0, at_h=22.0, vasodilation=1.0, start=14.0, end=32.0, dt=1 / 60):
    """Two compartments, core and shell (kJ, hours). Heat flows core -> shell through blood flow,
    whose conductance rises when the core is above its set point (the controller); the shell
    loses heat to the bedding. 'vasodilation' multiplies the conductance after 22:00."""
    C_CORE, C_SHELL = 160.0, 80.0          # kJ/K
    M = 288.0                              # kJ/h, about 80 W
    T_ENV, H = 29.0, 52.0                  # bedding microclimate °C; shell -> environment kJ/h/K
    core, shell = setpoint(start), 34.4
    out, h = {}, start
    while h <= end + 1e-9:
        k = max(20.0, 75.0 + 900.0 * (core - setpoint(h)))     # kJ/h/K, bounded below
        if h >= 22.0:
            k *= vasodilation
        q_cs = k * (core - shell)
        core += dt * (M - q_cs) / C_CORE
        shell += dt * (q_cs - H * (shell - T_ENV)) / C_SHELL
        if abs(h - at_h) < dt / 2:
            core += heat_kj / C_CORE                              # a drink or a bath adds heat to the core
        out[round(h * 60)] = core
        h += dt
    return out

def report(label, run):
    at = lambda hh: run[round(hh * 60)]
    slope = (at(24.0) - at(22.5)) / 1.5
    print(f"{label:41s} 22:00 {at(22.0):.2f}, 00:00 {at(24.0):.2f} °C, slope {slope:+.3f} K/h")

print("core temperature; slope from 22:30 to 00:00")
report("nothing", simulate())
report("drink 200 mL at 45 °C, 22:00 (heat only)", simulate(heat_kj=6.7))
report("drink + 30 % more blood flow to the skin", simulate(heat_kj=6.7, vasodilation=1.30))
report("warm bath at 20:30 storing 85 kJ", simulate(heat_kj=85.0, at_h=20.5))
# core temperature; slope from 22:30 to 00:00
# nothing                                   22:00 36.97, 00:00 36.78 °C, slope -0.095 K/h
# drink 200 mL at 45 °C, 22:00 (heat only)  22:00 37.01, 00:00 36.78 °C, slope -0.096 K/h
# drink + 30 % more blood flow to the skin  22:00 37.01, 00:00 36.75 °C, slope -0.098 K/h
# warm bath at 20:30 storing 85 kJ          22:00 37.00, 00:00 36.79 °C, slope -0.103 K/h