Thread: Wacky numbers in simple 'for' loop

  1. #1
    Registered User
    Join Date
    Sep 2020
    Posts
    31

    Wacky numbers in simple 'for' loop

    It seems like I'm not initializing the array to zero. Trying to get a list incremented from 1 by 0.2. Simple. What am I missing?

    Thanks.

    Edit: Never mind


    Code:
    #include <stdio.h>#include <stdlib.h>
    
    #define MAX 80
    
    int main()
    {
       double bp[MAX] = {0};
       bp[0] = 1;
    
    
       for ( int index = 1; index < MAX; index++)
        {
          bp[index] = bp[index-1] +0.2;
          printf("\nIndex = %d\t\tboil point = %d", index, bp[index]);
        }
    
    
        return 0;
    }
    Last edited by Buckeye Bing; 12-26-2021 at 10:38 AM. Reason: solved

  2. #2
    Registered User
    Join Date
    Sep 2020
    Posts
    150
    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]);

  3. #3
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    39,435
    Floating point numbers are inexact representations (aka approximations) of most real numbers.
    What Every Computer Scientist Should Know About Floating-Point Arithmetic
    Floating Point

    So repeated addition is going to accumulate a lot of errors.

    It might be better to do
    bp[index] = index * 0.2 + 1;
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  4. #4
    Registered User
    Join Date
    Sep 2020
    Posts
    31
    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.

  5. #5
    Registered User
    Join Date
    Sep 2020
    Posts
    31

    Repeated addition in functions

    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]);
          }
          */
    }

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. How to for() loop with hex numbers?
    By monk64 in forum C Programming
    Replies: 5
    Last Post: 04-08-2010, 11:25 PM
  2. Wacky Partition Layout
    By Yarin in forum A Brief History of Cprogramming.com
    Replies: 1
    Last Post: 01-19-2009, 04:14 AM
  3. while loop and numbers
    By rahulsk1947 in forum C Programming
    Replies: 5
    Last Post: 02-13-2006, 05:59 AM

Tags for this Thread