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)
	qDebug() << "ebc_to_srm(" << ebc << ") =" << srm;
    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))
	qDebug() << "srm_to_ebc(" << srm << ") =" << ebc;
    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)
{
    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);
}


double Utils::calc_svg(double og, double fg)
{
    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::toIBU(int Use, int Form, double SG, double Volume, double Amount, double Boiltime, double Alpha,
		    int Method, double Whirlpool9, double Whirlpool7, double Whirlpool6)
{
    double	fmoment = 1.0, pfactor = 1.0, ibu = 0, boilfactor;
    double	sgfactor, AddedAlphaAcids, Bigness_factor, BoilTime_factor, utiisation;

    double gravity = SG;
    double liters = Volume;
    double alpha = Alpha / 100.0;
    double mass = Amount * 1000.0;
    double time = Boiltime;

    if ((Use == 3) || (Use == 4) || (Use == 5)) {	// Aroma, Whirlpool or Dry hop.
	fmoment = 0.0;
    } else if (Use == 0) { // Mash
	fmoment += my_factor_mashhop / 100.0;		// Brouwhulp
    } else if (Use == 1) { // First wort
	fmoment += my_factor_fwh / 100.0;			// Brouwhulp, Louis, Ozzie
    }

    if (Form == 0) {				// Pellet
	pfactor += my_factor_pellet / 100.0;
    } else if (Form == 1) {			// Plug
	pfactor += my_factor_plug / 100.0;
    } else if (Form == 3) {			// Wet leaf
	pfactor += my_factor_wethop / 100.0;	// From https://github.com/chrisgilmerproj/brewday/blob/master/brew/constants.py
    } else if (Form == 4) {			// Cryo hop
	pfactor += my_factor_cryohop / 100.0;
    }

    // Ideas from Zymurgy March-April 2018. These are not exact formulas!
    double whirlibus = 0.0;
    if (Use == 3 || Use == 4) { // Flameout or any whirlpool

	if (Whirlpool9) {
	    // 20 mg/l/50 min
	    whirlibus += (alpha * mass * 20) / liters * (Whirlpool9 / 50.0);
	    qDebug() << "Whirlpool9" << alpha * mass * 20 << " liter:" << liters << " time:" << Whirlpool9 << " ibu" << (alpha * mass * 20) / liters * (Whirlpool9 / 50.0);
	} else {
	    if (Use == 3) { // Flameout hops are 2 minutes in this range.
		whirlibus += (alpha * mass * 20) / liters * (2.0 / 50.0);
	    }
	}
	if (Whirlpool7) {
	    // 6 mg/l/50 min
	    whirlibus += (alpha * mass * 6) / liters * (Whirlpool7 / 50.0);
	    qDebug() << "Whirlpool7" << alpha * mass * 6 << " liter:" << liters << " time:" << Whirlpool7 << " ibu" << (alpha * mass * 6) / liters * (Whirlpool7 / 50.0);
	} else {
	    if (Use == 3) { // Flameout hops are 4 minutes in this range.
		whirlibus += (alpha * mass * 6) / liters * (4.0 / 50.0);
	    }
	}
	if (Whirlpool6) {
	    // 2 mg/l/50 min
	    whirlibus += (alpha * mass * 2) / liters * (Whirlpool6 / 50.0);
	    //console.log('Whirlpool6:' + alpha * mass * 2 + ' liter:' + liters + ' time:' + Whirlpool6 + ' ibu' + (alpha * mass * 2) / liters * (Whirlpool6 / 50));
	}
    }

    if (Method == 0) { // Tinseth
	/* http://realbeer.com/hops/research.html */
	AddedAlphaAcids = (alpha * mass * 1000) / liters;
	Bigness_factor = 1.65 * pow(0.000125, gravity - 1);
	BoilTime_factor = ((1 - exp(-0.04 * time)) / 4.15);
	utiisation = Bigness_factor * BoilTime_factor;
	ibu = round((utiisation * AddedAlphaAcids * fmoment * pfactor + whirlibus) * 100) / 100;
    }
    if (Method == 2) { // Daniels
	if (Form == 2) // Leaf
	    boilfactor = -(0.0041 * time * time) + (0.6162 * time) + 1.5779;
	else
	    boilfactor = -(0.0051 * time * time) + (0.7835 * time) + 1.9348;
	if (gravity < 1050)
	    sgfactor = 0;
	else
	    sgfactor = (gravity - 1050) / 200;
	ibu = round((fmoment * ((mass * (alpha * 100) * boilfactor * 0.1) / (liters * (1 + sgfactor))) + whirlibus) * 100) / 100;
    }
    if (Method == 1) { // Rager
	boilfactor = fmoment * 18.11 + 13.86 * tanh((time * 31.32) / 18.27);
	if (gravity < 1050)
	    sgfactor = 0;
	else
	    sgfactor = (gravity - 1050) / 200;
	ibu = round(((mass * (alpha * 100) * boilfactor * 0.1) / (liters * (1 + sgfactor)) + whirlibus) * 100) / 100;
    }

    return ibu;
}


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

    if (use == 4 || use == 5) // Whirlpool or Dry-hop
	return 0;
    if (use == 1) {   // First wort
	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)
{
    double result = 0.0;

    if (use == 5) {         // Dry hop
	result = 1.33;
    } else if (use == 4) { // 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 == 2) {  // Boil
	result = 1;
    } else if (use == 3) {  // Aroma
	result = 1.2;
    }
    return (result * amount * 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::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;
}