/* bwrm1.c */

#include <math.h>
#include "bwrm.h"

/* -------------------------------------------------------------- */	
int
CalBWRM(/* 
       Calculates BWR coefficientsfor pure substance, 'c'.
*/
  BWRM    *c,
  double   tc,   /* critical temperature [K]  */
  double   rhoc, /* molar density [kg-mol/m3] */
  double   M,    /* molecular weight [kg/kg-mole] */
  double   w     /* omega */
)
{
#define		CF	(*c).cf
  double rhoc2, R_rhoc, R_rhoc2, tc2, tc3, tc4, tc5, tc9, tc18, w2;

  w2    = w * w ;
  rhoc2 = rhoc * rhoc ;
  R_rhoc  = Rgas / rhoc;
  R_rhoc2 = Rgas / rhoc2;
  tc2 = tc*tc;    tc3 = tc2*tc;   tc4 = tc3*tc;
  tc5 = tc4*tc;   tc9 = tc5*tc4;  tc18 = tc9 * tc9;
  
  CF.a     = (0.484011  + 0.754130*w )* tc * R_rhoc2;
  CF.b     = (0.528629  + 0.349261*w ) / rhoc2 ;
  CF.c     = (0.504087  + 1.32245 *w )* tc3 * R_rhoc2 ;
  CF.d     = (0.0732828 + 0.463492*w )* tc2 * R_rhoc2 ;
  CF.e     = (4.655930E-3 - 3.07393E-2*w + 5.58125E-2 * w2
              -3.40721E-3 * exp(-7.72753*w-45.3152*w2)) * tc5 * R_rhoc2;
  CF.f     = (0.697E-13 + 8.08E-13*w - 16.0E-13 * w2
              -0.363078E-13 * exp(30.9009*w-283.680*w2)) * (tc18*tc5*tc) * R_rhoc2;
  CF.g     = (2.20E-5 - 1.065E-4 * w 
				+1.09E-5 * exp(-26.024*w)) * tc9 * R_rhoc2 ;
  CF.h     = (-2.40E-11 + 11.8E-11 *w
				-2.05E-11 * exp(-21.52*w) ) * tc18 * R_rhoc2 ;
  CF.alpha = (0.0705233 - 0.044448*w ) / (rhoc2 * rhoc);
  CF.gamma = (0.544979  - 0.270896*w ) / rhoc2;
  
/*   CalBWRM_second( c, 1.0, tc, rhoc, w, tc, rhoc, w ); */
  CF.A0    = (1.28438   - 0.920731*w ) * tc * R_rhoc;
  CF.B0    = (0.443690  + 0.115449*w ) / rhoc;
  CF.C0    = (0.356306  + 1.70871 *w ) * tc3 * R_rhoc;
  CF.D0    = (0.0307452 + 0.179433*w ) * tc4 * R_rhoc;
  CF.E0    = (0.006450  - 0.022143*w*exp(-3.8*w) ) * tc5 * R_rhoc;

  (*c).n = 1;
  CF.A0_12 = 1.0;
  CF.B0_12 = 1.0;
  CF.C0_12 = 1.0;
  CF.D0_12 = 1.0;
  CF.E0_12 = 1.0;
  
  (*c).tmp.yprev = 1.0;
  (*c).tmp.Tprev = -1.0;
  (*c).tmp.rhoprev = -1.0;
  (*c).Pc = bwrmP(c, tc, rhoc);
  return 0;
#undef	CF
}









/* -------------------------------------------------------------- */	
int 
CalBWRM2(/* 
       Calculates BWR cficients'c' with consideration of composition.
       Caluculation of 'c.c1' and 'c.c2' by CalBWRM() and definetion of 
       c.mij are required before you call this function.
*/
  BWRM2   *c,
  char     mode, /* 
                       'l' or 'L' :  calulates (*c).b) (liquid)
                       'v' or 'V' :  calulates (*c).b) (vapor),
                       else       :  calulates (*c).b) (bulk)   */
  double   y1    /* molar fraction of component 1(more volatile) */
){
  double   y2, y1_2, y2_2, y1_y2mp2;
  double   m;
  double   rhoc, w, tc;
  double   tc3, tc4, tc5, R_rhoc;
  BWRM     *c0;
#define	CM      (*c0).cf
#define	C(i)    (*c).c[i].cf
#define RHOC(i) (*c).c[i].RHOc
#define TC(i)   (*c).c[i].Tc
#define W(i)    (*c).c[i].omega

  c0 = &((*c).b);                                  /* bulk */
  if (mode == 'l' || mode == 'L') c0 = &((*c).l);  /* liquid */
  if (mode == 'v' || mode == 'V') c0 = &((*c).v);  /* vapor */


  y2 = 1.0 - y1;  y2_2 = y2*y2;  y1_2 = y1*y1;  y1_y2mp2 = 2.0*y1*y2;

  (*c0).n = 2;
  (*c0).M = y1 * (*c).c[1].M + (1.0 - y1) * (*c).c[2].M ;

  CM.a = pow(y1*pow(C(1).a, (1.0/3.0) )  +  y2*pow(C(2).a , (1.0/3.0) ), 3.0);
  CM.b = pow(y1*pow(C(1).b, (1.0/3.0) )  +  y2*pow(C(2).b , (1.0/3.0) ), 3.0);
  CM.c = pow(y1*pow(C(1).c, (1.0/3.0) )  +  y2*pow(C(2).c , (1.0/3.0) ), 3.0);
  CM.d = pow(y1*pow(C(1).d, (1.0/3.0) )  +  y2*pow(C(2).d , (1.0/3.0) ), 3.0);
  CM.e = pow(y1*pow(C(1).e, (1.0/3.0) )  +  y2*pow(C(2).e , (1.0/3.0) ), 3.0);
  CM.f = pow(y1*pow(C(1).f, (1.0/3.0) )  +  y2*pow(C(2).f , (1.0/3.0) ), 3.0);
  CM.g = y1*C(1).g + y2*C(2).g;
  CM.h = y1*C(1).h + y2*C(2).h;
  CM.alpha = pow( y1*pow(C(1).alpha, (1.0/3.0)) + y2*pow(C(2).alpha, (1.0/3.0)), 3.0);
  CM.gamma = pow( y1*sqrt(C(1).gamma) + y2*sqrt(C(2).gamma ), 2.0);
  
  m =  ((*c).mij)(y1, 1.0);
  rhoc = pow((0.5 * (pow(RHOC(1), -(1.0/3.0)) + pow(RHOC(2), -(1.0/3.0)))), -3.0);
  tc   = m * sqrt(TC(1) * TC(2));
  w    = (W(1) + W(2)) / 2.0;
  tc3 = tc * tc * tc; tc4 = tc3 * tc; tc5 = tc4 * tc;
  R_rhoc = Rgas / rhoc;
  CM.A0_12 =  (1.28438   - 0.920731*w ) * tc * R_rhoc;
  CM.B0_12 =  (0.443690  + 0.115449*w ) / rhoc;
  CM.C0_12 =  (0.356306  + 1.70871 *w ) * tc3 * R_rhoc;
  CM.D0_12 =  (0.0307452 + 0.179433*w ) * tc4 * R_rhoc;
  CM.E0_12 =  (0.006450  - 0.022143*w*exp(-3.8*w) ) * tc5 * R_rhoc;

  CM.A0 = y1_2*C(1).A0 + y2_2*C(2).A0 + y1_y2mp2*CM.A0_12;
  CM.B0 = y1_2*C(1).B0 + y2_2*C(2).B0 + y1_y2mp2*CM.B0_12;
  CM.C0 = y1_2*C(1).C0 + y2_2*C(2).C0 + y1_y2mp2*CM.C0_12;
  CM.D0 = y1_2*C(1).D0 + y2_2*C(2).D0 + y1_y2mp2*CM.D0_12;
  CM.E0 = y1_2*C(1).E0 + y2_2*C(2).E0 + y1_y2mp2*CM.E0_12;


  (*c0).RHOc = rhoc; /* not exact, just for convenience */
  (*c0).Tc  = tc;    /* not exact, just for convenience */
  (*c0).tmp.yprev   =   y1;
  (*c0).tmp.Tprev   = -1.0;
  (*c0).tmp.rhoprev = -1.0;
  return 0;
#undef		CM
#undef		C
}






/* -------------------------------------------------------------- */	
int
InitBWRM(/* 
      Initialize BWRM structure for pure substance, 'c'.
*/
  BWRM    *c,
  double   tc,   /* critical temperature [K]  */
  double   rhoc, /* molar density [kg-mol/m3] */
  double   M,    /* molecular weight [kg/kg-mole] */
  double (*w_f)(double t),
		         /* function to calc. omega from T.
				    If omwga is constant, w_f = bwrmOMEGA_CONST. */
  double   w,    /* omega when w_f = bwrmOMEGA_CONST. */
  double   t,    /* temperature used for w_f. */
  double (*cp_id)(double t), /* function to calc. cp at ideal gas [J/kg-mol K]*/
  double (*h_id)(double t),  /* function to calc. h  at ideal gas [J/kg-mol]  */
  double (*s_id)(double t)   /* function to calc. s  at ideal gas [J/kg-mol K]*/
){
  (*c).n     = 1;
  (*c).Tc    = tc;
  (*c).M     = M;
  (*c).RHOc  = rhoc;

  (*c).omega = w;
  (*c).omega_f = w_f;
  if (w_f != bwrmOMEGA_CONST) (*c).omega = w_f(t);
  CalBWRM(c, tc, rhoc, M, (*c).omega);
  (*c).cp_id = cp_id;
  bwrmSetHid200SI(c, h_id);
  bwrmSetSid1SI(c, s_id);
/*   (*c).s_id = s_id; /* should make bwrmSetSid1SI(), later */
}
double bwrmOMEGA_CONST(double t){return 1.0;};
/* -------------------------------------------------------------- */	
int 
InitBWRM2(/* 
       Initializes BWR coefficients for binary mixtures, 'mx'.
       'pu1' and 'pu2' are required to be already caluculated
       using CalBWRM().
*/
  BWRM2  *mx,
  BWRM   *pu1,
  BWRM   *pu2,
  double  y1,
  double (*mij)(double y, double t)
)
{
  /* BWRMcpy(&((*mx).c[1]), pu1);  BWRMcpy(&((*mx).c[2]), pu2);/**/
  mx->c[1] = *pu1;
  mx->c[2] = *pu2;
  (*mx).c[1].tmp.yprev = 1.0;  (*mx).c[1].tmp.Tprev = -1.0;  (*mx).c[1].tmp.rhoprev = -1.0;
  (*mx).c[2].tmp.yprev = 1.0;  (*mx).c[2].tmp.Tprev = -1.0;  (*mx).c[2].tmp.rhoprev = -1.0;
  (*mx).mij = mij;
  CalBWRM2(mx, 'b', y1 );
  /* BWRMcpy(&((*mx).v), &((*mx).b));  BWRMcpy(&((*mx).l), &((*mx).b));/**/
  mx->v = mx->l = mx->b;
  return 0;
}












/* -------------------------------------------------------------- */	
int
bwrmCalTmp(/* Calculate c.tmp */
  BWRM *c,
  double    t,     /* [K] */
  double    rho    /* [kg-mole/m3] */
){
  bwrmCalTmpT(c, t);
  bwrmCalTmpRHO(c, rho);
  bwrmCalTmpOther(c, t, rho);
}

/* -------------------------------------------------------------- */	
int
bwrmCalTmpT(/* Calculate T^2, T^4, and so on,
               and substitute them into 'c.tmp'. */
  BWRM *c,
  double    t      /* [K] */
){
#define MAT	(*c).tmp
  MAT.Tprev = t;
  MAT.T   = t;          MAT.T2  = t*t;     MAT.T3   = MAT.T2*t ;
  MAT.T4  = MAT.T3*t;     MAT.T5  = MAT.T4*t;  MAT.T6   = MAT.T5*t ;
  MAT.T7  = MAT.T6*t;     MAT.T8  = MAT.T7*t;  MAT.T9   = MAT.T8*t ;
  MAT.T16 = MAT.T8*MAT.T8;  MAT.T17 = MAT.T16*t; MAT.T18  = MAT.T17*t;
  MAT.T22 = MAT.T16*MAT.T6; MAT.T23 = MAT.T22*t; MAT.T24  = MAT.T23*t;
#undef MAT
}

/* -------------------------------------------------------------- */	
int
bwrmCalTmpRHO(/* Calculate rho^2, rho^4, and so on,
                 and substitute them into 'c.tmp'. */
  BWRM *c,
  double    rho    /* [kg-mole/m3] */
){
#define MAT	(*c).tmp
  MAT.rhoprev = rho;
  MAT.rho  = rho;          MAT.rho2 = rho*rho;      MAT.rho3 = MAT.rho2*rho ;
  MAT.rho4 = MAT.rho3*rho; MAT.rho5 = MAT.rho4*rho; MAT.rho6 = MAT.rho5*rho ;
  MAT.rho7 = MAT.rho6*rho;
#undef MAT
}


/* -------------------------------------------------------------- */	
int
bwrmCalTmpOther(/* Calculate rho*R*T, B0*R*, and so on,
                  and substitute them into 'c.tmp'. */
  BWRM     *c,
  double    t,     /* [K] */
  double    rho    /* [kg-mole/m3] */
){
#define	C   (*c).cf
#define MAT (*c).tmp
  MAT.rho = rho ;
  MAT.rho2 = rho * rho ;
  MAT.rhoRT = MAT.rho * Rgas * MAT.T; MAT.B0RT  = C.B0 * Rgas * MAT.T;
  MAT.C0_T2 = C.C0 / MAT.T2;          MAT.D0_T3 = C.D0 / MAT.T3;
  MAT.E0_T4 = C.E0 / MAT.T4;          MAT.bRT   = C.b * Rgas * MAT.T;
  MAT.c_T2  = C.c / MAT.T2;           MAT.d_T   = C.d / MAT.T;
  MAT.e_T4  = C.e / MAT.T4;           MAT.f_T23 = C.f / MAT.T23;
  MAT.g_T8  = C.g / MAT.T8;           MAT.h_T17 = C.h / MAT.T17;
  MAT.alphaa= C.alpha * C.a;          MAT.alphad_T = C.alpha * MAT.d_T;
  MAT.alphae_T4 = C.alpha * MAT.e_T4; MAT.alphaf_T23= C.alpha * MAT.f_T23;
  MAT.gammarho2 = MAT.rho2 * C.gamma;
  MAT.EXP1      = exp( -C.gamma * MAT.rho2 );
  if ((*c).n == 2){
	MAT.B0_12RT  = C.B0_12 * Rgas * MAT.T	;
	MAT.C0_12_T2 = C.C0_12 / MAT.T2;
	MAT.D0_12_T3 = C.D0_12 / MAT.T3;
	MAT.E0_12_T4 = C.E0_12 / MAT.T4;
  }
#undef	C
#undef	MAT
}
/* -------------------------------------------------------------- */	
int BWRM2cpy(/* c'a' <- c'b'  */   
  BWRM2 *a,
  BWRM2 *b )
{
  int i;
  char  *p1, *p2;
  p1 = (char *)a;
  p2 = (char *)b;
  for(i = 0; i < sizeof(BWRM2); i++){
	*(p1+i) = *(p2+i);
  }
}
int BWRMcpy(/* c'a' <- c'b'  */   
  BWRM *a,
  BWRM *b )
{
  int i;
  char  *p1, *p2;
  p1 = (char *)a;
  p2 = (char *)b;
  for(i = 0; i < sizeof(BWRM); i++){
	*(p1+i) = *(p2+i);
  }
}


int  BWRctmpcpy(/* copy c.tmp of 'b' to 'a' */
  BWRM *a,
  BWRM *b )
{
  int i;
  char  *p1, *p2;
  p1 = (char *)a;
  p2 = (char *)b;
  for(i = 0; i < sizeof(struct bwrm_buffer); i++){
	*(p1+i) = *(p2+i);
  }
}