Here are the current hops utilization formulae that I have specs for.
For all of these the same inputs are used:
"minutes" is a number of minutes this hop is in the boil
"gravity" is a decimal-point specific gravity (e.g. "1.050")
"aa" is the alpha-acid content of the hops as a decimal (e.g 7.5% is "0.075")
"volume" is a decimal in liters (e.g. 5 gallons is "18.927").
"g" is a weight of the hops in grams (e.g. 1 oz is "28.35")
First, Rager. It definies hops utilization as a function of the number of minutes. tanh is a standard hyperbolic tangent function.
Then it defines a gravity-correction factor (here called "ga"), which is (gravity-1.050)/.2 if the gravity is over 1.050, or 0 if the gravity is 1.050 or lower.
Finally, it calculates the IBUs using these numbers:
Code:
def utilization(minutes):
return (18.11+13.86*tanh((minutes-32.32)/18.27))*0.01
def ga(gravity):
if gravity>1.050:
return 0
return (gravity-1.050)/.2
def ibu(g, minutes, aa, volume, gravity):
"""g: grams of hops"""
return (g*utilization(minutes)*aa*1000)/(volume*(1+ga(gravity)))
Next, Tinseth.
Tinseth defines utilization as the product of a "bigness" factor (dependent on the gravity of the boil) and a "btf" factor (dependent on the length of the boil). It then defines IBUs as the utilization times the milligrams per liter of alpha acids.
Code:
# Note: pow(x, y) is "x to the y power".
# e here is the mathematical constant e, approximately 2.71828...
def bigness(gravity):
return 1.65*pow(0.000125, (gravity-1))
def btf(minutes):
return (1-pow(e, (-0.04 *minutes)))/4.15
def utilization(minutes, gravity):
return btf(minutes)*bigness(gravity)
def mgl_aa(g, aa, volume):
return (aa*g*1000)/volume
def ibu(g, minutes, aa, volume, gravity):
return utilization(minutes, gravity) * mgl_aa(g, aa, volume)
Finally, Garetz. Garetz is a bit trickier since it requires on a large table lookup for utilization. I'll list it here, but I use a formula to get close results instead:
Code:
def chart_utilization(minutes):
if minutes <= 0:
return 0
elif minutes <= 15:
return 2
elif minutes <= 20:
return 5
elif minutes <= 25:
return 8
elif minutes <= 30:
return 11
elif minutes <= 35:
return 14
elif minutes <= 40:
return 16
elif minutes <= 45:
return 18
elif minutes <= 50:
return 19
elif minutes <= 60:
return 20
elif minutes <= 70:
return 21
elif minutes <= 80:
return 22
else:
return 23
If you do a curve fit, you can create a mathematical function that gives a very close result to the table lookup without discontinuities:
Code:
co_a = 2.5251990477100890E+01
co_b = 2.2060400680981388E+01
co_c = 3.3792756276592284E+00
co_d = 5.0643078483394244E-01
def utilization(x):
return co_a * (1.0 - pow(1.0+pow(x/co_b, co_c), -1.0*co_d))
Whichever method you want, it just gets more complicated from there There are a number of small factors calculated along the way, and then the very broadest level of the algorithm is an iterative process where you make a guess at the final IBUs (called "d_ibu" for desired IBUs), run the equations, and then if the answer isn't close to the guess you try again with a d_ibu closer to the answer you got:
Code:
def calc_ca(final_volume, boil_volume, gravity, d_ibu, elevation=0):
cf = final_volume/boil_volume
boil_gravity = cf*(gravity - 1) + 1
gf = (boil_gravity - 1.050)/.2 + 1
hf = (cf*d_ibu)/260 + 1
tf = (elevation/550)*.02 + 1
return gf*hf*tf
def ibu(g, minutes, aa, volume, gravity, boil_volume, elevation):
aa = 100 * aa
d_ibu = 50
calc_ibu = 0
# abs(x) means "the absolute value of x"
while abs(d_ibu-calc_ibu)>0.1:
correction = (d_ibu - calc_ibu)/5
d_ibu -= correction
ca = calc_ca(volume, boil_volume, gravity, d_ibu, elevation)
calc_ibu = (utilization(minutes)*aa*g*0.1)/(volume*ca)
return calc_ibu
As you can see, Garetz also takes elevation and boil volume into account.
For any readers who are programmers, note that the above is all well-formed Python code. If you head it up with the following:
Code:
from math import tanh, e, pow, log
you should be able to include any of those algorithms directly in your python program.