help with looping

[sai.aleenus]

if f(x,y) be a function (x^3 * y^-4) + (x^-2 *y^3) with 2<x<10 and 1.1<y<3.5.f(0,0) = f(9,1.2). how do you develop a c program of loop(x decremented by 1 unit) with in which another loop (y incremented by 0.8(K) initially and then as the y reaches a value where incrementation by 0.8 makes y>3.4 make 0.8 to 0.8/2 and follow this till K>0.002) and the values generated are continuous.i.e the output is in the form

iteration       x      y         f(x,y)evaluated.
          0     9      1.2
          1     9       2.0
          2     9       2.8       (note here 2.8 + 0.8 > 3.4 so K = 0.8/2)
          3     9       3.2        (note 3.2+0.4 > 3.4 so  K = 0.4/2)
          4     9       3.4
          ..................
          .................
          n     8      1.2
          .................
          ................
          n+j   8      3.4
          ................
          ...............
          .,..............
          n+m   3     1.2
          .................
          ....................
     the f is to be analysed and the min of the calculated values is to be given as output.
      can anyone please help me with the code for this..
      i have developed a code but its getting struck....

[MonsieurPointer]

Please show us what you have so far as code is concerned. Where are you getting stuck?

[sai.aleenus]

hey...i am posting the code for which i was actually facing the problem..
here the f(x,y) is the function bhpeqns()...x is fill_s[] and y is fill_h with ranges 12<=x<=19 and 0.6096<=y<=1.8288
please help me with the following code...

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<stdarg.h>
#include<graphics.h>
#include<iostream.h>

#define to 273.15
#define t9 273.16
#define spvolro 8.31432
#define xvmin 300
#define xvmax 750
#define ypmin 0
#define ypmax 1
#define pd_fill_n 0.8713
#define pd_fill_a (1.355 * pow(10,-4))
#define pd_fill_b (-3.718 * pow(10,-5))
#define pd_fill_m 0.3485
#define pd_fill_c 0.00589
#define pd_fill_o 0.874
#define pd_fill_d 0.119
#define pd_fill_p 0.893
#define pd_fill_q 0.461

float sat_vap(float twbi);
float sp_hum(float svp,float tdbi,float twbi);
float sat_vape(float twbo);
float sp_enthalpy(float tdbi,float sphum);
float sp_vole(float tdbi,float twbi,float svp);
float row(float sphum,float spvol);

float bhpeqns(float fh,float fs);
void kkt_optimize();

float thw=43.5,tcw=33.0f,twbi=28.82,tdbi=33.0f,app,r,fand=9.144,fann=0.865;
float aia,aih=5.5,edisa,edisa_f,cell=16.63,celw=16.63,cela,twbo,tdbo,gpm=24000,gpmpc,wl,spe_in,wl_f;
int ais = 3,celn = 7;
float kavlc[] = { 0.04503, 0.5467, 0.2516,0.637,0.429};// c,p,q,r,s
float svpp[] = { -2948.997118,  -2.1836674, -0.000150474, -0.0303738468,0.00042873,  4.76955, 25.83220018,24.72528008};
float sphup[] = {18.01534, 28.9645,  101325.0f,  0.000666,52.2542}; //m1,m2,po,c
float spep[] = {  1.00568, 1.84598, 2500.84};       //c2,c3,lo
float  spec2 = 1.00568, spec3 = 1.84598, spelo = 2500.84;float row_fill_f;

float fh,fs;





 int main()
{

     kkt_optimize();




    getch();
    return 0;
  //  getch();
 //   return 0;

}


float sat_vap(float twbi)
{
    float too,sat[7],fsvp;

    too = twbi + to;

    sat[0] = svpp[0]/too;
    sat[1] = svpp[1] * log(too);
    sat[2] = abs(svpp[3] * (twbi-0.01));
    sat[3] = svpp[2] * (pow(0.1,sat[2]));//
    sat[4] = svpp[5] * (1-(t9/too));
    sat[5] = svpp[4] * (pow(10,sat[4]));//
    sat[6] = sat[0] + sat[1] +sat[3] + sat[5] + svpp[6] ;
    sat[7] = pow(10,sat[6]);                //sat vap pressure calculated
    fsvp = sat[7];

return (fsvp);
}

float sp_hum(float svp,float tdbi,float twbi)
{

float hum[5],fsphum;

 hum[0] = sphup[0]/sphup[1];//
 hum[1] = sphup[2] * sphup[3] * (tdbi-twbi);
 hum[2] = svp - hum[1];
 hum[3] = (hum[0] * hum[2])/(sphup[2] - hum[2]); // sp humidity is calculated
 fsphum = hum[3];
 fsphum = fabs(fsphum);

return (fsphum);
}

float sat_vape(float twbo)
{
    float too,sat[7],fsvpe;

    too = twbo + to;

    sat[0] = svpp[0]/too;
    sat[1] = svpp[1] * log(too);
    sat[2] = fabs(svpp[3] * (twbi-0.01));
    sat[3] = svpp[2] * (pow(0.1,sat[2]));//
    sat[4] = svpp[5] * (1-(t9/too));
    sat[5] = svpp[4] * (pow(10,sat[4]));//
    sat[6] = sat[0] + sat[1] +sat[3] + sat[5] + svpp[7]  ;
    sat[7] = pow(10,sat[6]);                //sat vap pressure calculated
    fsvpe = sat[7];

return (fsvpe);
}

float sp_enthalpy(float tdbi,float sphum)
{
    float spe[5],fspe;
    spe[0] = spep[0] * tdbi ;
    spe[1] = (spep[1] * tdbi) + spep[2];
    spe[3] = sphum * spe[1] ;
    spe[4] = spe[0] + spe[3];
    fspe = spe[4];
    fspe = fabs(fspe);
return(fspe);
}

float sp_vole(float tdbi,float twbi,float svp)
{
    float vol[5],fspvol;
    vol[0] = spvolro * (tdbi + to);
    vol[1] = sphup[2] * sphup[3] * (tdbi - twbi);
    vol[2] = sphup[2] - svp + vol[1];
    vol[3] = sphup[4] * vol[2] ;
    vol[4] = (vol[0] * 1000)/vol[3] ;
    fspvol = vol[4];
    fspvol = fabs(fspvol);
return(fspvol);
}
float row(float sphum,float spvol)
{
    float row;
    row = (1 + sphum)/spvol;
    row = fabs(row);
return(row);
}


float bhpeqns(float fh, float fs)
{


    float fh_f,fs_f,lgf,fs_ff;
    float svp_in,svph,svpi,svpj,svpk, svpl, svpm, svpn,too;//sat vap pressure
    float sphu_in,sphua,sphub,sphud,hin_f;//sp humidity
    float too1,too2,too3,too4,twbi1,twbi2,twbi3,twbi4,svp1,svp2,svp3,svp4,t1,t2,t3,t4,tdbi1,tdbi2,tdbi3,tdbi4;
    float sphu1,sphu2,sphu3,sphu4,spe1,spe2,spe3,spe4,spe;
    float h1,h2,h3,h4,h11,h22,h33,h44;
    float speat1,speat2,speat3,speat4,ra,h1a,h1b,h2a,h2b,h3a,h3b,h4a,h4b;
    float kav1,kav2,lg;
    float kav2a,kav2b,kav2c,kav2d,edisa_f;
    float fvel,fvel_f,l_f,l_s,g_f,g_s,spvol_in_s,spvol_in_f,row_in_s,row_in_f;
    float too_o,twbo,tdbo,svp_out,sphu_out,spvol_out_s,spvol_out_f,hout_s,hout_f,tdbo_f,row_out_s,row_out_f;
    float vel_in_s,vel_in_f,vel_out_s,vel_out_f,vel_fill_s,vel_fill_f,row_fill_s,row_fill_f,spvol_fill_s,spvol_fill_f;
    float wl_f,pdf1,pdf2,pdf3,pdf4,pdf5,pdf6,pd_il_inh20,pd_il_pa,cfm_fan_s,cfm_fan_f,pd_de_inh20,pd_de_pa,pd_fill_inh20,pd_fill_pa;
    float tpd_s_inh20,tpd_s_pa,pd_v_f,pd_v_s,bhp,bhp_w,tpd_f,tpd_s;
    float sphure,sphuref,spvolref,spvolre,predp;
    float pdvf2,pdvf3,bhhp1,bhhp2;



           r = thw - tcw ;
            app = tcw - twbi ;
            cela = cell*celw ;
            aia = ais*celw*aih;
            edisa = (0.785*(fand*fand));
            edisa_f = edisa*3.2808*3.2808;
            printf("\n edisa %f",edisa_f);
            gpmpc = gpm/celn ;
            wl = gpmpc/cela;
            wl_f = (gpmpc*4.402867)/(cela*3.2808*3.2808);// fpm


            twbi = twbi;
            svp_in = sat_vap(twbi);                //sat vap pressure calculated
            printf("\n sat vap p inlet is %f Pa",svp_in);

            tdbi = tdbi;
            twbi = twbi;
            sphu_in = sp_hum(svp_in, tdbi, twbi); // sp humidity is calculated
            printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);

            spe_in = sp_enthalpy(tdbi,sphu_in);//sp enthalpy is calculated
            printf("\n inlet enthalpy is(kj/kg) %f",spe_in);
            spe = spe_in;

            hin_f = ((spe_in*1.8)/4.1858)+7.686 ;
            printf("\n inlet enthalpy is(btu/lb) %f",hin_f);

            lg = 1.0f;
            do
            {

 //         t1 = tcw + (0.1*r);
 //         t2 = tcw + (0.4*r);
 //         t3 = thw - (0.4*r);
 //         t4 = thw - (0.1*r);

            twbi1 = twbi + (0.1*r);
            twbi2 = twbi + (0.4*r);
            twbi3 = twbi + (0.6*r);
            twbi4 = twbi + (0.9*r);

            tdbi1 = twbi + (0.1*r);
            tdbi2 = twbi + (0.4*r);
            tdbi3 = twbi + (0.6*r);
            tdbi4 = twbi + (0.9*r);

            svp1 = sat_vap(twbi1);                //sat vap pressure calculated
         //   printf("\n sat vap p1 is %f",svp1);


            sphu1 = sp_hum(svp1, tdbi1, twbi1);  // sp humidity is calculated
          //  printf("\n sp humidity1 is %f",sphu1);


            spe1 = sp_enthalpy(tdbi1,sphu1);//sp enthalpy is calculated
          //  printf("\n enthalpy1 is %f",spe1);

            svp2 = sat_vap(twbi2);                //sat vap pressure calculated
         //   printf("\n sat vap p2 is %f",svp2);


            sphu2 = sp_hum(svp2, tdbi2, twbi2);  // sp humidity is calculated
          //  printf("\n sp humidity3 is %f",sphu2);


            spe2 = sp_enthalpy(tdbi2,sphu2);//sp enthalpy is calculated
          //  printf("\n enthalpy2 is %f",spe2);



            svp3 = sat_vap(twbi3);                //sat vap pressure calculated
         //   printf("\n sat vap p3 is %f",svp3);


            sphu3 = sp_hum(svp3, tdbi3, twbi3);  // sp humidity is calculated
          //  printf("\n sp humidity3 is %f",sphu3);


            spe3 = sp_enthalpy(tdbi3,sphu3);//sp enthalpy is calculated
          //  printf("\n enthalpy3 is %f",spe1);


            svp4 = sat_vap(twbi4);                //sat vap pressure calculated
         //   printf("\n sat vap p4 is %f",svp4);


            sphu4 = sp_hum(svp4, tdbi4, twbi4);  // sp humidity is calculated
          //  printf("\n sp humidity4 is %f",sphu4);


            spe4 = sp_enthalpy(tdbi4,sphu4);//sp enthalpy is calculated
          //  printf("\n enthalpy4 is %f",spe4);

            speat1 = spe + (0.1*lg*r);
            speat2 = spe + (0.4*lg*r);
            speat3 = spe + (0.6*lg*r);
            speat4 = spe + (0.9*lg*r);

            h1a = spe1-speat1;
            h1b = ((h1a*1.8)/4.1858);
            h1 = h1b + 7.686 ;
            h2a = spe2-speat2;
            h2b = ((h2a*1.8)/4.1858);
            h2 = h2b + 7.686 ;
            h3a = spe3-speat3;
            h3b = ((h3a*1.8)/4.1858);
            h3 = h3b + 7.686 ;
            h4a = spe4-speat4;
            h4b = ((h4a*1.8)/4.1858);
            h4 = h4b + 7.686 ;//btu/lb
            h11 = 1/h1;
            h22 = 1/h2;
            h33 = 1/h3;
            h44 = 1/h4;

            ra = (r*1.8)+32;//in deg f

            kav1 = ra*0.25*(h11+h22+h33+h44);//in fps ---- demand

 //    fvel = (wl*8.33)/(lg*0.07);//in m/s
            fvel_f = (wl*8.33)/(lg*0.07*0.227125);//in fpm conversion factor 0.227125 used

            fh_f = fh/0.3048;
            fs_f = fs*0.00328084;//in ft
            lgf = 1/lg;
            fs_ff = 1/fs_f;

            kav2a = kavlc[0] * pow(lgf,kavlc[1]);
            kav2b = pow(fvel_f,kavlc[2]);
            kav2c = pow(fh_f,kavlc[3]);
            kav2d = pow (fs_ff,kavlc[4]);
            kav2 = kav2a*kav2b*kav2c*kav2d; //supply ntu

            lg = lg+0.001;
            continue;


            }while(kav1 <=kav2 && lg<=3.0f && lg>=1.0f);


               printf("\n the l/g ratio required for the performance is %f ",lg);
               kav2 = kav2a*kav2b*kav2c*kav2d;
               printf("\n the KAV/L_2  for the performance is %f ",kav2);
               kav1 = ra*0.25*(h11+h22+h33+h44);
               printf("\n the KAV/L_1  for the performance is %f ",kav1);

 // fvel = (wl*8.33)/(lg*0.07);//in m/s
               fvel_f = (wl*8.33)/(lg*0.07*0.00508);
               //printf("\n the velocity at the fill is %f m/s or %f fpm \n",fvel,fvel_f);

               l_f = (gpmpc * 4.402 * 8.33);//lb/min
               l_s = l_f * 0.0075598;// kg/sec
               printf("\n the mass flow rate of water is %f kg/s or %f lb/min",l_s,l_f);

              // lg = lg + 0.5;
               g_f = l_f/lg;
               g_s = l_s/lg;//g_f = g_s/(0.0075598);
               printf("\n the mass flow rate of air is %f kg/s or %f lb/min",g_s,g_f);

              printf("\n sat vap p inlet is %f Pa",svp_in);
               printf("\n sp humidity inlet is %f kg/kg of dry air",sphu_in);
               spvol_in_s = sp_vole(tdbi,twbi, svp_in);;//cum/kg
               spvol_in_f = spvol_in_s / 0.0624279606;//cuft/lb
               printf("\n the inlet specific volume is %f cum/kg or %f cu ft/lb ",spvol_in_s,spvol_in_f);//

               row_in_s = row(sphu_in,spvol_in_s);;
               row_in_f = row_in_s/16.0184634;
               printf("\n the inlet density is %f kg/cum or %f lb/cu ft",row_in_s,row_in_f);

               hout_s = spe + (lg*r);
               hout_f = ((hout_s*1.8)/4.1858)+7.686 ;
               printf("\n the inlet enthalpy is %f kj/kg or %f btu/lb",spe_in,hin_f);
               printf("\n the exit enthalpy is %f kj/kg or %f btu/lb",hout_s,hout_f);

               tdbo = hout_s/1.0062;
               printf("\n tdbo is %f",tdbo);
           //    tdbo_f = hout_f/0.240;//deg f
               twbo = tdbo;


            svp_out = sat_vape(twbo);               //sat vap pressure calculated
            printf("\n sat vap p exit is %f Pa  ",svp_out);



            sphu_out = sp_hum(svp_out,tdbo,twbo);

            printf("\n sp humidity exit is %f kg/kg of satureated air",sphu_out);


                spvol_out_s = sp_vole(tdbo,twbo,svp_out);//cum/kg
               spvol_out_f = spvol_out_s / 0.0624279606;//cuft/lb
               printf("\n the exit specific volume is %f cum/kg or %f cu ft/lb ",spvol_out_s,spvol_out_f);

               row_out_s = row(sphu_out,spvol_out_s);

               row_out_f = row_out_s/16.0184634;
               printf("\n the exit density is %f kg/cum or %f lb/cu ft",row_out_s,row_out_f);

               vel_in_s = (g_s * spvol_in_s)/aia ;//m/s
               vel_in_f = vel_in_s/0.00508;//fpm
               vel_out_s = (g_s*spvol_out_s)/cela ;
               vel_out_f =  vel_out_s/0.00508;

               row_fill_s = (row_in_s + row_out_s)/2 ;
               row_fill_f = (row_in_f + row_out_f)/2 ;
               spvol_fill_s = (spvol_in_s+spvol_out_s)/2;
               spvol_fill_f = (spvol_in_f+spvol_out_f)/2;   //

               vel_fill_s = (g_s*spvol_fill_s)/cela ;
               vel_fill_f = vel_fill_s/0.00508;

               printf("\n the specific volume at fill is %f cum/kg or %f cu ft/lb ",spvol_fill_s,spvol_fill_f);
               printf("\n the density at fill is %f kg/cum or %f lb/cu ft",row_fill_s,row_fill_f);
               printf("\n the velocity of air at inlet is %f m/s or %f fpm ",vel_in_s,vel_in_f);
               printf("\n the velocity of air at fill is %f m/s or %f fpm ",vel_fill_s,vel_fill_f);
               printf("\n the velocity of air at exit is %f m/s or %f fpm ",vel_out_s,vel_out_f);

               cfm_fan_s = g_s * spvol_out_s;//cum/s
               cfm_fan_f = g_f * spvol_out_f ;//cfm_fan_s * pow(3.2808,3) * 60;or cfm_fan_s * 2118.89;
               printf("\n the volume flow rate of air at fan CFAN is %f cum/s or %f cu ft/min \n",cfm_fan_s,cfm_fan_f);



               pd_il_inh20 = (3 * ((2*pow(10,-7)*vel_in_f*vel_in_f)+(3*pow(10,-6)*vel_in_f))*(row_in_f/0.07));
        //       pd_il_pa = pd_il_inh20*249;//pascal
               printf("\n pd il %f in ",pd_il_inh20);

               pd_de_inh20 = (1*((3*pow(10,-7)*vel_out_f*vel_out_f)-(9*pow(10,-6)*vel_out_f) + 0.003 )*(row_out_f/0.07));
        //       pd_de_pa = pd_de_inh20*249;//pascal
               printf("\n pd de %f in ",pd_de_inh20);



               pdf1 = (fh_f/fs_f);
               pdf2 =pow(pdf1,pd_fill_n);
               pdf3 = ((pd_fill_a)+(pd_fill_b *pow(fh_f,pd_fill_m)));
               pdf4 = (pd_fill_c * (pow(vel_fill_f,pd_fill_o)));
               pdf5 = (pd_fill_d * pow(wl_f,pd_fill_q));
               pdf6 = (pow(vel_fill_f,pd_fill_p));

               pd_fill_inh20 = ((pdf2 * pdf3 * pdf4 * pdf5 * pdf6)*(row_fill_f/0.07));
               printf("\n pd fill %f in ",pd_fill_inh20);
         //      pd_fill_pa = pd_fill_inh20*249;

               tpd_s_inh20 = (pd_il_inh20 + pd_de_inh20 + pd_fill_inh20);
               printf("\n hs %f in ",tpd_s_inh20);
               tpd_s_pa = tpd_s_inh20*249;

               pd_v_f =  (row_out_f/1202312)*(cfm_fan_f/edisa_f)*(cfm_fan_f/edisa_f);
               printf("\n hv %f in",pd_v_f);
               pd_v_s = pd_v_f * 249 ;
               tpd_f = pd_v_f + tpd_s_inh20;
               tpd_s = tpd_s_pa + pd_v_s;
               printf("\n the overall pressure drop across the cooling tower is %f pa or %f inches of water col",tpd_s,tpd_f);

               bhp = ((cfm_fan_f*(tpd_s_inh20+pd_v_f))/(6356*fann));//hp
               bhp_w = bhp*0.745699872;//kw
               printf("\n the required brake horse power of the fan is %f Bhp or %f kw",bhp,bhp_w);

return (bhp_w);
}


void kkt_optimize()
{
    float ofh,ofs,obhp,bhp_max,bhp_min,fill_h[40],fill_s[40],bh_p[40],bhp_dif[40];
    float fhi = 0.6096,/*fsi = 12.0f,fhf = 1.8288,*/ fsf = 19.0f,lam = 0.8;
    int p,i,j,q,z;
    int x0[10],xe[10],location = 0;
    float y0[10],ye[10];
    int gd = DETECT, gm;
     initgraph(&gd,&gm,"c:\\turbo\\tc\\bgi");

    fill_h[0] = fhi;
    fill_s[0] = fsf;
    bh_p[0] = bhpeqns(fill_h[0],fill_s[0]);
    bhp_dif[0] = 0;
    p=1;
    printf("at the initial condition of the fill height of %f m and  flute spacing of %f mm the fan poer is %f kw",fill_h[0],fill_s[0],bh_p[0]);
    q = 0;

    for(i = 1 ; i < 8 ; i++)
    {
    while(fill_s[i-1] > 11.0f)
    {
        p = q + 1 ;
        for(j = 1 ; lam > 0.002 ; j++)
        {
            while(fill_h[j-1] < 1.8289)
            {
            fill_h[j] = fill_h[j-1] + lam;
            fill_s[j] = fill_s[i-1];
            bh_p[j] = fabs(bhpeqns(fill_h[j],fill_s[j]));
            bhp_dif[j] = fabs(((bh_p[j] - bh_p[j-1])));
            if(bhp_dif[j] < 0.5)
            {
                fill_h[j] = fill_h[j-1] + (lam + 0.01);
                fill_s[j] = fill_s[i-1];
                bh_p[j] = fabs(bhpeqns(fill_h[j],fill_s[j]));
            }
            else
            {
                fill_h[j] = fill_h[i-1] + lam;
                fill_s[j] = fill_s[i-1];
                bh_p[j] = fabs(bhpeqns(fill_h[j],fill_s[j]));
            }
            fill_h[p] = fill_h[j];
            fill_s[p] = fill_s[j];
            bh_p[p] = bh_p[j];
            if(fill_h[j] > 1.4)
            {
                lam = lam/2;
                 j++;
            }
            else
            {
                lam = lam;
                j++;
            }
            p = p +1 ;   q = p;
           }
        }
        fill_s[i] = fill_s[i - 1] - 1;
        i++;
      }
    }
    z = q;
printf("\n the total number of iterations done are %d",q);
for(p = 0; p < z ; p++)
{
printf("\n %d  fh %f  fs  %f  bhp  %f",p,fill_h[p],fill_s[p],bh_p[p]);
}
bhp_max =  bh_p[0];
bhp_min =  bh_p[0];

for(i=1;i< z ;i++)
{
if(bhp_max < (bh_p[i]))
{   bhp_max =  bh_p[i]; }
if(bhp_min > (bh_p[i]))
{   bhp_min = bh_p[i]; location = i; }
}
printf("\n max is %f \n min is %f",bhp_max ,bhp_min);
printf("\n the optimum is at condition %d having  bhp %f at fh  %f  fs  %f",location,bh_p[location],fill_h[location],fill_s[location]);

getch();

}

i am facing the problem with kkt_optimize function

[Schol-R-LEA]

I'll assume from a previous related thread that you are using Turbo C++, which is unfortunate but I know there's nothing to be done about it (except protest to the government University Standards board, which I understand is probably futile). Therefore I'll limit my comments about the headers to a reminder that you don't actually need several of the ones you have included.

In fact, I am dismayed at how little you actually followed my earlier advice, in general. I do see a lot of improvements, but the code is still pretty much a disaster. However, if the bhpeqns() equations function now works, then that is what is most important, I suppose.

As for the ktt_optimize() function, what exactly is it doing, and where is it going wrong?

[sai.aleenus]

i also have some graphs to be generated for which i am using the graphics header and when i am removing any of thge other header files,the program is throwing an error..
bhpeqns is working and is also giving quite accurate results..i am basically a mechanical student and my grasp of c is not that extensive so do excuse the lengthy code..
well coming to point , the objective of optimize function is to evaluate bhp for varying values of fill_h and fill_s and the output format is to be in the form given in the initial thread of this discussion..
but it is evaluating for only one fill_h for a particular fill_s and if even the lam value is changed from 0.8 to 0.6,it shows error's of the kind(pow : error, / : error, etc).to be noted i have made all the base as well as exponents of every pow function positive in this program to negate the error..