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(double SG, double Volume, double Amount, double Time, double Alpha)
{
    double alpha = Alpha / 100.0;
    double mass = Amount * 1000.0;

    /*
     * 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_factor = ((1 - exp(-0.04 * Time)) / 4.15);
    double ibu = Bigness_factor * BoilTime_factor * AddedAlphaAcids;

    qDebug() << "boilIBU" << SG << Volume << Amount << Time << Alpha << "IBU:" << ibu;
    return ibu;
}


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 Fulltime, int Cooltype, double Coolparm1, double Coolparm2)
{
    double steep_time = 0;		/* Total time a hop in the kettle.				*/
    double loss_boiltemp = 1.0;		/* Loss due to the lower boil temperature at higher altitude.	*/

    double ibu = 0.0, whirlibus = 0.0;
    double alpha = Alpha / 100.0;
    double mass = Amount * 1000.0;

    // Ideas from Zymurgy March-April 2018. These are not exact formulas!
    if (Use == HOP_USEAT_AROMA) {
	if (Whirlpool9) {	// Flameout hops are 2 minutes in this range.
	    whirlibus += (alpha * mass * 20) / Volume * (2.0 / 50.0);
	}
	if (Whirlpool7) {	// Flameout hops are 4 minutes in this range.
	    whirlibus += (alpha * mass * 6) / Volume * (4.0 / 50.0);
	}
	// Experiment.
//	double wibu = boilIBU(Form, SG, Volume, Amount, 6, Alpha, Method);	// IBU's for 6 minutes full
//	double fibu = wibu * 0.067 * IBU_reduction(94);		// Add timed segments
//	fibu += wibu * 0.1 * IBU_reduction(84);
//	fibu += wibu * 0.167 * IBU_reduction(74);
//	fibu += wibu * 0.250 * IBU_reduction(64);
//	fibu += wibu * 0.417 * IBU_reduction(54);
	//qDebug() << "  94" << wibu * 0.067 * IBU_reduction(94);
	//qDebug() << "  84" << wibu * 0.1 * IBU_reduction(84);
	//qDebug() << "  74" << wibu * 0.167 * IBU_reduction(74);
	//qDebug() << "  64" << wibu * 0.250 * IBU_reduction(64);
	//qDebug() << "  54" << wibu * 0.417 * IBU_reduction(54);
//	qDebug() << "flamout" << wibu << fibu;
    }
    if (Use == HOP_USEAT_WHIRLPOOL) {	// Flameout or any whirlpool
	if (Whirlpool9) {
	    // 20 mg/l/50 min
	    whirlibus += (alpha * mass * 20) / Volume * (Whirlpool9 / 50.0);
	    //qDebug() << "Whirlpool9" << alpha * mass * 20 << " liter:" << liters << " time:" << Whirlpool9 << " ibu" << (alpha * mass * 20) / liters * (Whirlpool9 / 50.0);
	}
	if (Whirlpool7) {
	    // 6 mg/l/50 min
	    whirlibus += (alpha * mass * 6) / Volume * (Whirlpool7 / 50.0);
	    //qDebug() << "Whirlpool7" << alpha * mass * 6 << " liter:" << liters << " time:" << Whirlpool7 << " ibu" << (alpha * mass * 6) / liters * (Whirlpool7 / 50.0);
	}
	if (Whirlpool6) {
	    // 2 mg/l/50 min
	    whirlibus += (alpha * mass * 2) / Volume * (Whirlpool6 / 50.0);
	}
//	double wibu = boilIBU(Form, SG, Volume, Amount, Boiltime, Alpha, Method);
//	qDebug() << "whirpool" << wibu << wibu * IBU_reduction(74);
    }

    if ((Use == HOP_USEAT_MASH) || (Use == HOP_USEAT_FWH) || (Use == HOP_USEAT_BOIL)) {
	steep_time += Boiltime;
    }

    /*
     * IBU's from hops during Mash, FWH and boil.
     */
    if ((Use == HOP_USEAT_MASH) || (Use == HOP_USEAT_FWH) || (Use == HOP_USEAT_BOIL)) {
	ibu = TinsethIBU(SG, Volume, Amount, Fulltime, Alpha);
	/*
	 * Corrections for Mash and FWH
	 */
	if (Use == HOP_USEAT_MASH) {
	    ibu *= (1 + my_factor_mashhop / 100.0);
	}
	if (Use == HOP_USEAT_FWH) {
            ibu *= (1 + my_factor_fwh / 100.0);
	}

	/*
	 * 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);
	} else if (Form == HOP_FORMS_EXTRACT) {
	    // Nothing for now.
	}

	if (Method > 0) {
	    loss_boiltemp = IBU_reduction(boilPoint());
	    ibu *= loss_boiltemp;
	    qDebug() << "ibu" << ibu << "loss_boiltemp" << loss_boiltemp;
	}

//    } else if (Use == HOP_USEAT_AROMA) {
	/*
	 * At flameout. The cooling method is important.
	 * Emersion chiller, Counterflow chiller, Au bain marie or natural.
	 * Assume the hop is removed for all methods except Emersion chilling.
	 */
    } else {
//	qDebug() << "whirlibus" <<  whirlibus << Use;
    }

    return round((ibu + whirlibus) * 100.0) / 100.0;
}


double Utils::hopFlavourContribution(double bt, double vol, int use, double amount)
{
    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)
{
    double result = 0.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 * 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;
}


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