src/Utils.cpp

/**
 * Utils.cpp is part of bmsapp.
 *
 * bmsapp is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * bmsapp is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include "Utils.h"
#include "MainWindow.h"
#include "global.h"

#include <QDebug>
#include <math.h>


double Utils::lintner_to_kolbach(double lintner)
{
    double wk = (3.5 * lintner) - 16;
    if (wk < 0)
	return 0.0;
    return wk;
}


double Utils::kolbach_to_lintner(double kolbach)
{
    return round(((kolbach + 16) / 3.5) * 1000.0) / 1000.0;
}


/**
 * Often used formulas divide or multiply with 1.97 to convert between EBC and SRM.
 * Almost all software in the world use this '1.97' formula.
 * The only alternative I have seen is "srm = (ebc * 0.375 + 0.46)", and that has
 * almost the same results as the formulas used in this program.
 * These formulas come from the Dutch 'brouwhulp' program written by Adrie Otten.
 *
 * Seems that the 'brouwhulp' version use the old EBC conversion (see below) and
 * not the new 1.97 formula. Make it selectable?
 */
double Utils::ebc_to_srm(double ebc)
{
    double srm = -0.00000000000132303 * pow(ebc, 4) - 0.00000000291515 * pow(ebc, 3) + 0.00000818515 * pow(ebc, 2) + 0.372038 * ebc + 0.596351;
    if (ebc < 0 || srm < 0)
	return 0;
    return srm;
}


double Utils::srm_to_ebc(double srm)
{
    double ebc = round( 0.000000000176506 * pow(srm, 4) + 0.000000154529 * pow(srm, 3) - 0.000159428 * pow(srm, 2) + 2.68837 * srm - 1.6004 );
    if ((ebc < 0) || (srm < 0))
	return 0;
    return ebc;
}


QString Utils::hours_to_string(int hours)
{
    int dd, hh;

    if (hours == 1)
	return QObject::tr("1 hour");
    if (hours < 24)
	return QString("%1 ").arg(hours) + QString(QObject::tr("hours"));

    dd = hours / 24;
    hh = hours % 24;
    if (dd == 1) {
	if (hh == 0)
	    return QString(QObject::tr("1 day"));
	else if (hh == 1)
	    return QString(QObject::tr("1 day, ")) + QString("%1 ").arg(hh) + QString(QObject::tr("hour"));
	else
	    return QString(QObject::tr("1 day, ")) + QString("%1 ").arg(hh) + QString(QObject::tr("hours"));
    } else {
	if (hh == 0)
	    return QString("%1 ").arg(dd) + QString(QObject::tr("days"));
	else if (hh == 1)
	    return QString("%1 ").arg(dd) + QString(QObject::tr("days, ")) + QString("%1 ").arg(hh) + QString(QObject::tr("hour"));
	else
	    return QString("%1 ").arg(dd) + QString(QObject::tr("days, ")) + QString("%1 ").arg(hh) + QString(QObject::tr("hours"));
    }
    return QString("hours_to_string error");
}


QColor Utils::srm_to_color(int srm)
{
    int		i;
    QColor	result;

    i = round(srm * 10);
    if (i < 0)
	i = 0;
    if (i > 299)
	i = 299;

    // A well known table for SRM to RGB conversion, range 0.1 to 30 SRM.
    const int R[] {
 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, //0
 250, 250, 250, 250, 250, 249, 248, 247, 246, 245, 244, 243, 242, 241, 240, 239, 238, 237, 236, 235, //2
 234, 233, 232, 231, 230, 229, 228, 227, 226, 225, 224, 223, 222, 221, 220, 219, 218, 217, 216, 215, //4
 214, 213, 212, 211, 210, 209, 208, 207, 206, 205, 204, 203, 202, 201, 200, 200, 199, 199, 198, 198, //6
 197, 197, 196, 196, 195, 195, 194, 194, 193, 193, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, //8
 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, 192, //10
 192, 192, 192, 192, 192, 192, 192, 192, 191, 190, 189, 188, 187, 186, 185, 184, 183, 182, 181, 180, //12
 179, 178, 177, 175, 174, 172, 171, 169, 168, 167, 195, 164, 162, 161, 159, 158, 157, 155, 154, 152, //14
 151, 149, 148, 147, 145, 144, 142, 141, 139, 138, 137, 135, 134, 132, 131, 129, 128, 127, 125, 124, //16
 122, 121, 119, 118, 117, 115, 114, 112, 111, 109, 108, 107, 105, 104, 102, 101, 99, 98, 97, 95, //18
 94, 92, 91, 89, 88, 87, 85, 84, 82, 81, 79, 78, 77, 75, 74, 72, 71, 69, 68, 67, //20
 65, 64, 62, 61, 59, 58, 57, 55, 54, 52, 51, 49, 48, 47, 45, 44, 43, 41, 39, 38, //22
 37, 37, 36, 36, 35, 35, 34, 34, 33, 33, 32, 32, 31, 31, 30, 30, 29, 29, 28, 28, //24
 27, 27, 26, 26, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 19, 18, 18, //26
 17, 17, 16, 16, 15, 15, 14, 14, 13, 13, 12, 12, 11, 11, 10, 10, 9, 9, 8, 8};

 const int G[] {
 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250, 250,
 250, 250, 250, 250, 250, 250, 249, 248, 247, 246, 245, 244, 242, 240, 238, 236, 234, 232, 230, 228,
 226, 224, 222, 220, 218, 216, 214, 212, 210, 208, 206, 204, 202, 200, 198, 196, 194, 192, 190, 188,
 186, 184, 182, 180, 178, 176, 174, 172, 170, 168, 166, 164, 162, 160, 158, 156, 154, 152, 150, 148,
 146, 144, 142, 141, 140, 139, 139, 138, 137, 136, 136, 135, 134, 133, 133, 132, 131, 130, 130, 129,
 128, 127, 127, 126, 125, 124, 124, 123, 122, 121, 121, 120, 119, 118, 118, 117, 116, 115, 115, 114,
 113, 112, 112, 111, 110, 109, 109, 108, 107, 106, 106, 105, 104, 103, 103, 102, 101, 100, 100, 99,
 98, 97, 97, 96, 95, 94, 94, 93, 92, 91, 91, 90, 89, 88, 88, 87, 86, 85, 85, 84,
 83, 82, 82, 81, 80, 79, 78, 77, 76, 75, 75, 74, 73, 72, 72, 71, 70, 69, 69, 68,
 67, 66, 66, 65, 64, 63, 63, 62, 61, 60, 60, 59, 58, 57, 57, 56, 55, 54, 54, 53,
 52, 51, 51, 50, 49, 48, 48, 47, 46, 45, 45, 44, 43, 42, 42, 41, 40, 39, 39, 38,
 37, 36, 36, 35, 34, 33, 33, 32, 31, 30, 30, 29, 28, 27, 27, 26, 25, 24, 24, 23,
 22, 22, 22, 21, 21, 21, 20, 20, 20, 19, 19, 19, 18, 18, 18, 17, 17, 17, 16, 16,
 16, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 12, 11, 11, 11, 10, 10, 10, 9,
 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3};

 const int B[] {
 210, 204, 199, 193, 188, 182, 177, 171, 166, 160, 155, 149, 144, 138, 133, 127, 122, 116, 111, 105,
 100, 94, 89, 83, 78, 72, 67, 61, 56, 50, 45, 45, 45, 46, 46, 46, 46, 47, 47, 47,
 47, 48, 48, 48, 48, 49, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 51, 52, 52, 52,
 52, 53, 53, 53, 53, 54, 54, 54, 54, 55, 55, 55, 55, 56, 56, 56, 56, 56, 56, 56,
 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56,
 56, 56, 56, 55, 55, 55, 55, 54, 54, 54, 54, 53, 53, 53, 53, 52, 52, 52, 52, 51,
 51, 51, 51, 50, 50, 50, 50, 49, 49, 48, 47, 47, 46, 45, 45, 44, 43, 43, 42, 41,
 41, 40, 39, 39, 38, 37, 37, 36, 35, 34, 33, 32, 31, 29, 28, 27, 26, 25, 24, 23,
 21, 20, 19, 18, 17, 16, 15, 13, 12, 11, 10, 9, 8, 9, 9, 10, 10, 11, 11, 12,
 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22,
 21, 21, 21, 20, 20, 20, 19, 19, 19, 18, 18, 18, 17, 17, 17, 17, 16, 16, 15, 15,
 15, 14, 14, 14, 13, 13, 13, 12, 12, 12, 11, 11, 11, 10, 10, 10, 9, 9, 9, 8,
 8, 8, 7, 7, 7, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2};

    result = QColor::fromRgb(R[i], G[i], B[i]);
    return result;
}


QColor Utils::ebc_to_color(int ebc)
{
    return srm_to_color(ebc_to_srm(ebc));
}


QString Utils::srm_to_style(int srm)
{
    QColor color = srm_to_color(srm);
    return QString("background-color: %1; color: %2;").arg(color.name()).arg((srm > 15) ? "#E0E1E3" : "#19232D");
}


QString Utils::ebc_to_style(int ebc)
{
    return srm_to_style(ebc_to_srm(ebc));
}


/*
 * Return incremented color by the boil and yeast.
 * https://www.hobbybrouwen.nl/forum/index.php/topic,19020.msg281132.html#msg281132
 */
double Utils::get_kt(int ebc)
{
    double kt = 1;

    if (ebc < 3)
	kt = 3.5;
    else if (ebc < 6)
	kt = 3;
    else if (ebc < 8)
	kt = 2.75;
    else if (ebc < 10)
	kt = 2.5;
    else if (ebc < 20)
	kt = 1.8;
    else if (ebc < 30)
	kt = 1.6;
    else if (ebc < 60)
	kt = 1.3;
    else if (ebc < 100)
	kt = 1.2;
    else if (ebc < 300)
	kt = 1.1;
    return kt;
}


double Utils::sg_to_plato(double sg)
{
    // On stChiellus: return ((135.997 * sg - 630.272) * sg + 1111.14) * sg - 616.868;
    return -668.962 + (1262.45 * sg) - (776.43 * sg * sg) + (182.94 * sg * sg * sg);
}


double Utils::plato_to_sg(double plato)
{
    return 1.00001 + (0.0038661 * plato) + (1.3488e-5 * plato * plato) + (4.3074e-8 * plato * plato * plato);
}


double Utils::sg_to_brix(double sg)
{
    return sg_to_plato(sg) * my_brix_correction;
}


double Utils::brix_to_sg(double brix)
{
    if (my_brix_correction > 0)
	return plato_to_sg(brix / my_brix_correction);
    return plato_to_sg(brix);
}


/*
 * @brief Calculate Final Gravity.
 *        Formula by Petr Novotny, Zymurgy July/August 2017.
 * @param o_plato Original Plato.
 * @param refracto the refractometer reading.
 * @return Final gravity.
 */
double Utils::brix_to_fg(double o_plato, double refracto)
{
    double FBc = refracto / my_brix_correction;

    double rc = round((1 + (0.006276 * FBc) - (0.002349 * o_plato)) * 10000.0) / 10000.0;
//    qDebug() << "brix_to_plato" << o_plato << refracto << FBc << "rc" << rc;
    return rc;
}


double Utils::calc_svg(double og, double fg)
{
    if (og == 0.0 || fg == 0.0)
	return 0;

    double oe = sg_to_plato(og);
    double ae = sg_to_plato(fg);

    return (oe - ae) / oe * 100;
}


double Utils::estimate_sg(double sugars, double batch_size)
{
    double plato = 100 * sugars / batch_size;
    double sg = plato_to_sg(plato);
    for (int i = 0; i < 20; i++) {
	if (sg > 0)
	    plato = 100 * sugars / (batch_size * sg);
	sg = plato_to_sg(plato);
    }

    return round(sg * 10000) / 10000;
}


/*
 * Returns the log base b of y.
 */
double Utils::logbase(double y, int b)
{
    double lg;

    lg = log10(y) / log10(b);
    return lg;
}


double Utils::estimate_fg(double psugar, double pcara, double wgratio, double mashtime, double mashtemp, double svg, double og, bool sta1)
{
    double BD;
    double att_mashtime;
    double att_mashtemp;

    if (psugar > 40)
	psugar = 0;
    if (pcara > 50)
	pcara = 0;

    if (wgratio > 0 && mashtime > 0) {
	BD = wgratio;
	if (BD < 2)
	    BD = 2;
	if (BD > 6)
	    BD = 6;
    } else {
	BD = 3.5;
	mashtemp = 67;
	mashtime = 75;
    }
    if (svg < 30)
	svg = 77;

    /*
     * Original from brouwhulp:
     *  0.00825 Attenuation factor yeast (Real to apparant ??)
     *  0.00817 Attenuation factor water/grain ratio
     * -0.00684 Attenuation factor mash temperature
     *  0.00026 Attenuation factor total mash time
     * -0.00356 Attenuation factor percentage crystal malt
     *  0.00553 Attenuation factor percentage simple sugars
     *  0.547   Attenuation factor constant
     *  0.597   Attenuation factor constant when STA1 gen is true.
     */
    double top_mashtemp = 66.11;	/* Highest fermentable at 151 degrees fahrenheit */
    /*
     * Derived from Karl Troester's data and Matt Kahn.
     * Original was just: -0.00684 * mashtemp
     */
    if (mashtemp > top_mashtemp) {
	/* Above optimum mash temperature decrease attenuation */
        att_mashtemp = -0.00684 * mashtemp;
    } else {
	/* Below optimum mash temperature decrease slowly attenuation */
        att_mashtemp = -0.00684 * (top_mashtemp + (top_mashtemp - mashtemp) / 4);
    }

    /*
     * Reference is set at 60 minutes equals old formula
     * Use log base 5 to create a good realistic response.
     * With mashtime < 12 minutes, the result will be negative, good.
     */
    att_mashtime = logbase(mashtime / 12, 5) * 0.0156;
    //qDebug() << "estimate_fg temp" << att_mashtemp << "time" << att_mashtime;

    double AttBeer = 0.00825 * svg + 0.00817 * BD + att_mashtemp + att_mashtime - 0.00356 * pcara + 0.00553 * psugar;
    AttBeer += (sta1) ? 0.597:0.547;
    qDebug() << "estimate_fg(" << psugar << pcara << BD << mashtime << mashtemp << svg << og << sta1 << ") AttBeer:" << AttBeer;
    return round((1 + (1 - AttBeer) * (og -1)) * 10000) / 10000;
}


/*
 * Kleurwerking to SRM. Not for Halberstadt, Naudts.
 */
double Utils::kw_to_srm(int colormethod, double c)
{
    if (colormethod == 0)
	return 1.4922 * pow(c, 0.6859);	//Morey
    if (colormethod == 1)
	return 0.3 * c + 4.7;		//Mosher
    if (colormethod == 2)
	return 0.2 * c + 8.4;		//Daniels
    return 0;				//Halberstadt,Naudts
}


double Utils::kw_to_ebc(int colormethod, double c)
{
    return srm_to_ebc(kw_to_srm(colormethod, c));
}


double Utils::kw_to_newebc(int colormethod, double c)
{
    return 1.97 * kw_to_srm(colormethod, c);
}


double Utils::abvol(double og, double fg)
{
    if (((og - fg) < 0) || (fg < 0.9))
	return 0;

    double factor = og * 3157 * pow(10, -5) + 9.716 * pow(10, -2);
    return round((og * 1000 - fg * 1000) * factor * 100) / 100;
}


double Utils::brewery_hPa()
{
    return Seapressure * exp( - MolMassAir * Gravacc * my_height / (Gasconst * (20 + Kelvin)));
}


double Utils::boilPoint()
{
    double P2 = brewery_hPa();

    return (1 / (1/(100 + Kelvin) - Gasconst * log(P2 / Seapressure) / EoVwater)) - Kelvin;
}


/*
 * Formula is from the 'Mash Made Easy' spreadsheet.
 * https://mashmadeeasy.yolasite.com/
 * https://www.homebrewtalk.com/threads/a-rather-simplified-whirlpool-hop-ibu-computation-method.701093/
 *
 * Source of the formula Mark G. Malowicki.
 */
double Utils::IBU_reduction(double Tc)
{
    /*
     * Original formula plus a small correction factor.
     */
    return 2.39 * pow(10, 11) * exp(-(9773 / (Tc + Kelvin))) * 1 / 1.009231743647;
}


double Utils::TinsethIBU(int Form, double SG, double Volume, double Amount, double T1, double T2, double Alpha, double Utilisation, double BU_factor)
{
    double alpha = Alpha / 100.0;
    double mass = Amount * 1000.0;

    /*
     * Basic Glenn Tinseth formula.
     * http://realbeer.com/hops/research.html
     *
     * kUtilisation was hardcoded as 4.15.
     * It is now a formula since we need variable Utilisation values per hop.
     * Default "normal" hops Utilisation is 20.
     */
    double kUtilisation = 83.0 / Utilisation;
    double AddedAlphaAcids = (alpha * mass * 1000) / Volume;
    double Bigness_factor = 1.65 * pow(0.000125, SG - 1);
    double BoilTime_factor1 = ((1 - exp(-0.04 * T1)) / kUtilisation);
    double BoilTime_factor2 = ((1 - exp(-0.04 * T2)) / kUtilisation);
    double ibu = Bigness_factor * (BoilTime_factor2 - BoilTime_factor1) * AddedAlphaAcids;

#ifdef DEBUG_IBU
    qDebug() << "TinsethIBU" << SG << Amount << Alpha << Volume << kUtilisation << Bigness_factor * (BoilTime_factor2 - BoilTime_factor1) << "ibu:" << ibu;
    qDebug() << "   boilIBU" << Form << SG << Volume << Amount << BoilTime_factor2 << BoilTime_factor1 << Alpha << "IBU:" << ibu;
#endif
    return ibu;
}


double Utils::toIBU(int Use, int Form, double preSG, double postSG, double Volume, double Amount, double Boiltime, double Alpha,
		    int Method, double Whirlpool9, double Whirlpool7, double Whirlpool6, double Fulltime,
		    int Cooltype, double CoolTo79, double CoolLPM, double Utilisation, double BU_factor)
{
    double ibu = 0.0;
    double bp = boilPoint();

    if (Use == HOP_USEAT_MASH) {
	/*
	 * Mash hops. About -30% to -90% utilization. Is a global setting.
	 * They count for 60 minutes, this is a fixed value.
	 * Almost all these hops will be gone after removing the malt. From
	 * pellets there may be some dust left, but that has minor effects.
	 *
	 * http://scottjanish.com/the-locksmith-utilizing-bioengineered-yeast-and-high-bound-thiol-precersour-hops-and-phantasm-powder-to-thiol-drive-beer/
	 */
	ibu = TinsethIBU(Form, preSG, Volume, Amount, 0, 60, Alpha, Utilisation, BU_factor) * (1 + my_factor_mashhop / 100.0);

    } else if ((Use == HOP_USEAT_FWH) || (Use == HOP_USEAT_BOIL)) {
	/*
	 * IBU's from hops during FWH and boil.
	 */
	double boil_time = Fulltime;
	if (Use == HOP_USEAT_BOIL)
	    boil_time = Boiltime;
	double fromSG = postSG - ((boil_time / Fulltime) * (postSG - preSG));	/* SG when this hop addition starts */
	double avgSG = (postSG + fromSG) / 2;					/* Average SG during this addition */
	ibu = TinsethIBU(Form, avgSG, Volume, Amount, 0, boil_time, Alpha, Utilisation, BU_factor);

	/*
	 * Correction for FWH
	 */
	if (Use == HOP_USEAT_FWH) {
            ibu *= (1 + my_factor_fwh / 100.0);
	}

	if (Method > 0) {
	    /*
	     * Extra IBU's during flameout, chilling, hopstands for hops added during boil or first wort.
	     */
	    double nibu = ibu;
	    if (Use == HOP_USEAT_BOIL)	/* Correct for brewery height above sealevel */
	    	nibu *= IBU_reduction(bp);

	    /*
	     * Flameout, currently fixed 1 minute.
	     */
	    double flameout_time = 1;
	    double fibu = TinsethIBU(Form, postSG, Volume, Amount, boil_time, boil_time + flameout_time, Alpha, Utilisation, BU_factor);
	    fibu *= IBU_reduction(bp - 0.5);	/* Boilpoint minus 0.5 degree */
#ifdef DEBUG_IBU
	    qDebug() << "during flameout" << fibu << bp;
#endif
	    nibu += fibu;

	    if ((Cooltype == CHILLER_TYPE_IMMERSION || Cooltype == CHILLER_TYPE_AUBAINMARIE || Cooltype == CHILLER_TYPE_NOCHILL) && CoolTo79 > 0) {
		/*
		 * Direct wort chilling, calculate the IBU's when chilling down to 79 degrees.
		 */
		double cibu = 0, tibu;
		for (int i = 1; i < int(CoolTo79); i++) {
		    tibu = TinsethIBU(Form, postSG, Volume, Amount, boil_time + flameout_time + i, boil_time + flameout_time + i + 1, Alpha, Utilisation, BU_factor);
		    tibu *= IBU_reduction((bp - 1) - (i * ((bp - 80) / CoolTo79)));
#ifdef DEBUG_IBU
		    qDebug() << "  chill" << i << "time" << boil_time + flameout_time + i << "temp" << (bp - 1) - (i * ((bp - 80) / CoolTo79)) << "ibu" << tibu;
#endif
		    cibu += tibu;
		}
#ifdef DEBUG_IBU
		qDebug() << "  immersion" << cibu;
#endif
		nibu += cibu;
	    }

	    if (Cooltype == CHILLER_TYPE_COUNTERFLOW && CoolLPM > 0) {
		/*
		 * Counterflow chilling, calculate the IBU's isomerized in the kettle during transfer the hot wort to the chiller.
		 */
		double cibu = 0, tibu, V, A;
		int steps = trunc(Volume / CoolLPM);
		for (int i = 1; i <= steps; i++) {
		    V = Volume - (i * CoolLPM);
		    A = Amount - (i * (Amount / steps));
		    tibu = TinsethIBU(Form, postSG, V, A, boil_time + flameout_time + i, boil_time + flameout_time + i + 1, Alpha, Utilisation, BU_factor);
		    /*
		     * Correction for the natural cooling of the wort in the kettle.
		     * Asume 0.1 degree per minute.
		     */
		    tibu *= IBU_reduction((bp - 1) - (i * 0.1));
#ifdef DEBUG_IBU
		    qDebug() << "  chill" << i << steps << "left" << V << A << "time" << boil_time + flameout_time + i << "ibu" << tibu;
#endif
		    cibu += tibu;
		}
#ifdef DEBUG_IBU
		qDebug() << "  counterflow" << cibu;
#endif
                nibu += cibu;
	    }

	    /*
	     * Hopstands, this boil hop adds some IBU's too.
	     */
	    if (Whirlpool9) {
		double wibu9 = TinsethIBU(Form, postSG, Volume, Amount, boil_time + flameout_time, boil_time + flameout_time + Whirlpool9, Alpha, Utilisation, BU_factor);
		wibu9 *= IBU_reduction(87.0);
#ifdef DEBUG_IBU
		qDebug() << "during whirlpool9" << wibu9;
#endif
		nibu += wibu9;
	    }
	    if (Whirlpool7) {
		double wibu7 = TinsethIBU(Form, postSG, Volume, Amount, boil_time + flameout_time, boil_time + flameout_time + Whirlpool7, Alpha, Utilisation, BU_factor);
		wibu7 *= IBU_reduction(74.0);
#ifdef DEBUG_IBU
		qDebug() << "during whirlpool7" << wibu7;
#endif
		nibu += wibu7;
	    }
	    if (Whirlpool6) {
                double wibu6 = TinsethIBU(Form, postSG, Volume, Amount, boil_time + flameout_time, boil_time + flameout_time + Whirlpool6, Alpha, Utilisation, BU_factor);
                wibu6 *= IBU_reduction(63.0);
#ifdef DEBUG_IBU
                qDebug() << "during whirlpool6" << wibu6;
#endif
                nibu += wibu6;
            }
#ifdef DEBUG_IBU
	    qDebug() << "Old IBU" << ibu << "New IBU" << nibu;
#endif
	    ibu = nibu;
	}

    } else if ((Use == HOP_USEAT_AROMA) && (Method > 0)) {
	/*
	 * At flameout, and only using extended calculation.
	 */
	double flameout_time = 1;
        ibu = TinsethIBU(Form, postSG, Volume, Amount, 0, flameout_time, Alpha, Utilisation, BU_factor);
        ibu *= IBU_reduction(bp - 0.5);

	if ((Cooltype == CHILLER_TYPE_IMMERSION || Cooltype == CHILLER_TYPE_AUBAINMARIE || Cooltype == CHILLER_TYPE_NOCHILL) && CoolTo79 > 0) {
	    /*
	     * Direct wort chilling, calculate the IBU's when chilling down to 79 degrees.
	     */
	    double cibu = 0, tibu;
	    for (int i = 1; i < int(CoolTo79); i++) {
		tibu = TinsethIBU(Form, postSG, Volume, Amount, flameout_time + i, flameout_time + i + 1, Alpha, Utilisation, BU_factor);
		tibu *= IBU_reduction((bp - 1) - (i * ((bp - 80) / CoolTo79)));
#ifdef DEBUG_IBU
		qDebug() << "  chill" << i << "time" << flameout_time + i << "temp" << (bp - 1) - (i * ((bp - 80) / CoolTo79)) << "ibu" << tibu;
#endif
		cibu += tibu;
	    }
#ifdef DEBUG_IBU
	    qDebug() << "  immersion" << cibu;
#endif
	    ibu += cibu;
	}

	if (Cooltype == CHILLER_TYPE_COUNTERFLOW && CoolLPM > 0) {
	    /*
	     * Counterflow chilling, calculate the IBU's isomerized in the kettle during transfer the hot wort to the chiller.
	     */
	    double cibu = 0, tibu, V, A;
	    int steps = trunc(Volume / CoolLPM);
	    for (int i = 1; i <= steps; i++) {
		V = Volume - (i * CoolLPM);
		A = Amount - (i * (Amount / steps));
		tibu = TinsethIBU(Form, postSG, V, A, flameout_time + i, flameout_time + i + 1, Alpha, Utilisation, BU_factor);
		/*
		 * Correction for the natural cooling of the wort in the kettle.
		 * Asume 0.1 degree per minute.
		 */
		tibu *= IBU_reduction((bp - 1) - (i * 0.1));
#ifdef DEBUG_IBU
		qDebug() << "  chill" << i << steps << "left" << V << A << "time" << flameout_time + i << "ibu" << tibu;
#endif
		cibu += tibu;
	    }
#ifdef DEBUG_IBU
	    qDebug() << "  counterflow" << cibu;
#endif
	    ibu += cibu;
	}

	/*
         * Hopstands, this flameout hop adds some IBU's too.
         */
        if (Whirlpool9) {
            double wibu9 = TinsethIBU(Form, postSG, Volume, Amount, flameout_time, flameout_time + Whirlpool9, Alpha, Utilisation, BU_factor);
            wibu9 *= IBU_reduction(87.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool9" << wibu9;
#endif
            ibu += wibu9;
        }
        if (Whirlpool7) {
            double wibu7 = TinsethIBU(Form, postSG, Volume, Amount, flameout_time, flameout_time + Whirlpool7, Alpha, Utilisation, BU_factor);
            wibu7 *= IBU_reduction(74.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool7" << wibu7;
#endif
            ibu += wibu7;
        }
        if (Whirlpool6) {
            double wibu6 = TinsethIBU(Form, postSG, Volume, Amount, flameout_time, flameout_time + Whirlpool6, Alpha, Utilisation, BU_factor);
            wibu6 *= IBU_reduction(63.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool6" << wibu6;
#endif
            ibu += wibu6;
        }

    } else if ((Use == HOP_USEAT_WHIRLPOOL) && (Method > 0)) {
	/*
         * Hopstands.
         */
        if (Whirlpool9) {
            double wibu9 = TinsethIBU(Form, postSG, Volume, Amount, 0, Whirlpool9, Alpha, Utilisation, BU_factor);
            wibu9 *= IBU_reduction(87.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool9" << wibu9;
#endif
            ibu = wibu9;
        }
        if (Whirlpool7) {
            double wibu7 = TinsethIBU(Form, postSG, Volume, Amount, 0, Whirlpool7, Alpha, Utilisation, BU_factor);
            wibu7 *= IBU_reduction(74.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool7" << wibu7;
#endif
            ibu = wibu7;
        }
        if (Whirlpool6) {
            double wibu6 = TinsethIBU(Form, postSG, Volume, Amount, 0, Whirlpool6, Alpha, Utilisation, BU_factor);
            wibu6 *= IBU_reduction(63.0);
#ifdef DEBUG_IBU
            qDebug() << "during whirlpool6" << wibu6;
#endif
            ibu = wibu6;
        }
    } else if (Use == HOP_USEAT_BOTTLING) {
	/*
	 * Isomerized hop extracts to use at bottling.
	 * Assume 10% volume is lost during fermentation and transfers.
	 */
	if (Form == HOP_FORMS_ISOEXTRACT) {
	    /*
	     * ISO-alpha-extract hops, Tetra etc.
	     */
	    double dosageHl = ((100.0 / Utilisation) * (100.0 / Alpha) * 0.001) / BU_factor;        // IBU per liter
	    ibu = (Amount * 1000.0) / (dosageHl * Volume * 0.9);
	    qDebug() << "ISO extract" << Amount << Alpha << Volume * 0.9 << Utilisation << BU_factor << "ibu:" << ibu;
	}
    }

    double rc = round(ibu * 1000.0) / 1000.0;

#ifdef DEBUG_IBU
    qDebug() << "toIBU" << Use << Form << preSG << postSG << Volume << Amount << Boiltime << Alpha << Method << Whirlpool9 << Whirlpool7 << Whirlpool6 << Fulltime << Cooltype << CoolTo79 << CoolLPM << Utilisation << BU_factor << "rc:" << rc;
#endif
    return rc;
}


double Utils::hopFlavourContribution(double bt, double vol, int use, double amount, int form)
{
    double result;

    if (use == HOP_USEAT_WHIRLPOOL || use == HOP_USEAT_DRY_HOP)
	return 0;
    if (use == HOP_USEAT_FWH) {
	result = 0.15;   // assume 15% flavourcontribution for fwh
    } else if (bt > 50) {
	result = 0.10;   // assume 10% flavourcontribution as a minimum
    } else {
	result = 15.25 / (6 * sqrt(2 * 3.1416)) * exp(-0.5 * pow((bt - 21.0) / 6.0, 2.0));
	if (result < 0.10)
	    result = 0.10;  // assume 10% flavourcontribution as a minimum
    }
    return (result * amount * 1000.0) / vol;
}


double Utils::hopAromaContribution(double bt, double vol, int use, double amount, int form)
{
    double result = 0.0;
    double factor = 1.0;

    if (form == HOP_FORMS_CRYO)
	factor = 2.0;

    if (use == HOP_USEAT_DRY_HOP) {
	result = 1.33;
    } else if (use == HOP_USEAT_WHIRLPOOL) {
	if (bt > 30)
	    bt = 30; // Max 30 minutes
	result = 0.62 * bt / 30.0;
    } else if (bt > 20) {
	result = 0.0;
    } else if (bt > 7.5) {
	result = 10.03 / (4 * sqrt(2 * 3.1416)) * exp(-0.5 * pow((bt - 7.5) / 4.0, 2.0));
    } else if (use == HOP_USEAT_BOIL) {  // Boil
	result = 1;
    } else if (use == HOP_USEAT_AROMA) {  // Aroma
	result = 1.2;
    }
    return (result * amount * factor * 1000.0) / vol;
}


double Utils::mix(double v1, double v2, double c1, double c2)
{
    if ((v1 + v2) > 0) {
        return ((v1 * c1) + (v2 * c2)) / (v1 + v2);
    }
    return 0;
}


double Utils::Hardness(double calcium, double magnesium)
{
    return 2.497 * calcium + 4.164 * magnesium;
}


double Utils::Bicarbonate(double total_alkalinity, double ph)
{
    return (total_alkalinity / (1 + 2*pow(10, ph - 10.33)) * MMHCO3 /*61.016*/ / (MMCaCO3 / 2) /*50.043*/);
}


double Utils::RA_CaCO3(double bicarbonate, double carbonate, double calcium, double magnesium)
{
    return ((bicarbonate / MMHCO3) + (2*carbonate / MMCO3) - (2*calcium / MMCa)/3.5 - (2*magnesium / MMMg)/7) * 50;
}


double Utils::ResidualAlkalinity(double total_alkalinity, double calcium, double magnesium)
{
    return total_alkalinity - (calcium / 1.4 + magnesium / 1.7);
}


double Utils::PartCO3(double pH)
{
    double H = pow(10.0, -pH);
    return 100.0 * Ka1 * Ka2 / (H * H + H * Ka1 + Ka1 * Ka2);
}


double Utils::PartHCO3(double pH)
{
    double H = pow(10.0, -pH);
    return 100.0 * Ka1 * H / (H * H + H * Ka1 + Ka1 * Ka2);
}


double Utils::Charge(double pH)
{
    return (-2.0 * PartCO3(pH) - PartHCO3(pH));
}


double Utils::CalcFrac(double TpH, double pK1, double pK2, double pK3)
{
    double r1d = pow(10.0, TpH - pK1);
    double r2d = pow(10.0, TpH - pK2);
    double r3d = pow(10.0, TpH - pK3);
    double dd = 1.0 / (1.0 + r1d + r1d * r2d + r1d * r2d * r3d);
    double f2d = r1d * dd;
    double f3d = r1d * r2d * dd;
    double f4d = r1d * r2d * r3d * dd;
    return f2d + 2.0 * f3d + 3.0 * f4d;
}


double Utils::kettle_cm(double volume, double kettle_volume, double kettle_height)
{
    if ((volume > 0) && (kettle_volume > 0) && (volume <= kettle_volume))
	return round(100 * ((1 - volume / kettle_volume) * kettle_height) * 10.0) / 10.0;
    return 0;
}


double Utils::kettle_vol(double cm, double kettle_volume, double kettle_height)
{
    if ((cm >= 0) && (kettle_volume > 0) && (cm <= (kettle_height * 100)))
	return round(((kettle_height - (cm / 100)) / kettle_height) * kettle_volume * 10.0) / 10.0;
    return 0;
}


double Utils::ResCO2(double T)
{
    double F = T * 1.8 + 32;
    return round((3.0378 - 0.050062 * F + 0.00026555 * F * F) * 1000000.0) / 1000000.0;
}


double Utils::CarbCO2toS(double CO2, double T, double SFactor)
{
    double sugar = round((SFactor * (CO2 - ResCO2(T)) * 4.014094) * 1000000.0) / 1000000.0;
    if (sugar < 0)
        sugar = 0;
    return sugar;
}


double Utils::GetPressureBar(double gl, double T)
{
    if (gl <= 0)
        return 0;

    double P = ((gl / 10.0) / exp(-10.73797 + (2617.25 / (T + 273.15)))) - 1.013;
    if (P < 0)
        P = 0;

    P = round(P * 100.0) / 100.0;
    //qDebug() << "  GetPressureBar(" << gl << "," << T << ") P:" << P;
    return P;
}


double Utils::HydroCorrection(double mg, double tr, double tc)
{
    double trf = tr * 1.8 + 32;
    double tcf = tc * 1.8 + 32;

    return mg * ((1.00130346 - 1.34722124E-4 * trf + 2.04052596E-6 * trf * trf - 2.32820948E-9 * trf * trf * trf) /
	         (1.00130346 - 1.34722124E-4 * tcf + 2.04052596E-6 * tcf * tcf - 2.32820948E-9 * tcf * tcf * tcf));
}