/*
 * Copyright (c) 2003 by Domenick Venezia 
 * 
 *  This software is distributed under the so-called ``Berkeley 
 *  License'': 
 * 
 * Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions 
 * are met: 
 * 
 * 1. Redistributions of source code must retain the above copyright 
 *    notice, this list of conditions and the following disclaimer. 
 * 2. Redistributions in binary form must reproduce the above copyright 
 *    notice, this list of conditions and the following disclaimer in the 
 *    documentation and/or other materials provided with the distribution. 
 * 
 * This software is provided ``as is'', and any express or implied 
 * warranties, including, but not limited to, the implied warranties of 
 * merchantability and fitness for a particular purpose are disclaimed. 
 * In no event shall Domenick Venezia be liable for any 
 * direct, indirect, incidental, special, exemplary, or consequential 
 * damages (including, but not limited to, procurement of substitute 
 * goods or services; loss of use, data, or profits; or business 
 * interruption) however caused and on any theory of liability, whether 
 * in contract, strict liability, or tort (including negligence or 
 * otherwise) arising in any way out of the use of this software, even if 
 * advised of the possibility of such damage. 
 *
 * Beer library functions.
 * Adapted by Greg Lehey
 *
 * $Id: beer.c,v 1.4 2004/01/28 02:41:40 grog Exp $
 */

#include "beer.h"

float density_water[BEER_NUM_DENS] =
{
  0.999868,			/*   0 C */
  0.999927,
  0.999968,
  0.999992,
  1.000000,
  0.999992,
  0.999968,
  0.999930,
  0.999877,
  0.999809,
  0.999728,			/*  10 C */
  0.999634,
  0.999526,
  0.999406,
  0.999273,
  0.999129,
  0.998972,
  0.998804,
  0.998625,
  0.998435,
  0.998234,			/*  20 C */
  0.998022,
  0.997801,
  0.997569,
  0.997327,
  0.997075,
  0.996814,
  0.996544,
  0.996264,
  0.995976,
  0.995678,			/*  30 C */
  0.995372,
  0.995057,
  0.994734,
  0.994403,
  0.994063,
  0.993716,
  0.993360,
  0.992997,
  0.992626,
  0.992247,			/*  40 C */
  0.991861,
  0.991467,
  0.991067,
  0.990659,
  0.990244,
  0.989822,
  0.989393,
  0.988957,
  0.988515,
  0.988066,			/*  50 C */
  0.987610,
  0.987148,
  0.986680,
  0.986205,
  0.985723,
  0.985236,
  0.984743,
  0.984243,
  0.983737,
  0.983226,			/*  60 C */
  0.982708,
  0.982185,
  0.981655,
  0.981120,
  0.980580,
  0.980034,
  0.979482,
  0.978924,
  0.978361,
  0.977793,			/*  70 C */
  0.977219,
  0.976640,
  0.976056,
  0.975466,
  0.974871,
  0.974271,
  0.973665,
  0.973055,
  0.972439,
  0.971819,			/*  80 C */
  0.971193,
  0.970562,
  0.969926,
  0.969286,
  0.968640,
  0.967990,
  0.967335,
  0.966674,
  0.966009,
  0.965340,			/*  90 C */
  0.964665,
  0.963986,
  0.963302,
  0.962613,
  0.961920,
  0.961222,
  0.960519,
  0.959812,
  0.959100,
  0.958384,			/* 100 C */
  0.957662,
  0.956937,
  0.956207,
  0.955472,
  0.954733,
  0.953989,
  0.953240,
  0.952488,
  0.951730,
  0.950968};

/*
 * Calculate the SG at 14° C (it seems; it's a strange value for a
 * Fahrenheit-based conversion), given the SG at temperature temp.
 */

float
  corrected_SG (float SG, float temp, char temptype)
{
  float tempF, correction;

  if ((temptype == 'C') || (temptype == 'c'))
    tempF = C_to_F_degrees (temp);
  else
    tempF = temp;

  correction = A0 + A1 * tempF + A2 * tempF * tempF + A3 * tempF * tempF * tempF;
  return SG + (correction * 0.001);
  }

/*
 * Calculate the SG at 15° C (59° F), given the SG at temperature
 * temp.
 */

float
  corrected_SG_at_15C (float SG, float temp, char temptype)
{
  float C, correction;

  if ((temptype == 'F') || (temptype == 'f'))
    C = F_to_C_degrees (temp);
  else
    C = temp;

  /*
   * Use a linear interpolation for 0-3.98 C
   */
  if (C < 3.98)
    correction = -0.000032692 * C - 0.000740644;
  else if (C < 50)
    correction = -0.0008031922 - 0.0000473773 * C + 0.000007231263 * C * C
      - 0.00000003078278 * C * C * C;
  else
    correction = -0.005431719 + 0.0001963596 * C + 0.000002661056 * C * C;

  return SG + correction;
  }

/* Reference: Manning, M.P., Understanding Specific Gravity and Extract,
 *            Brewing Techniques, 1,3:30-35 (1993)
 *
 * Equation is:  P = -676.67 + 1286.4SG - 800.47SG**2 + 190.74SG**3
 */

float
  SG_to_Plato (float SG)
{
  return -676.67 + 1286.4 * SG - 800.47 * SG * SG + 190.74 * SG * SG * SG;
  }

float Plato_to_SG (float Plato)
{
  return Plato / (258.6 - ((Plato / 258.2) * 227.1)) + 1.0;
  }

/* Reference: Manning, M.P., Understanding Specific Gravity and Extract,
 *            Brewing Techniques, 1,3:30-35 (1993)
 *
 * Equation is:  We = P*SG*R/100
 *       where:  P = degrees Plato @ 60F
 *              SG = Specific Gravity @ 60F
 *               R = density of water @ 60F = 8.338lb/USgal = 0.999kg/L
 */

float
  lbs_extract_per_gal_SG (float SG_raw, float temp, char temptype)
{
  float SG, Plato;

  SG = corrected_SG (SG_raw, temp, temptype);
  Plato = SG_to_Plato (SG);

  return SG * Plato * 8.338 * 0.01;
  }

/* Reference: Manning, M.P., Understanding Specific Gravity and Extract,
 *            Brewing Techniques, 1,3:30-35 (1993)
 *
 * Equation is:  SG2pts = SG1pts * (V1/V2)
 *       where:  SGpts  = SG - 1
 */

float
  SG_change_with_vol (float SG0, float vol0, float vol1)
{
  return ((SG0 - 1.0) * (vol0 / vol1)) + 1.0;
  }

/*
 * Equation is:  V1 = V0 * (SG0pts/SG1pts)
 *       where:  SGpts  = SG - 1
 */

float
  vol_change_with_SG (float vol0, float SG0, float SG1)
{
  return ((SG0 - 1.0) / (SG1 - 1.0)) * vol0;
  }


float
  temp_to_F_degrees (float temp, char temptype)
{
  if ((temptype == 'C') || (temptype == 'c'))
    return C_to_F_degrees (temp);
  else
    return temp;
  }


float
  temp_to_C_degrees (float temp, char temptype)
{
  if ((temptype == 'F') || (temptype == 'f'))
    return F_to_C_degrees (temp);
  else
    return temp;
  }


float
  C_to_F_degrees (float temp)
{
  return (temp * 1.8) + 32.0;
  }


float
  F_to_C_degrees (float temp)
{
  return (temp - 32.0) / 1.8;
  }


float
  lbs_to_kgs (float lbs)
{
  return lbs * 0.453592;
  }


float
  kgs_to_lbs (float kgs)
{
  return kgs * 2.204623;
  }


float
  gals_to_liters (float gals)
{
  return gals * 3.785412;
  }


float
  liters_to_gals (float liters)
{
  return liters * 0.264272;
  }


float
  qts_to_liters (float qts)
{
  return qts * 0.946353;
  }


float
  liters_to_qts (float liters)
{
  return liters * 1.056688;
  }

/* Calories due to alcohol per 12 ounces of beer */

float
  alc_calories_per_12_oz (float OSG, float FSG)
{
  return (1881.22 * FSG) * ((OSG - FSG) / (1.775 - OSG));
  }

/* Calories due to extract per 12 ounces of beer */

float
  ext_calories_per_12_oz (float OSG, float FSG)
{
  return 3550.0 * FSG * ((0.1808 * OSG) + (0.8192 * FSG) - 1.0004);
  }

/* Note temp MUST be fahrenheit
 */

float
  CO2_pressure (float vol, float tempF)
{
  float KP;						    /* keg pressure */

  KP = -16.6999 - (0.0101059 * tempF)
    + (0.00116512 * tempF * tempF)
      + (0.173354 * tempF * vol)
	+ (4.24267 * vol)
	  - (0.0684226 * vol * vol);
  return KP;
  }

/* Find volume of decoction to raise mash temp from initial temp to target temp. */

float decoct (float vt, float tt, float t0, float td)

/*       vt,   Total Volume */
/*       tt,   Target temp */
/*       t0,   Initial temp */
/*       td;   decoction temp */

{
  return (vt * (tt - t0)) / (td - t0);
  }

/* Find volume addition at given temp to raise mash temp to target.
 */

float
  beer_step_add_vol (float vm, float tm, float ta, float tf)

/*       vm,   Volume mash */
/*       tm,   Temperature mash */
/*       ta,   temperature addition */
/*       tf;   Temperature final */

{
  float va;						    /* Volume addition - what we are calculating */

  va = (vm * (tm - tf)) / (tf - ta);
  return va;
  }

/* Returns the Relative density difference between water at the reference
 * temp, t0, and the target temp, t1.
 */

float
  beer_relative_density_of_water (float t0, float t1)
{
  int tlow, thigh;

  float d0, d1, dlow, dhigh;

  extern float density_water[BEER_NUM_DENS];

  /*
   * First find the density of water at the reference temperature by interpolation
   */
  if (t0 == 0.0)
    d0 = density_water[0];
  else
    {
    tlow = floor ((double) t0);
    thigh = ceil ((double) t0);

    dlow = density_water[tlow];
    dhigh = density_water[thigh];

    d0 = dlow;
    d0 += (t0 - tlow) * (dhigh - dlow);
    }

  /*
   * Next find the density of water at the target temperature by interpolation
   */
  if (t1 == 0.0)
    d1 = density_water[0];
  else
    {
    tlow = floor ((double) t1);
    thigh = ceil ((double) t1);

    dlow = density_water[tlow];
    dhigh = density_water[thigh];

    d1 = dlow;
    d1 += (t1 - tlow) * (dhigh - dlow);
    }

  return d1 / d0;
  }

/* Returns the Relative density difference between water at the reference
 * temp, t0, and the target temp, t1.
 */

float
  beer_relative_volume_of_water (float t0, float v0, float t1)
{
  float rd1;						    /* relative density of target volume at target temp */

  rd1 = beer_relative_density_of_water (t0, t1);
  return v0 / rd1;
  }


/* The following are Glenn Tinseth's hopping functions */
float
  tinseth_hop_gravity_adjustment (float SGB)

/* SGB   Specific Gravity of the boil */

{
  return 1.65 * pow (0.000125, (double) SGB - 1.0);
  }


float
  tinseth_hop_boil_time_factor (float boil_time)
{
  return (1.0 - exp (-0.04 * (double) boil_time)) / 4.15;
  }


float
  tinseth_hop_utilization (float SGB, float boil_time)
{
  float time_factor, SG_factor;

  time_factor = tinseth_hop_boil_time_factor (boil_time);
  SG_factor = tinseth_hop_gravity_adjustment (SGB);
  return SG_factor * time_factor;
  }

float
  tinseth_hop_ibu_from_weight (float Util, float Acid,
			       float Wt, float Vol, float GA, float WtFactor)

/*    Util,      Percent Hop utilization as fraction ( Util < 1.0)
      Acid,      Percent Alpha acid content of hops as fraction
      Wt,        Weight of Hops in grams
      Vol,       Volume of batch size in Liters or Gallons
      GA,        Gravity adjustment
      WtFactor;  Weight conversion value for metric or english measures
      */
{
  return (Util * Acid * Wt * WtFactor) / Vol;
  }


float
  tinseth_hop_weight_from_ibu (float Util, float Acid,
			       float IBU, float Vol, float GA, float WtFactor)

/*       Util,     Percent Hop utilization as fraction ( Util < 1.0) */
/*       Acid,     Percent Alpha acid content of hops as fraction */
/*       IBU,      Desired IBU */
/*       Vol,      Volume of batch size in Liters */
/*       GA,       Gravity adjustment */
/*       WtFactor; Weight conversion value for metric or english measures */

{
  return (Vol * IBU) / (Util * Acid * WtFactor);
  }

/* Jackie Rager's hopping functions */


float
  rager_hop_gravity_adjustment (float SGB)

/* SGB;    Specific Gravity of the boil */
{
  if (SGB <= RAGER_HIGH_SPECIFIC_GAVITY)
    return 0.0;

  return (SGB - RAGER_HIGH_SPECIFIC_GAVITY) / (0.2);
  }


float
  rager_hop_util (int boil_minutes)
{
  if (boil_minutes == 0)
    return 0.0;
  else if (boil_minutes <= 5)
    return 0.05;
  else if (boil_minutes <= 10)
    return 0.06;
  else if (boil_minutes <= 15)
    return 0.08;
  else if (boil_minutes <= 20)
    return 0.101;
  else if (boil_minutes <= 25)
    return 0.121;
  else if (boil_minutes <= 30)
    return 0.153;
  else if (boil_minutes <= 35)
    return 0.188;
  else if (boil_minutes <= 40)
    return 0.228;
  else if (boil_minutes <= 45)
    return 0.269;
  else if (boil_minutes <= 50)
    return 0.281;
  else if (boil_minutes <= 55)
    return 0.291;
  else
    return 0.30;
  }


float
  rager_hop_ibu_from_weight (float Util, float Acid,
			     float Wt, float Vol, float GA, float WtFactor)

/*       Util,     Percent Hop utilization as fraction ( Util < 1.0) */
/*       Acid,     Percent Alpha acid content of hops as fraction */
/*       Wt,       Weight of Hops in grams */
/*       Vol,      Volume of batch size in Liters or Gallons */
/*       GA,       Gravity adjustment */
/*       WtFactor; Weight conversion value for metric or english measures */

{
  return (Util * Acid * Wt * WtFactor) / (Vol * (1.0 + GA));
  }


float
  rager_hop_weight_from_ibu (float Util, float Acid,
			     float IBU, float Vol, float GA, float WtFactor)

/*       Util,     Percent Hop utilization as fraction ( Util < 1.0) */
/*       Acid,     Percent Alpha acid content of hops as fraction */
/*       IBU,      Desired IBU */
/*       Vol,      Volume of batch size in Liters */
/*       GA,       Gravity adjustment */
/*       WtFactor; Weight conversion value for metric or english measures */

{
  return (Vol * IBU * (1.0 + GA)) / (Util * Acid * WtFactor);
  }

/* .39661 + 0.0017091*Po + (1.0788E-5)*Po^2  AJ DeLange
 */

float
  ballings_factor (float Plato)
{
  return 0.39661 + (0.0017091 * Plato) + (1.0788e-5 * Plato * Plato);
  }


float
  alcohol_by_weight_fix (float OSG, float FSG)
{
  float Alc, OE, AE, RE;
  
  OE = SG_to_Plato (OSG);
  AE = SG_to_Plato (FSG);
  
  RE = (0.1808 * OE) + (0.8192 * AE);
  Alc = (OE - RE) / (2.0665 - (0.010665 * OE));
  
  return Alc;
  }


float
  alcohol_by_weight_AJ (float OSG, float FSG)
{
  float PO,						    /* Plato original (corrected) */
  PF,							    /* Plato finished (corrected) */
  BF,							    /* Balling factor */
  ABW;
  
  PO = SG_to_Plato (OSG);
  PF = SG_to_Plato (FSG);
  BF = ballings_factor (PO);
  ABW = (PO - PF) * BF;
  return ABW;
  }

/* Alcohol by weight to alcohol by volume */

float
  abw_to_abv (float abw)
{
  return abw / 0.793573;
  }

/*
 * RI = 1.33302 + 0.001427193(B) + 0.000005791157(B^2)
 *
 * where:
 *
 * B  = measured refractivity in Brix
 * RI = calculated Refractive Index
 */

float
  brix_to_ri (float brix)
{
  return (float) 1.33302 + 0.001427193 * brix + 0.000005791 * brix * brix;
  }

/*
 * SG = 1.001843 - 0.002318474(OG) - 0.000007775(OG^2) - 0.000000034(OG^3) +
 *      0.00574(AG) + 0.00003344(AG^2) + 0.000000086(AG^3)
 *
 * where:
 *
 * SG = estimated specific gravity of the sample
 * OG = Original Gravity of the batch (in Brix)
 * AG = Apparent Gravity of the sample (in Brix)
 */

float
  brix_to_fg (float ob, float ab)
{
  double OB, AB, SG;

  OB = (double) ob;
  AB = (double) ab;

  SG = 1.001843 - (0.002318474 * OB) - (0.000007775 * OB * OB)
    - (0.000000034 * OB * OB * OB) + (0.00574 * AB)
      + (0.00003344 * AB * AB) + (0.000000086 * AB * AB * AB);

  return (float) SG;
  }

/* Alcohol by Weight from Specific Gravity and Degrees Brix of finished beer */

float
  brix_to_abw (float sg, float ab)
{
  double RI;

  RI = brix_to_ri (ab);

  return (float) (1017.5596 - (277.4 * sg) + RI * ((937.8135 * RI) - 1805.1228));
  }

#if 0
/* Alcohol by Weight from Specific Gravity and Degrees Brix of finished beer */
/* This function is redundant: it's an inline version of brix_to_abw() with the
 * parameters reversed. */

float
  brix_sg_to_abw (float ab, float sg)
{
  double SG, AB, RI;

  SG = sg;
  AB = ab;

  RI = 1.33302 + 0.001427193 * AB + 0.000005791 * AB * AB;

  return (float) (1017.5596 - (277.4 * SG) + RI * ((937.8135 * RI) - 1805.1228));
  }
#endif

/* Specific Gravity @ 15C from Degrees Brix @ 20C */

float
  brix_to_sg (float brix)
{
  return 1.000898
    + 0.003859118 * brix
    + 0.00001370735 * brix * brix
    + 0.00000003742517 * brix * brix * brix;
  }