/* bwrm5.c */
#include <stdio.h>
#include <math.h>
#include "bwrm.h"


double
bwrmRTlnF(/*  Calculates RTln(f) of pure substances, where f is fugacity. */
  BWRM    *m,
  double   t,   /* temperature [K]  */
  double   rho /* molar density [kg-mole/m3] */
){
#define C	(*m).cf
#define MAT	(*m).tmp
  double f;

  RecalBWRM(m, t, rho);

  f  = -2.0*(MAT.c_T2 + MAT.g_T8 + MAT.h_T17 )/C.gamma
	*( 1.0 - ( 1.0 + MAT.gammarho2 + 0.5*MAT.gammarho2*MAT.gammarho2)*MAT.EXP1 ) ;
  f += 3.0*MAT.rho2*( MAT.c_T2 + MAT.g_T8 + MAT.h_T17 )
	*( 1.0/MAT.gammarho2 - ( 1.0/MAT.gammarho2 + 0.5 )*MAT.EXP1 );
  f += 6.0/5.0*MAT.rho5
	*( MAT.alphaa + MAT.alphad_T + MAT.alphae_T4 + MAT.alphaf_T23 );
  f += 3.0/2.0*MAT.rho2*( MAT.bRT - C.a - MAT.d_T - MAT.e_T4 - MAT.f_T23 );
  f += 2.0*rho*( MAT.B0RT - C.A0 - MAT.C0_T2 + MAT.D0_T3 - MAT.E0_T4 );
  f += Rgas*t*log( MAT.rhoRT );
  
  return(f);
#undef C
#undef MAT
}


double
bwrmF(/*  Calculates fugacity of pure substances  */
  BWRM    *coef,
  double   t,   /* temperature [K]  */
  double   rho  /* molar density [kg-mole/m3] */
){
  return( exp(bwrmRTlnF(coef, t, rho)/Rgas/t) );
}



double
bwrmF2(/*  Calculates fugacity of binary mixtures */
  BWRM2   *coef, /* BWR coefficients of mixtures */
  char     mode1,/* 
                   '1' or 1 :  component 1
                   '2' or 2 :  component 2  */
  char     mode2, /* 
                    'l' or 'L' :  calulates (*c).b) (liquid)
                    'v' or 'V' :  calulates (*c).b) (vapor),
                    else       :  calulates (*c).b) (bulk)   */
  double   t,    /* temperature [K]  */
  double   rho,  /* molar density [kg-mole/m3] */
  double   y     /* molar fraction [mol/mol]   */
){
  BWRM   *m, *p;
  struct bwrm_buffer *tmp_m, *tmp_p;
  double    RTlnf_x ;
  double    a,b,c,d,e,f,alpha,gamma,dpt,ept4,fpt23;
  double    buf;

  /* Component 1 or 2 ?? */
  if      (mode1 == '1' || mode1 == 1 ) p = &(coef->c[1]);
  else if (mode1 == '2' || mode1 == 2 ) p = &(coef->c[2]);
  else{
#ifndef __MS_C
	fprintf(stderr, "bwrmF_bi FATAL ERROR: mode1 must be '1' or '2'\n");
	fprintf(stderr, "                    : Check your program!\n");
	exit(1);
#else
	return(-1.0e+30);
#endif
  }
  /* bulk, vapor or liquid  ?? */
  m = &(coef->b);                                    /* bulk */
  if (mode2 == 'l' || mode2 == 'L') m = &(coef->l);  /* liquid */
  if (mode2 == 'v' || mode2 == 'V') m = &(coef->v);  /* vapor */
  tmp_m = &(m->tmp);
  tmp_p = &(p->tmp);

  if ( y == 0.0 ) return(0.0);


  RecalBWRM(p, t, rho);
  if (y != m->tmp.yprev){
	bwrmCalTmp(m, t, rho); /* cal. tmp by force */
  }else{
	RecalBWRM(m, t, rho);
  }

#define C_MX	(*m).cf
#define C_PU	(*p).cf
#define M_MX	(*tmp_m)
#define M_PU	(*tmp_p)
  a = pow( C_MX.a * C_MX.a * C_PU.a , (1.0/3.0));
  b = pow( C_MX.b * C_MX.b * C_PU.b , (1.0/3.0));
  c = pow( C_MX.c * C_MX.c * C_PU.c , (1.0/3.0));
  d = pow( C_MX.d * C_MX.d * C_PU.d , (1.0/3.0));
  e = pow( C_MX.e * C_MX.e * C_PU.e , (1.0/3.0));
  f = pow( C_MX.f * C_MX.f * C_PU.f , (1.0/3.0));
  alpha = pow(C_MX.alpha * C_MX.alpha * C_PU.alpha , (1.0/3.0));
  gamma = pow(C_MX.gamma * C_MX.gamma * C_PU.gamma , (1.0/3.0));
  
  buf = -2.0*(M_MX.c_T2 + M_MX.g_T8 + M_MX.h_T17)*sqrt(C_PU.gamma)/pow(C_MX.gamma, 1.5 );
  
  RTlnf_x = buf*(1.0 - (1.0 + M_MX.gammarho2 + 0.5*M_MX.gammarho2*M_MX.gammarho2)
                  * M_MX.EXP1 ) ;
  buf = M_MX.rho2
         * (3.0*c/M_MX.T2 + M_PU.g_T8 + 2.0*M_MX.g_T8 + M_PU.h_T17 + 2.0*M_MX.h_T17 );
  RTlnf_x += buf * ( 1.0/M_MX.gammarho2 - ( 1.0/M_MX.gammarho2 + 0.5 )*M_MX.EXP1 );
  
  RTlnf_x += 0.6*M_MX.rho5
	          * (alpha*( C_MX.a + C_MX.d/M_MX.T + C_MX.e/M_MX.T4 + C_MX.f/M_MX.T23));
  RTlnf_x += 0.6*M_MX.rho5*( C_MX.alpha*( a + d/M_MX.T + e/M_MX.T4 + f/M_MX.T23 ));
  RTlnf_x  += 1.5*M_MX.rho2*( b*Rgas*M_MX.T - a - d/M_MX.T - e/M_MX.T4 - f/M_MX.T23);

  RTlnf_x  += 2.0 * rho * (1.0-y)
             * (M_MX.B0_12RT - C_MX.A0_12 - M_MX.C0_12_T2 + M_MX.D0_12_T3 - M_MX.E0_12_T4);
  RTlnf_x  += 2.0*rho*y*( M_PU.B0RT - C_PU.A0 - M_PU.C0_T2 + M_PU.D0_T3 - M_PU.E0_T4);
  RTlnf_x  += Rgas*t*log( M_MX.rhoRT );
  
  return( exp( RTlnf_x/Rgas/t ) * y );
}
#undef C_MX
#undef C_PU
#undef M_MX
#undef M_PU