Thread: Wacky numbers in simple 'for' loop

You need to compile with warnings:
main.c:18:45: warning: format ‘%d’ expects argument of type ‘int’, but argument 3 has type ‘double’ [-Wformat=]
18 | printf("\nIndex = %d\t\tboil point = %d", index, bp[index]);

thmm-- Thanks. I did see it eventually. I'm a total noob but still...

Salem-- Appreciate it. I'd gotten that advice here before. Different context. Have to keep it in mind.

Taking to heart the aforementioned cumulative addition errors; I'd like to improve the way my functions iterate by adding lots of small numbers. Functions calc_bp, calc_xw and calc_yd are the same equation solved in different terms.
I guess the question is how to converge on the problem condition (psat sums = atm_p) in the shortest amount of steps.
Advice appreciated. I haven't seen any coding examples of this type before.

Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


#define W_DENS 0.998       //Water density at 20 degrees C
#define E_DENS 0.789        //Ethanol density at 20 degrees C
#define MM_W 18.0153       //Molar mass of water - grams per mol. Source NIST
#define MM_E 46.0684       //Molar mass of ethanot grams per mol. Source NIST
#define MAX_NUM 81


#define GAL2L 3.785411784  //Gallons to liters factor




/*********FUNCTION PROTOYPES*************************/
double get_input(void);
double abv_to_abw(double abv);  /**convert alcohol by volume to alcohol by weight*/
double abw_to_mf_e(double abw);  /**convert alcohol by weight to mol fraction*/


double mf_to_abw(double mf_e);   /**convert mol fraction to alcohol by weight*/
double abw_to_abv(double abw);  /**convert alcohol by weight to alcohol by volume*/


double calc_total_mol(double wash_l, double abv);


/**calc mol fraction converts from abv to abw then from abw to ethanol molar fraction*/
double calc_mol_fract(double abv, double (*abv_to_abw)(double abv), double (*abw_to_mf_e)(double abw));


/**van Laar activity coefficient for ethanol and water to model for non-ideal mix*/
double calc_gamma_e(double x1, double x2);
double calc_gamma_w(double x1, double x2);


/**Find saturated vapor pressure*/
double calc_psat(double A, double B, double C, double bp, double gamma, double x);


double calc_boil_pt(double abv, double atm_p, double mol_fract,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));


/**calc xw finds boiler mol fraction given boil point and Antoines Eq and Raoults law.*/
double calc_xw(double atm_p, double mol_fract, double boil_pt,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));


/**calc yd finds distllate purity given boil point and boiler mol fraction */
double calc_yd(double atm_p, double mol_fract, double boil_pt,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x));


void write_vle(double *boil_pt, double *xw, double *yd, double *W, double *D, double *dist_avg);


//Antoine Coefficients (C)
double eA = 8.13484;
double eB = 1662.48;
double eC = 238.131;


double wA = 8.05573;
double wB = 1723.64;
double wC = 233.076;


int main()
{
    //printf("\nEnter Wash Gallons: ");
    double wash_l = 100; //get_input() * GAL2L;


    //printf("\nEnter ABV (between 0 and 1): ");
    double abv = 0.40; //get_input();


    //printf("\nEnter Atmospheric pressure (mmHg): ");
    double atm_p = 760; //get_input();




    double boil_pt[MAX_NUM] = {0};
    double xw[MAX_NUM] = {0};/**Boiler mol fraction (Purity)*/
    double W[MAX_NUM] = {0};/**boiler molar volume*/


    double D[MAX_NUM] = {0};/**Distillate molar volume*/
    double yd[MAX_NUM] = {0};/**Molar fraction of ethanol vapor (Purity)*/


    double recip[MAX_NUM] = {0};
    double integral[MAX_NUM] = {0};
    double dist_avg[MAX_NUM] = {0};


    /**Determine initial conditions for VLE calculations*/
    xw[0] = calc_mol_fract(abv, abv_to_abw, abw_to_mf_e);      /**Purity to start*/
    boil_pt[0] = calc_boil_pt(abv, atm_p, xw[0], calc_gamma_e, calc_gamma_w, calc_psat); /**Boil point to start*/
    W[0] = calc_total_mol(wash_l, abv);                        /**Molar volume to start*/


    printf("\nInitial boil point = %.3lf", boil_pt[0]);
    printf("\nInitial wash purity = %.3lf", xw[0]);
    printf("\nInitial wash molar vol = %.3lf", W[0]);


    yd[0] = calc_yd(atm_p, xw[0], boil_pt[0], calc_gamma_e, calc_gamma_w, calc_psat);
    /**initial values of integral[0] and D[0] distillate volume in mols are zero*/


    /**for loop is Rayliegh Eq.*/
    for ( int index = 1; index < MAX_NUM; index++)
    {
      boil_pt[index] = index * 0.2 + boil_pt[0];


      /**xw[0] = Initial boiler wash purity*/
      xw[index] = calc_xw(atm_p, xw[0], boil_pt[index], calc_gamma_e, calc_gamma_w, calc_psat);/**Boiler liquid purity at instant*/


      /**distillate purity = vapor %ethanol in mols- condensed vapor = distillate*/
      yd[index] = calc_yd(atm_p, xw[0], boil_pt[index], calc_gamma_e, calc_gamma_w, calc_psat); /**distillate purity at instant*/


      recip[index] = 1 / (yd[index] - xw[index]);


      integral[index] = integral[index-1]+ (xw[index]- xw[index-1])* ((recip[index] + recip[index-1])/2);


      W[index] = W[0] * pow(10, integral[index]); /**Boiler liquid mols at instant*/


      D[index] = W[0] - W[index];/**Distillate molar volume*/


      dist_avg[index] = (W[0] *xw[0] -xw[index] *yd[index])/D[index];


      write_vle(boil_pt, xw, yd, W, D, dist_avg);
    }
    system("PAUSE");
    return 0;
}
/*******************************/




/*********FUNCTION DEFINITIONS*************************/
double get_input(void)
{
double input = 0;
scanf("%lf", &input);
return input;
}
/**************************************************/
//alcohol conversion functions abv to abw, and abw to abv results are between 0 and 1 and only valid for input between 0 and 1
//From "On the Conversion of Alcohol by Edwin Croissant -  R1: 2014*02-16
double abw_to_abv(double abw)
{
    const double a = -0.000039705486746795932;
    const double b =  1.2709666849144778;
    const double c = -0.40926819348115739;
    const double d =  2.0463351302912738;
    const double f = -7.8964816507513707;
    const double g =  15.009692673927390;
    const double h = -15.765836469736477;
    const double i = 8.8142267038252680;
    const double j = -2.0695760421183493;


    double temp = j;
    temp = temp * abw +i;
    temp = temp * abw +h;
    temp = temp * abw +g;
    temp = temp * abw +f;
    temp = temp * abw +d;
    temp = temp * abw +c;
    temp = temp * abw +b;
    double abv = temp * abw +a;


    return abv;
}
/*******************************/


double abv_to_abw(double abv)
{
    const double a = 0.00018684999875047631;
    const double b = 0.77602465132552556;
    const double c = 0.41803095099103116;
    const double d = -2.5221614925275091;
    const double f = 9.5827123045656251;
    const double g = -19.928886159385002;
    const double h = 24.165120890385651;
    const double i = -15.830262207383321;
    const double j = 4.3390473620304988;


    double temp = j;
    temp = temp * abv +i;
    temp = temp * abv +h;
    temp = temp * abv +g;
    temp = temp * abv +f;
    temp = temp * abv +d;
    temp = temp * abv +c;
    temp = temp * abv +b;
    double abw = temp * abv +a;


    return abw;
}
/*******************************/


double abw_to_mf_e(double abw)
{
   double mf_e = abw / (abw + MM_E / MM_W * (1-abw));
   return mf_e;
}
/*******************************/


double calc_total_mol(double wash_l, double abv)
{
    double water_mol = wash_l *(1-abv) *1000 *W_DENS /MM_W;
    double eth_mol = wash_l *1000 *abv *E_DENS /MM_E;


    return water_mol + eth_mol;
}
/***********************************/
double calc_mol_fract(double abv, double (*abv_to_abw)(double abv), double (*abw_to_mf_e)(double abw))
{
    double abw = abv_to_abw(abv);
    double eth_mol_fract = abw_to_mf_e(abw);
    return eth_mol_fract;
}


/***********************************/
   double calc_boil_pt(double abv, double atm_p, double mol_fract,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
   double bp = 70;                //Seed temperature at loop beginning
   double x1 = mol_fract;         //ethanol component of mixture in mols
   double x2 = 1-x1;              //water component of mixture in mols
   double psat1 = 0;
   double psat2 = 0;
   double gamma_e = 0;
   double gamma_w = 0;


   gamma_e = calc_gamma_e(x1,x2);
   gamma_w = calc_gamma_w(x1,x2);
    do
    {
        psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
        psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);
        bp = bp + 0.0001;
    }
    while (psat1 + psat2 < atm_p);
    return bp;
}
/************************************************/
/*Calculate Van Laar activity coefficient for ethanol*/
double calc_gamma_e(double x1, double x2)
{
double const a12 = 1.68811;/*a12 and a21 are Van Laar activity model constants*/
double const a21 = 0.95268;
double result = exp((a12 * pow(x2, 2)) / (pow(((a12 * x1 / a21) + x2), 2)));
return (result);
}
/***************************************************/
/*Calculate Van Laar activity coefficient for water*/
double calc_gamma_w(double x1, double x2)
{
double const a12 = 1.68811;
double const a21 = 0.95268;
double result = exp((a21 * pow(x1, 2)) /  (pow(((a21 * x2 / a12) + x1), 2)));
return (result);
}
/***************************************************/
/* Antoine Eq. to calculate saturated vapor pressure for water and ethanol */
double calc_psat(double A, double B, double C, double bp, double gamma, double x)
{
double psat = x*gamma*pow(10, (A-((B)/(C + bp))));
return psat;
}
/***************************************************/


/***************************************************/
double calc_xw(double atm_p, double mol_fract, double bp,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
   double x1 = mol_fract;         //ethanol component of mixture in mols
   double x2 = 1-x1;              //water component of mixture in mols
   double psat1 = 0;
   double psat2 = 0;
   double gamma_e = 0;
   double gamma_w = 0;


   do
    {
       gamma_e = calc_gamma_e(x1,x2);
       gamma_w = calc_gamma_w(x1,x2);


       psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
       psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);


       x1 -= 0.000001;
       x2 = 1-x1;
   }
    while (psat1 + psat2 > atm_p);
    return x1;
}


/**********************************************/
double calc_yd(double atm_p, double mol_fract, double bp,
         double (*calc_gamma_e)(double x1, double x2),
         double (*calc_gamma_w)(double x1, double x2),
         double (*calc_psat)(double A, double B, double C, double bp, double gamma, double x))
{
   double x1 = mol_fract;         //ethanol component of mixture in mols
   double x2 = 1-x1;              //water component of mixture in mols
   double psat1 = 0;
   double psat2 = 0;
   double gamma_e = 0;
   double gamma_w = 0;


   do
    {
       gamma_e = calc_gamma_e(x1,x2);
       gamma_w = calc_gamma_w(x1,x2);


       psat1 = calc_psat(eA, eB, eC, bp, gamma_e, x1);
       psat2 = calc_psat(wA, wB, wC, bp, gamma_w, x2);
       x1 -= 0.000001;
       x2 = 1-x1;
   }
    while (psat1 + psat2 > atm_p);
    double y_d = psat1 / atm_p;
    return y_d;
}


/*****************************************************************/
void write_vle(double *boil_pt, double *xw, double *yd, double *W, double *D, double *dist_avg)
{
   //Write to screen
   for(int i = 0; i < MAX_NUM; i++)
   {
     printf("\nindex\tboil_pt\t\tboiler mol\tdist purity\twash mols\tdist mols\tdist avg\n");
     printf("%d\t%.2lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\t\t%.3lf\n",
            i, boil_pt[i], xw[i], yd[i], W[i], D[i], dist_avg[i]);
   }


 //Write to file
 /*FILE *pWrite;
 pWrite = fopen("vle.dat", "w");
    if (pWrite == NULL)
    printf("\nFile not opened\n");
    else
    fprintf(pWrite, "\nboil_pt\t\txw\t\tyd\tW\tD\tdist_avg\n");


    for(int i = 0; i < MAX_NUM; i++)
      {
         fprintf(pWrite,"%.2lf\t\t%.3lf\t%.3lf\t%.3lf\t\t%.3lf\t%.3lf\t%.3lf\t%.3lf\n",
            boil_pt[i], xw[i], yd[i], W[i], D[i], D[i]);
      }
      */
}