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

#define MESSAGE (*m).b.tmp.errorcode




double
bwrmH2(/*
      calculate molar enthalpy 'h'[J/kg-mol] of binary mixtures.
      The function 'bwrmSetHid200SI' must be already done for BWRM strucure
	  of each component,c1 and c2, respectively. */
  BWRM2 *m,
  char  mode, /* 'l' or 'L': BWRM structure,l for liquid phase is used.
                 'v' or 'V': BWRM structure,v for liquid phase is used.
                  else     : BWRM structure,b for bulk fluid   is used.  */
  double t,   /* [K] */
  double rho, /* [kg-mol/m3] */
  double y1   /* molar fraction of component 1 */
){
  double h1_id, h2_id, h_id;
  BWRM *cf;
  cf = &((*m).b);                                  /* bulk */
  if (mode == 'l' || mode == 'L') cf = &((*m).l);  /* liquid */
  if (mode == 'v' || mode == 'V') cf = &((*m).v);  /* vapor */
  
  RecalBWRM2(m, mode, t, rho, y1);

  h1_id = (*m).c[1].h_id(t) + (*m).c[1].h_adj;
  h2_id = (*m).c[2].h_id(t) + (*m).c[2].h_adj;
  h_id = y1 * h1_id  +  (1.0-y1) * h2_id;

  (*cf).h = h_id + bwrmdH(cf, t, rho);
  return( (*cf).h  );
}




























double
bwrmT_Phy(/*
      Determines temperature of binary mixtures from pressure,P[Pa]
	  and molar  enthalpy 'h'[J/kg-mol].   And the quality is set on 
	  (*m).b.x.  You have to set  (*m)c[12].h_id(t) and (*m)c[12].h_adj with  
      bwrmSetHid200SI() before you use this function.  */
  BWRM2   *m,   /* structure of BWR coefficients */
  double  p,    /* pressure [Pa] */
  double  h,    /* molar enthalpy [J/kg-mol]        */
  double  yb)   /* bulk molar fraction [mol/mol]    */
{
  int    n;
  double x, T, dT, dhdT, yv, yl, F;
  double Tv, Tl, hv, hl; /* Tv_dew, Tl_bubble, hv_dew, hl_bubble */
  #define EPS  1e-5

  MESSAGE = BWRMERR_NOPROBLEM;
  n = 0; (*m).b.y[1] = yb;

  /* if (yb <= 1e-6 || yb >= 1.0 - 1e-6) goto RETURN;   /* pure substance */

  /* Liquid ?   */
  Tl = bwrmTyv_Pyl(m, p, yb);  hl = (*m).l.h = bwrmH2(m, 'l', Tl, (*m).l.rho, yb);
  if (h == hl){T = Tl; goto RETURN;}
  if (h <  hl){T = bwrmT_Ph2(m, 'l', p, h, yb); goto RETURN;}

  /* Gas ?      */
  Tv = bwrmTyl_Pyv(m, p, yb);  hv = (*m).v.h = bwrmH2(m, 'v', Tv, (*m).v.rho, yb);
  if (h == hv){T = Tv; goto RETURN;}
  if (h >  hv){T = bwrmT_Ph2(m, 'v', p, h, yb); goto RETURN;}

  /* Saturated  */
  x = (*m).b.x = (h - hl) / (hv - hl);
  T = (*m).b.T = x * Tv + (1.0 - x) * Tl;
  if ((*m).c[1].Tc  < (*m).c[2].Tc){
	yl = 0.9 * yb;                yv = yb + (1.0 - yb) * 0.3;
  }else{
	yl = yb + (1.0 - yb) * 0.3;   yv = 0.9 * yb;
  }	
  do{
	/* begin -- calculation of h and x */
	bwrmYvYl_PT(m,  p, T, yv, yl); yv = (*m).v.y[1]; yl = (*m).l.y[1];
	if (yb == 0.0 || yb == 1.0)	  goto RETURN;
	x = (yl - yb)/ (yl - yv);
	(*m).b.x = x;
	(*m).v.h = bwrmH2(m, 'v', T, (*m).v.rho, (*m).v.y[1]);
	(*m).l.h = bwrmH2(m, 'l', T, (*m).l.rho, (*m).l.y[1]);
	(*m).b.h = x*(*m).v.h + (1.0-x)*(*m).l.h;
	/* end   -- calculation of h and x */
	F = (*m).b.h - h;
	dhdT = 0.5*( ((*m).b.h  - hl)/(T - Tl) + (hv - (*m).b.h)/(Tv - T) );
	dT = F / dhdT;
	T -= dT;
	if (T <= Tl) {T += dT; T = Tl + (T - Tl) * 0.5;}
	if (T >= Tv) {T += dT; T = Tv - (Tv - T) * 0.5;}
	if (n++ > 100){
  /* printf("bwrmT_Ph2 warning: not convergenced(h=%lf[kJ/kg])\n", h/1.0e+3/(*m).b.M);/**/
	  MESSAGE = BWRMERR_CANT_CONVERGE ;
	  T = -1e+10;
	  goto RETURN;
	}
  }while(fabs(dT) > EPS );
  bwrmYvYl_PT(m,  p, T, yv, yl); yv = (*m).v.y[1]; yl = (*m).l.y[1];
  (*m).v.h = bwrmH2(m, 'v', T, (*m).v.rho, (*m).v.y[1]);
  (*m).l.h = bwrmH2(m, 'l', T, (*m).l.rho, (*m).l.y[1]);
  (*m).b.x = (yl -  yb) / (yl - yv);

 RETURN:
  (*m).b.T = T;   (*m).b.tmp.n_iteration = n;
  return(T);
}
#undef EPS



double
bwrmT_Ph2(/* 
      Determines temperature of binary mixtures from  molar enthalpy,
	  'h'[J/kg-mol].   You have to set  (*c).h_id(t) and  (*c).h_adj
      with  bwrmSetHid200SI() before you use this function.  */
  BWRM2  *m,    /* structure of BWR coefficients */
  char    mode,/* 
                 'l' or 'L' :  calulates (*c).b) (liquid)
                 'v' or 'V' :  calulates (*c).b) (vapor) */
  double  p,    /* pressure [Pa] */
  double  h,    /* molar enthalpy [J/kg-mol]     */
  double  y)    /* molar fraction [mol/mol]    */
#define LIQUID  (mode == 'l' || mode == 'L') 
#define VAPOR   (mode == 'v' || mode == 'V') 
{
  double t, dt;
  double rho, rhop, rhom;
  int    n = 0;
  BWRM *cf;
  #define EPS     1e-4
  #define D       1e-3

  if (mode == 'l' || mode == 'L') cf = &((*m).l);  /* liquid */
  if (mode == 'v' || mode == 'V') cf = &((*m).v);  /* vapor */

  if( LIQUID ){
   	cf = &((*m).l); rho = 30.0;    t = T0bwrm -10.0;
  }else if( VAPOR ){
   	cf = &((*m).v); rho = 0.00001; t = y * (*m).c[1].Tc + (1.0-y)*(*m).c[1].Tc;
  }else{
	/*	fprintf(stderr, "bwrmT_Ph2 FATAL ERROR: mode1 must be 'l' or 'v'\n");
	fprintf(stderr, "                     : Check your program!\n");
	exit(1);/**/
	(*cf).tmp.errorcode = BWRMERR_CANT_CONVERGE ;
	return -1.0;
  }
  do{
	if (LIQUID)	  rho  = bwrmRHO(cf, p, t, 1.1 * rho);
	else          rho  = bwrmRHO(cf, p, t, 0.9 * rho);
	rhop = bwrmRHO(cf, p, t+D, rho);
	rhom = bwrmRHO(cf, p, t-D, rho);
	RecalBWRM2(m, mode, t, rho, y);

	dt = (bwrmH2(m, mode, t+D, rhop, y) - bwrmH2(m, mode, t-D, rhom, y)) / (2.0*D);
	if (dt == 0.0) dt = 1e+10;
	dt = (bwrmH2(m, mode, t, rho, y) - h) / dt;

	if (fabs(dt) > 10.0)
	  if (dt > 0.0) dt = 10.0;else dt = -10.0; 
	t -= dt;
	if (n++ > 100){
	  /* printf("bwrmT_Ph2 warning: not convergenced(h=%lf[kJ/kg])\n",
			 h/1.0e+3 / (*cf).M);/**/
	  (*cf).tmp.errorcode = BWRMERR_CANT_CONVERGE ;
	  goto RETURN;
	}
  }while(fabs(dt) > EPS );
  (*cf).tmp.errorcode = BWRMERR_NOPROBLEM;
  goto RETURN;
  
 RETURN:
  (*cf).rho = rho;
  if( LIQUID ) (*cf).rhol = rho;
  else         (*cf).rhov = rho;
  (*cf).T = t;
  (*cf).tmp.n_iteration = n;
  return((*cf).T);
}
#undef EPS