/* heat index formula */
hi: c1 + c2*T + c3*h + c4*T*h + c5*T^2 + c6*h^2 + c7*T^2*h + c8*T*h^2 + c9*T^2*h^2;

/* parameter values */
vals: [c1=-42.379, c2=2.04901523, c3=10.1433127, c4=-0.22475541,
c5=-6.83783e-3, c6=-5.481717e-2, c7=1.22874e-3, c8=8.5282e-4,
c9=-1.99e-6]$

/* heat index calculation may be extended to the point where "heat index equals
temperature"
There are two solutions. */
[T1, T2]: solve([hi = T], [T]);

/* The limiting cases are where the discriminant becomes zero.
Pick the discriminant. */
disc: part(T2,2,1,1,1);

/* It's too complex calculating the humidity for a vanishing discriminant algebraically.
Calculate it numerically instead. There are two zeros. First the upper solution: */
hupper: find_root(ev(disc, vals) = 0, h, 70, 90);

/* Then the second, lower solution */
hlower: find_root(ev(disc, vals) = 0, h, 30, 50);

/* Pick the second, upper solution for "heat index eqals temperature" */
Tmatch: part(T2,2);

/* In the range (hlower, hupper) there is no solution to "heat index equals temperature".
Find the minimum difference instead. To calculate this minimum, take the first derivation of
hi wrt. T and find it's zero. */
Tnear: part(solve(diff(hi - T, T, 1) = 0, [T])[1],2);

/* Now there are three distict parts for the least temperature where the heat index may be
calculated: where "heat index equals temperature" for low humidities, where
"heat index comes closest to temperature" for medium humidities and where
"heat index equals temperature" for high humidities */
for plotval in [[Tmatch, 10, hlower], [Tnear, hlower, hupper], [Tmatch, hupper+0.00001, 100]] do (
[Tfunc, h0, h1]: plotval,
wxplot2d(ev(Tfunc, vals), [h, h0, h1], [ylabel, T]))$

/* heat index as a function of temperature for some humidities */
plotvals: [[30, Tmatch], [hlower, Tmatch], [(hupper+hlower)/2, Tnear], [hupper+0.0000001, Tmatch],
[90, Tmatch]]$

[Tmin, Tmax]: [65, 105]$

for plotval in plotvals do (
[h0, Tfunc]: plotval,
T0: ev(Tfunc, append(vals, [h=h0])),
T1: ev(hi,append(vals, [h=h0, T=T0])),
wxplot2d([ev(hi, append(vals, [h=h0])), T, [discrete, [[T0, T1]]]], [T, Tmin, Tmax],
[ylabel, HI], [style, lines, lines, points], [legend, h=h0, T, "Tlim"], [point_type, bullet]))$

limits: []$

/* calculate the limiting temperatures for all three regions,
collecting them in limits */
for hx from 10 step 1 thru hlower-1 do
(tl: ev(Tmatch, append(vals, [h=hx])),
limits: append(limits, [[hx, tl]]))$

for hx from floor(hlower) step 1 thru hupper do
(tl: ev(Tnear, append(vals), [h=hx]),
limits: append(limits, [[hx,tl]]))$

for hx from ceiling(hupper) step 1 thru 100 do
(tl: ev(Tmatch, append(vals, [h=hx])),
limits: append(limits, [[hx,tl]]))$

/* Plot the result */
wxplot2d([discrete, limits], [xlabel, h], [ylabel, T])$

/* generate a C header file declaring a static array that contains the calculated values */
heatindex_h: openw("heatindex.h")$

printf(heatindex_h, "/* file generated from wxMaxima~%")$

printf(heatindex_h, " *~% * this array contains limiting temperatures for wind chill calculations~%")$

[hmin, hmax]: [first(limits)[1], last(limits)[1]];

printf(heatindex_h, " * for humidities from ~d to ~d~% */~%~%", hmin, hmax)$

printf(heatindex_h, "static float heatindex[] = {~%")$

for pair in limits do (
[hx, Tx]: pair,
printf(heatindex_h, "    ~f,~%", Tx))$

printf(heatindex_h, "};~%")$

close(heatindex_h)$