rendered paste body#include <stdio.h>
#include <math.h>
#include <conio.h>
#define R 8.314
int main()
{
int comp,i,j,k,a,b;
double tempc,tempk,pressdp,pressbp,vi,vio,est,f,fo,df,**mtrx,pressorg,**invv,**tcrit,**zcrit,**wcrit,**vcrit,*gammab,tempb,tempd,*gammad,**pcrit,*phid,*phib,**bij,**delta,*x,*y,*ki,*z,*xo,*yo,*phi,pjsat,tempco,pjsato,*tsat,*gamma,*psat,frac=0.0,press=0.0,sum1=0.0,sum2=0.0,sum3=0.0,presso,*gammao,b0,b1;
FILE *input;
input=fopen("data.txt","r");
fscanf(input,"%d\n",&comp);
x=(double*)malloc(comp*sizeof(double));
y=(double*)malloc(comp*sizeof(double));
phi=(double*)malloc(comp*sizeof(double));
gamma=(double*)malloc(comp*sizeof(double));
psat=(double*)malloc(comp*sizeof(double));
tsat=(double*)malloc(comp*sizeof(double));
mtrx=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
mtrx[i]=(double *) malloc ((9+comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<9+comp;j++)
fscanf(input,"%lf",&mtrx[i][j]);
}
invv=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
invv[i]=(double *) malloc ((comp) * sizeof (double));
tcrit=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
tcrit[i]=(double *) malloc ((comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
tcrit[i][j]=pow((mtrx[i][1]*mtrx[j][1]),0.5);
}
zcrit=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
zcrit[i]=(double *) malloc ((comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
zcrit[i][j]=(mtrx[i][0]+mtrx[j][0])/2.0;
}
wcrit=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
wcrit[i]=(double *) malloc ((comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
wcrit[i][j]=(mtrx[i][3]+mtrx[j][3])/2.0;
}
vcrit=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
vcrit[i]=(double *) malloc ((comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
vcrit[i][j]=pow(((pow(mtrx[i][4],(1.0/3.0))+pow(mtrx[j][4],(1.0/3.0)))/2.0),3.0);
}
pcrit=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
pcrit[i]=(double *) malloc ((comp) * sizeof (double));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
pcrit[i][j]=(zcrit[i][j]*R*1000.0*tcrit[i][j])/vcrit[i][j];
}
bij=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
bij[i]=(double *) malloc ((comp) * sizeof (double));
delta=(double**)malloc(comp * sizeof (double *));
for(i=0;i<comp;i++)
delta[i]=(double *) malloc ((comp) * sizeof (double));
printf("*************************************************************************\n");
printf("* Bogazici University *\n");
printf("* Department of Chemical Engineering *\n");
printf("* ChE 321 - Introduction to Thermodynamics *\n");
printf("* Fall 2005 - Term Project *\n* *\n");
printf("* prepared by Ali Riza Tanriverdi - Murat Sezgin *\n");
printf("*************************************************************************\n\n");
printf("This program calculates Dew/Bubble Temperatures and Dew/Bubble Pressures for the given T/P and x/y values in addition the user can perform a flash calculation. The program can calculate these values by using three types of equations: Raoult's Law, Modified Raoult's Law and Gamma/Fi Calculation(the real case).\n");
printf("\n Which computation do you want to make ?\n");
printf(" (1) Bubble Pressure Calculation\n (2) Dew Pressure Calculation\n (3) Bubble Temperature Calculation\n (4) Dew Temperature Calculation\n (5) P,T Flash Calculation\n\n");
printf(" Please, enter your choice:");
scanf("%d",&a);
printf("\n Which correalation formula do you want to use ?\n");
printf(" (1) Raoult's Law\n (2) Modified Raoult's Law with Wilson Equation\n (3) Gamma-Fi Formulation\n\n");
printf(" Please, enter your choice:");
scanf("%d",&b);
switch(a)
{
case 1:
switch(b)
{
case 3:
//BUBBLE P
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[comp-1]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
printf("Pressure:%f kPa\n",press);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
case 2:
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[comp-1]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
printf("Pressure:%f kPa\n",press);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
case 1:
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[comp-1]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
printf("Pressure:%f kPa\n",press);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
}
break;
case 2:
switch(b)
{
case 3:
// DEW P
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
printf("Pressure is %f kPa\n",press);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
case 2:
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
printf("Pressure is %f kPa\n",press);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
case 1:
printf("Please, enter temperature(in Celcius):");
scanf("%lf",&tempc);
tempk=tempc+273.15;
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
printf("Pressure is %f kPa\n",press);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
}
break;
case 3:
switch(b)
{
case 3:
//BUBBLE T
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f\n",tempc);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
case 2:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f\n",tempc);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
case 1:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter x%d:",i+1);
scanf("%lf",&x[i]);
frac=x[i]+frac;
}
x[i]=1.0-frac;
for(i=0;i<comp;i++)
{
phi[i]=1.0;
gamma[i]=1.0;
}
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
sum1=0.0;
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f\n",tempc);
for(i=0;i<comp;i++)
printf("y%d:%f gamma%d:%f fi%d:%f\n",i+1,y[i],i+1,gamma[i],i+1,phi[i]);
break;
}break;
// DEW T
case 4:
switch(b)
{
case 3:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f C\n",tempc);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
case 2:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f C\n",tempc);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
case 1:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
for(i=0;i<comp-1;i++)
{
printf("Please, enter y%d:",i+1);
scanf("%lf",&y[i]);
frac=y[i]+frac;
}
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
printf("Temperature is %f C\n",tempc);
for(i=0;i<comp;i++)
printf("x%d:%f gamma%d:%f fi%d:%f\n",i+1,x[i],i+1,gamma[i],i+1,phi[i]);
break;
}
break;
//FLASH
case 5:
switch(b) {
case 3:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
pressorg=press;
z=(double*)malloc(comp * sizeof (double));
ki=(double*)malloc(comp * sizeof (double));
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter z%d:",i+1);
scanf("%lf",&z[i]);
frac=z[i]+frac;
}
z[i]=1.0-frac;
for(i=0;i<comp;i++)
{
y[i]=z[i];
}
//dew t
y[i]=1.0-frac;
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
tempd=tempc;
//DEW T
//buble t
for(i=0;i<comp;i++)
{
x[i]=z[i];
}
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
tempb=tempc;
printf("Tdew is %f C\n",tempd);
printf("Tbubble is %f C\n",tempb);
tempc=(tempb+tempd)/2.0;
printf("Temperature value is %f C\n",tempc);
tempk=tempc+273.15;
//bubble t END
frac=0.0;
for(i=0;i<comp;i++)
{
y[i]=z[i];
}
////dewp
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
pressdp=press;
phid=(double*)malloc(comp*sizeof(double));
gammad=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phid[i]=phi[i];
gammad[i]=gamma[i];
}
//bubble p
for(i=0;i<comp;i++)
{
x[i]=z[i];
}
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
pressbp=press;
phib=(double*)malloc(comp*sizeof(double));
gammab=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phib[i]=phi[i];
gammab[i]=gamma[i];
}
//bubble p END
printf("Pdew is %f kPa at %f K\n",pressdp,tempc);
printf("Pbubble is %f kPa at %f K\n\n",pressbp,tempc);
if(pressorg<pressdp||pressorg>pressbp)
{
printf("Your pressure value is not between Pdew and Pbuble!\n");
return -1;
}
vi=(pressbp-pressorg)/(pressbp-pressdp);
est=(pressorg-pressdp)/(pressbp-pressdp);
for(i=0;i<comp;i++)
{
phi[i]=(est*(phib[i]-phid[i]))+phid[i];
gamma[i]=(est*(gammab[i]-gammad[i]))+gammad[i];
}
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
// flash 1st loop
while(1)
{
for(i=0;i<comp;i++)
ki[i]=gamma[i]*psat[i]/(phi[i]*pressorg);
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
fo=f;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
vio=vi;
// calculating V by Newton's Method
while(1)
{
vi=vi-f/df;
fo=f;
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
if(fo-f<0.001&&fo-f>-0.001)
break;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
}
xo=(double*)malloc(comp*sizeof(double));
yo=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
xo[i]=x[i];
yo[i]=y[i];
x[i]=z[i]/(1.0+vi*(ki[i]-1.0));
y[i]=ki[i]*x[i];
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
b0=0.083-(0.422/(pow((tempk/tcrit[i][j]),(1.6))));
b1=0.139-(0.172/(pow((tempk/tcrit[i][j]),(4.2))));
bij[i][j]=((R*1000.0*tcrit[i][j])/pcrit[i][j])*(b0+wcrit[i][j]*b1);
}
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
delta[j][i]=2.0*bij[j][i]-bij[j][j]-bij[i][i];
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
{
for(k=0;k<comp;k++)
sum1+=y[j]*y[k]*(2.0*delta[j][i]-delta[j][k]);
}
phi[i]=exp((bij[i][i]*(press-psat[i])+0.5*press*sum1)/(R*1000.0*tempk));
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(xo[i]-x[i]<0.001&&xo[i]-x[i]>-0.001)
j++;
if(yo[i]-y[i]<0.001&&yo[i]-y[i]>-0.001)
j++;
}
if(vio-vi<0.001&&vio-vi>-0.001)
j++;
if(j>comp+comp)
break;
}
printf("Vi = %f\n",vi);
for(i=0;i<comp;i++)
{
printf("x%d:%f y%d:%f gamma%d:%f fi%d:%f \n",i+1,x[i],i+1,y[i],i+1,gamma[i],i+1,phi[i]);
}
break;
case 2:
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
pressorg=press;
z=(double*)malloc(comp * sizeof (double));
ki=(double*)malloc(comp * sizeof (double));
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter z%d:",i+1);
scanf("%lf",&z[i]);
frac=z[i]+frac;
}
z[i]=1.0-frac;
for(i=0;i<comp;i++)
{
y[i]=z[i];
}
//dew t
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*1.0/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
x[i]=y[i]*1.0*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*1.0/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*1.0*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
tempd=tempc;
//DEW T
//buble t
for(i=0;i<comp;i++)
x[i]=z[i];
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
tempb=tempc;
printf("Tdew is %f C\n",tempd);
printf("Tbubble is %f C\n",tempb);
tempc=(tempb+tempd)/2.0;
printf("Temperature value is %f C\n",tempc);
tempk=tempc+273.15;
//bubble t
frac=0.0;
for(i=0;i<comp;i++)
{
y[i]=z[i];
}
////dewp
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
pressdp=press;
phid=(double*)malloc(comp*sizeof(double));
gammad=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phid[i]=phi[i];
gammad[i]=gamma[i];
}
//bubble p
for(i=0;i<comp;i++)
x[i]=z[i];
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
pressbp=press;
phib=(double*)malloc(comp*sizeof(double));
gammab=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phib[i]=1.0;
gammab[i]=gamma[i];
}
//bubble p
printf("Pdew is %f kPa at %f K\n",pressdp,tempc);
printf("Pbubble is %f kPa at %f K\n\n",pressbp,tempc);
if(pressorg<pressdp||pressorg>pressbp)
{
printf("Your pressure value is not between Pdew and Pbuble!\n");
return -1;
}
vi=(pressbp-pressorg)/(pressbp-pressdp);
est=(pressorg-pressdp)/(pressbp-pressdp);
for(i=0;i<comp;i++)
{
phi[i]=1.0;
gamma[i]=(est*(gammab[i]-gammad[i]))+gammad[i];
}
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
// flash 1st loop
while(1)
{
for(i=0;i<comp;i++)
ki[i]=gamma[i]*psat[i]/(phi[i]*pressorg);
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
fo=f;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
vio=vi;
// calculating V by Newton's Method
while(1)
{
vi=vi-f/df;
fo=f;
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
if(fo-f<0.001&&fo-f>-0.001)
break;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
}
xo=(double*)malloc(comp*sizeof(double));
yo=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
xo[i]=x[i];
yo[i]=y[i];
x[i]=z[i]/(1.0+vi*(ki[i]-1.0));
y[i]=ki[i]*x[i];
}
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
invv[i][j]=(mtrx[j][5]/mtrx[i][5])*exp((-1.0*mtrx[i][j+9])/(R*tempk));
}
sum1=0.0;
for(i=0;i<comp;i++)
{
for(j=0;j<comp;j++)
sum1+=x[j]*invv[i][j];
sum2=0.0;
for(k=0;k<comp;k++)
{
sum3=0.0;
for(j=0;j<comp;j++)
sum3+=x[j]*invv[k][j];
sum2+=(x[k]*invv[k][i])/sum3;
}
gamma[i]=exp(1.0-log(sum1)-sum2);
sum1=0.0;
}
j=0;
for(i=0;i<comp;i++)
{
if(xo[i]-x[i]<0.001&&xo[i]-x[i]>-0.001)
j++;
if(yo[i]-y[i]<0.001&&yo[i]-y[i]>-0.001)
j++;
}
if(vio-vi<0.001&&vio-vi>-0.001)
j++;
if(j>comp+comp)
break;
}
printf("Vi = %f\n",vi);
for(i=0;i<comp;i++)
{
printf("x%d:%f y%d:%f gamma%d:%f fi%d:%f \n",i+1,x[i],i+1,y[i],i+1,gamma[i],i+1,phi[i]);
}
break;
case 1: /////////////////// //////////////////////////////////////////////////////////// case 11cse11111case11111111111111case111111
printf("Please, enter pressure(in kPa):");
scanf("%lf",&press);
pressorg=press;
z=(double*)malloc(comp * sizeof (double));
ki=(double*)malloc(comp * sizeof (double));
frac=0.0;
for(i=0;i<comp-1;i++)
{
printf("Please, enter z%d:",i+1);
scanf("%lf",&z[i]);
frac=z[i]+frac;
}
z[i]=1.0-frac;
for(i=0;i<comp;i++)
y[i]=z[i];
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=y[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
sum1=0.0;
pjsat=psat[0];
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
{
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
}
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
pjsat=psat[0];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=(y[i]*phi[i]/gamma[i])*(pjsat/psat[i]);
pjsat=press*sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
} tempd=tempc;
//DEW T
//buble t
for(i=0;i<comp;i++)
x[i]=z[i];
for(i=0;i<comp;i++)
{
phi[i]=1.0;
gamma[i]=1.0;
}
for(i=0;i<comp;i++)
tsat[i]=(-1.0*mtrx[i][8]+(mtrx[i][7]/(-1.0*log(press)+mtrx[i][6])));
tempc=0.0;
for(i=0;i<comp;i++)
tempc+=x[i]*tsat[i];
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
tempk=tempc+273.15;
pjsat=psat[0];//the first entry of the data file is selected for the calculation of Pjsat
sum1=0.0;
pjsato=psat[0];
for(i=0;i<comp;i++)
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
pjsat=press/sum1;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
while(1)
{
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
sum1=0.0;
for(i=0;i<comp;i++)
{
sum1+=(x[i]*gamma[i]/phi[i])*(psat[i]/pjsat);
}
pjsat=press/sum1;
tempco=tempc;
tempc=(-1.0*mtrx[0][8]+(mtrx[0][7]/(-1.0*log(pjsat)+mtrx[0][6])));
tempk=tempc+273.15;
if(tempco-tempc<0.1&&tempco-tempc>-0.1)
break;
}
tempb=tempc;
printf("Tdew is %f C\n",tempd);
printf("Tbubble is %f C\n",tempb);
tempc=(tempb+tempd)/2.0;
printf("Temperature value is %f C\n",tempc);
tempk=tempc+273.15;
//bubble t
frac=0.0;
for(i=0;i<comp;i++)
{
y[i]=z[i];
}
////dewp
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
while(1)
{
while(1)
{
for(i=0;i<comp;i++)
x[i]=y[i]*phi[i]*press/(gamma[i]*psat[i]);
gammao=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
gammao[i]=gamma[i];
sum1=0.0;
for(i=0;i<comp;i++)
sum1+=x[i];
for(i=0;i<comp;i++)
x[i]=x[i]/sum1;
j=0;
for(i=0;i<comp;i++)
{
if(gammao[i]-gamma[i]<0.1&&gammao[i]-gamma[i]>-0.1)
j++;
}
if(j==comp)
break;
}
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=((y[i]*phi[i])/(gamma[i]*psat[i]));
press=1.0/press;
if(press-presso<1.0&&presso-press>-1.0)
break;
}
pressdp=press;
phid=(double*)malloc(comp*sizeof(double));
gammad=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phid[i]=1.0;
gammad[i]=1.0;
gamma[i]=1.0;
phi[i]=1.0;
}
//bubble p
for(i=0;i<comp;i++)
x[i]=z[i];
for(i=0;i<comp;i++)
phi[i]=1.0;
for(i=0;i<comp;i++)
gamma[i]=1.0;
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
while(1)
{
for(i=0;i<comp;i++)
y[i]=x[i]*gamma[i]*psat[i]/(phi[i]*press);
presso=press;
press=0.0;
for(i=0;i<comp;i++)
press+=x[i]*gamma[i]*psat[i]/phi[i];
if(press-presso<0.1&&presso-press>-0.1)
break;
}
pressbp=press;
phib=(double*)malloc(comp*sizeof(double));
gammab=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
phib[i]=1.0;
gammab[i]=1.0;
gamma[i]=1.0;
}
//bubble p
printf("Pdew is %f kPa at %f C\n",pressdp,tempc);
printf("Pbubble is %f kPa at %f C,\n\n",pressbp,tempc);
if(pressorg<pressdp||pressorg>pressbp)
{
printf("Your pressure value is not between Pdew and Pbuble!\n");
return -1;
}
vi=(pressbp-pressorg)/(pressbp-pressdp);
est=(pressorg-pressdp)/(pressbp-pressdp);
for(i=0;i<comp;i++)
{
phi[i]=1.0;
gamma[i]=1.0;
}
for(i=0;i<comp;i++)
psat[i]=exp(mtrx[i][6]-(mtrx[i][7]/(tempc+mtrx[i][8])));
// flash 1st loop
while(1)
{
for(i=0;i<comp;i++)
ki[i]=gamma[i]*psat[i]/(phi[i]*pressorg);
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
fo=f;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
vio=vi;
// calculating V by Newton's Method
while(1)
{
vi=vi-f/df;
fo=f;
f=0.0;
for(i=0;i<comp;i++)
f+=(z[i]*(ki[i]-1.0))/(1.0+vi*(ki[i]-1.0));
if(fo-f<0.001&&fo-f>-0.001)
break;
df=0.0;
for(i=0;i<comp;i++)
df+=-1.0*(z[i]*(pow((ki[i]-1.0),2.0))/(pow((1.0+vi*(ki[i]-1.0)),2.0)));
}
xo=(double*)malloc(comp*sizeof(double));
yo=(double*)malloc(comp*sizeof(double));
for(i=0;i<comp;i++)
{
xo[i]=x[i];
yo[i]=y[i];
x[i]=z[i]/(1.0+vi*(ki[i]-1.0));
y[i]=ki[i]*x[i];
}
sum1=0.0;
j=0;
for(i=0;i<comp;i++)
{
if(xo[i]-x[i]<0.001&&xo[i]-x[i]>-0.001)
j++;
if(yo[i]-y[i]<0.001&&yo[i]-y[i]>-0.001)
j++;
}
if(vio-vi<0.001&&vio-vi>-0.001)
j++;
if(j>comp+comp)
break;
}
printf("Vi = %f\n",vi);
for(i=0;i<comp;i++)
{
printf("x%d:%f y%d:%f gamma%d:%f fi%d:%f \n",i+1,x[i],i+1,y[i],i+1,gamma[i],i+1,phi[i]);
}
break; }
break;
}
getch();
fclose(input);
return 7;
}