/* bwrm11.c by Satoru MOMOKI, Jan.1996 */
#include <stdio.h>
#include <math.h>
#include "bwrm.h"


double
bwrmYvYl_PT(/* 
      calculates yl[mol/mol]  & yv[mol/mol] from P[Pa] and T[K] of bibary 
      bibary mixtures at saturated condition. Return value is yv[mol/mol].
      And yl[mol/mol], yv,  rhol[kg-mol/m3] and  rhov are set on (*c).l.y[1], 
      (*c).v.y[1], (*c).prop_l.rho, and (*c).v.rho.  */
  BWRM2   *m,
  double   p,    /* pressure [Pa] */
  double   t,    /* temperature [K] */
  double   yv0,  /* initial value on y_v [mol/mol] */
  double   yl0   /* initial value on y_l [mol/mol] */
)
#define yv     X[0]
#define yl     X[1]
#define EPS0   1.0e-006  /* convergence limit  */
#define EPS1   1.0e-006  
#define DELTA0 1.0e-006
#define DELTA1 1.0e-006
{
  #define  N 2
  BWRM2    Xp[N], Xm[N]; /* Xp = x + DELTAXi, Xm = x - DELTAXi */
  double   A[N][N], B[N], X[N], dX[N];
  double   absDX[N], f;
  int      nit=0; /* nit : iterative count  */
  int      i;
  int      mode;
  #define MESSAGE (*m).b.tmp.errorcode

  if (yv0 > yl0) mode=1; else mode = 0;
  yv = yv0;  yl = yl0;
  (*m).v.rho = 0.0;    (*m).l.rho = 30.0;
  for(i = 0; i < N; i++)   Xp[i] = Xm[i] = *m;
  for(i = 0; i < N; i++)   dX[i] = 1.2 * EPS0;
  do{/* -----------------  BEGIN : Newton method  ------------------ */
	/* printf("%f %f %f\n", f, yv, yl);/**/
	CALC_F0F1_BWRM2(m,      p, t, yv,        yl );
	CALC_F0F1_BWRM2(&Xp[0], p, t, yv+DELTA0, yl );
	CALC_F0F1_BWRM2(&Xm[0], p, t, yv-DELTA0, yl );
	CALC_F0F1_BWRM2(&Xp[1], p, t, yv,        yl+DELTA1);
	CALC_F0F1_BWRM2(&Xm[1], p, t, yv,        yl-DELTA1);
	B[0] = F0_BWRM2(m);
 	B[1] = F1_BWRM2(m);
	for(i = 0; i < N; i++){
	  A[0][i] = (F0_BWRM2(&Xp[i]) - F0_BWRM2(&Xm[i])) / (2.0*DELTA0);/* dF0/dXi */
	  A[1][i] = (F1_BWRM2(&Xp[i]) - F1_BWRM2(&Xm[i])) / (2.0*DELTA1);/* dF1/dXi */
	}
	MESSAGE = bwrmSolveA22x_b(A, B, dX);
	if ( MESSAGE != 0 ){MESSAGE = BWRMERR_CANT_RESOLVE_MATRIX; break;}
	if ( nit++ > 100 ) {MESSAGE = BWRMERR_CANT_CONVERGE;       break;}

	/* cal f which is a factor to avoid y<0 or y>1 */
	f = 1.0;
	for(i=0; i < N; i++){
	  double z0, z1, zf;
	  zf = 1.0;
	  z0 = X[i];
	  z1 = X[i] - dX[i];;
	  if ( z1 <= 0.0) zf = (z0)/2.0       / (z0 - z1);
	  if ( z1 >= 1.0) zf = (1.0 - z0)/2.0 / (z1 - z0);
	  if ( zf < f) f = zf;
	}
	for(i=0; i < N; i++){
	  absDX[i] = fabs(dX[i]);
	  X[i] -= f * dX[i];
	}
	/*	printf("%f %f %f\n", f, yv, yl);/**/
	if (f < EPS0){ MESSAGE = BWRMERR_CANT_CONVERGE;       break;}
	/* if ((mode == 1 && yv <= yl)  || (mode == 0 && yv >= yl)){
	  double z;
	  z =yl; yl = yv; yv = z;
	}/**/
  }while( absDX[0] > EPS0 || absDX[1] > EPS1 );
  /* ---------------------  END : Newton method  -------------------- */
  if (!(0.0 <= yv && yv <= 1.0 && 0.0 <= yl && yl <= 1.0) ){
	MESSAGE = BWRMERR_CANT_CONVERGE;
  }
  switch(MESSAGE){
  case BWRMERR_CANT_CONVERGE:
	/* bwrm2ERROR(stdout, "bwrmYvYl_PT warning: not convergenced",p,t,yv,yl); /**/
	t = -1e+10;
	break;
  case BWRMERR_CANT_RESOLVE_MATRIX:
	/*
	  bwrm2ERROR(stderr, "Fatal error : bwrmYvYl_PT",  p, t, yv, yl);
	  bwrm2ERROR(stdout, "Fatal error : bwrmYvYl_PT",  p, t, yv, yl);
	/**/
	t = -1e+10;
	break;
  default:
	MESSAGE = BWRMERR_NOPROBLEM;
	CALC_F0F1_BWRM2(m, p, t, yv, yl); /* sets yv, yl, rhov,  rhol e.t.c. on 'm'  */
  }
  (*m).b.tmp.n_iteration = nit;
  return(p);
  #undef MESSAGE
}



int CALC_F0F1_BWRM2(
		BWRM2 *m,
		double P, double T, double yv1, double yl1
){
  (*m).v.y[1] = yv1;
  (*m).v.y[2] = 1.0-(yv1);
  if ((*m).v.rho < 0.0)         (*m).v.rho = 0.0000000001;
  if ((*m).v.rho > (*m).v.RHOc) (*m).v.rho = (*m).v.RHOc * 0.8;
  RecalBWRM2(m, 'v', T, (*m).v.rho, yv1);
  (*m).v.rho = bwrmRHO(&((*m).v), P, T, (*m).v.rho * 0.9);
  /* ----------------------- liquid ---------------------------  */
  (*m).l.y[1] = yl1;
  (*m).l.y[2] = 1.0-(yl1);
  if ((*m).l.rho < (*m).l.RHOc) (*m).l.rho = 30.0;/**/
  RecalBWRM2(m, 'l', T, (*m).l.rho, yl1);
  (*m).l.rho = bwrmRHO(&((*m).l), P, T, (*m).l.rho + 1.0);
  
  /*  printf("%5.3f %6.2f : yv%.5f yl%.5f rhov%6.4f rhol%6.3f\n",    
	 P/1.0e+6, T, yv1, yl1, (*m).v.rho, (*m).l.rho );/**/
  /* -------------------- f1v,f1l,f2v,f2l-------------------  */
  (*m).v.f[1] = bwrmF2(m, 1, 'v', T, (*m).v.rho, yv1);
  (*m).l.f[1] = bwrmF2(m, 1, 'l', T, (*m).l.rho, yl1);
  (*m).v.f[2] = bwrmF2(m, 2, 'v', T, (*m).v.rho, 1.0 - (yv1));
  (*m).l.f[2] = bwrmF2(m, 2, 'l', T, (*m).l.rho, 1.0 - (yl1));

  if ((*m).v.rho >= (*m).l.rho ) return -1;
  return 0;
}