/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 0.8.5 ] */ /* [wxMaxima: input start ] */ /* 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; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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]$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* heat index calculation may be extended to the point where "heat index equals temperature" There are two solutions. */ [T1, T2]: solve([hi = T], [T]); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* The limiting cases are where the discriminant becomes zero. Pick the discriminant. */ disc: part(T2,2,1,1,1); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* Then the second, lower solution */ hlower: find_root(ev(disc, vals) = 0, h, 30, 50); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* Pick the second, upper solution for "heat index eqals temperature" */ Tmatch: part(T2,2); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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]))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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]]$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ [Tmin, Tmax]: [65, 105]$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ 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]))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ limits: []$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* 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]]))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ for hx from floor(hlower) step 1 thru hupper do (tl: ev(Tnear, append(vals), [h=hx]), limits: append(limits, [[hx,tl]]))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ for hx from ceiling(hupper) step 1 thru 100 do (tl: ev(Tmatch, append(vals, [h=hx])), limits: append(limits, [[hx,tl]]))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* Plot the result */ wxplot2d([discrete, limits], [xlabel, h], [ylabel, T])$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ /* generate a C header file declaring a static array that contains the calculated values */ heatindex_h: openw("heatindex.h")$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ printf(heatindex_h, "/* file generated from wxMaxima~%")$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ printf(heatindex_h, " *~% * this array contains limiting temperatures for wind chill calculations~%")$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ [hmin, hmax]: [first(limits)[1], last(limits)[1]]; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ printf(heatindex_h, " * for humidities from ~d to ~d~% */~%~%", hmin, hmax)$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ printf(heatindex_h, "static float heatindex[] = {~%")$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ for pair in limits do ( [hx, Tx]: pair, printf(heatindex_h, " ~f,~%", Tx))$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ printf(heatindex_h, "};~%")$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ close(heatindex_h)$ /* [wxMaxima: input end ] */ /* Maxima can't load/batch files which end with a comment! */ "Created with wxMaxima"$