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;
}


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

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

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

    /*
     * From brouwhulp:
     * 0.00825 Attenuation factor yeast
     * 0.00817 Attenuation factor water/grain ratio
     * -0.00684 Attenuation factor mash temperature
     * 0.00026 Attenuation factor total mash time  (at some places this is 0.0026 this is wrong!)
     * -0.00356 Attenuation factor percentage crystal malt
     * 0.00553 Attenuation factor percentage simple sugars
     * 0.547 Attenuation factor constant
     */
    double AttBeer = 0.00825 * svg + 0.00817 * BD - 0.00684 * mashtemp + 0.00026 * mashtime - 0.00356 * pcara + 0.00553 * psugar + 0.547;
    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;
    double ibu = 0.0;

    /*
     * Most suppliers of CO2 hop extract (not isomerized) use formulas
     * like this. For now, manual set a utilisation factor of 35%.
     */
    if (Form == HOP_FORMS_CO2EXTRACT) {
	double utilisation = 0.35;	// 35%
	double factor_sg = 1.0;
	double factor_bt = 1.0;
	double bt = T2 - T1;

	/* Table from Wildabouthops */
	if (SG >= 1.150)
	    factor_sg = 1.3;
	else if (SG >= 1.110)
	    factor_sg = 1.2;
	else if (SG >= 1.080)
	    factor_sg = 1.1;

	if (bt > 60) {
	    factor_bt = ((bt - 60) / 300) + 1.0;
	} else {
	    factor_bt = (bt / 60.0);
	}
	qDebug() << "factor_sg" << factor_sg << "factor_bt" << factor_bt;

	ibu = ((Utilisation / 100.0) * alpha * mass * 1000.0) / Volume;
	ibu = (ibu * factor_bt) / factor_sg;
	qDebug() << "TinsethIBU CO2" << Amount << Alpha << Volume << Utilisation << BU_factor << "SG" << SG << "T" << T1 << T2 << "ibu:" << ibu;
	return ibu;
    }

    if (Form == HOP_FORMS_ISOEXTRACT) {

	double dosageHl = ((100.0 / Utilisation) * (100.0 / Alpha) * 0.001) / BU_factor;	// IBU per liter
	qDebug() << (100 / Utilisation) << (100 / Alpha) << dosageHl << dosageHl * Volume << mass / (dosageHl * Volume);
	ibu = mass / (dosageHl * Volume);
	qDebug() << "TinsethIBU ISO" << Amount << Alpha << Volume << Utilisation << BU_factor << "ibu:" << ibu;
	return ibu;
    }

    /*
     * Basic Tinseth formula.
     * http://realbeer.com/hops/research.html
     */
    double AddedAlphaAcids = (alpha * mass * 1000) / Volume;
    double Bigness_factor = 1.65 * pow(0.000125, SG - 1);
    double BoilTime_factor1 = ((1 - exp(-0.04 * T1)) / 4.15);
    double BoilTime_factor2 = ((1 - exp(-0.04 * T2)) / 4.15);
    ibu = Bigness_factor * (BoilTime_factor2 - BoilTime_factor1) * AddedAlphaAcids;
    qDebug() << "TinsethIBU nor" << SG << Amount << Alpha << Volume << Bigness_factor * (BoilTime_factor2 - BoilTime_factor1) << "ibu:" << ibu;

    /*
     * Correction for hop forms
     */
    if (Form == HOP_FORMS_PELLET) {
	ibu *= (1 + my_factor_pellet / 100.0);
    } else if (Form == HOP_FORMS_PLUG) {
	ibu *= (1 + my_factor_plug / 100.0);
    } else if (Form == HOP_FORMS_LEAF_WET) {
	ibu *= (1 + my_factor_wethop / 100.0);      // From https://github.com/chrisgilmerproj/brewday/blob/master/brew/constants.py
    } else if (Form == HOP_FORMS_CRYO) {
	ibu *= (1 + my_factor_cryohop / 100.0);
    }

    qDebug() << "boilIBU" << Form << SG << Volume << Amount << BoilTime_factor2 << BoilTime_factor1 << Alpha << "IBU:" << ibu;
    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 Coolparm1, double Coolparm2, double Utilisation, double BU_factor)
{
    double ibu = 0.0;

    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) {
	    double nibu = ibu;
	    nibu *= IBU_reduction(boilPoint());

	    /*
	     * 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(98.0);
	    //qDebug() << "during flameout" << fibu;
	    nibu += fibu;

	    // Add this hop during cooling
	    /*
	     * 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);
		//qDebug() << "during whirlpool9" << wibu9;
		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);
		//qDebug() << "during whirlpool7" << wibu7;
		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);
                //qDebug() << "during whirlpool6" << wibu6;
                nibu += wibu6;
            }
	    //qDebug() << "Old IBU" << ibu << "New IBU" << nibu;
	    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(98.0);
	/*
         * 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);
            //qDebug() << "during whirlpool9" << wibu9;
            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);
            //qDebug() << "during whirlpool7" << wibu7;
            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);
            //qDebug() << "during whirlpool6" << wibu6;
            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);
            //qDebug() << "during whirlpool9" << wibu9;
            ibu = wibu9;
        }
        if (Whirlpool7) {
            double wibu7 = TinsethIBU(Form, postSG, Volume, Amount, 0, Whirlpool7, Alpha, Utilisation, BU_factor);
            wibu7 *= IBU_reduction(74.0);
            //qDebug() << "during whirlpool7" << wibu7;
            ibu = wibu7;
        }
        if (Whirlpool6) {
            double wibu6 = TinsethIBU(Form, postSG, Volume, Amount, 0, Whirlpool6, Alpha, Utilisation, BU_factor);
            wibu6 *= IBU_reduction(63.0);
            //qDebug() << "during whirlpool6" << wibu6;
            ibu = wibu6;
        }
    } else if (Use == HOP_USEAT_BOTTLING) {
	/*
	 * Isomerized hop extracts to use at bottling.
	 * Assume 10% volume is lost during fermentation and transfers.
	 */
	ibu = TinsethIBU(Form, postSG, Volume * 0.9, Amount, 0, 0, Alpha, Utilisation, BU_factor);
    }

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

    qDebug() << "toIBU" << Use << Form << preSG << postSG << Volume << Amount << Boiltime << Alpha << Method << Whirlpool9 << Whirlpool7 << Whirlpool6 << Fulltime << Cooltype << Coolparm1 << Coolparm2 << Utilisation << BU_factor << "rc:" << rc;
    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;
}