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

#define MESSAGE (*m).b.tmp.errorcode
#define D    10.0
#define EPS  1.0e-1

double
bwrmPyv_Tyl(/* 
      Calculates P[Pa]  & yv[mol/mol] from T[K] and yl[mol/mol] of binary 
      bibary mixtures at saturated condition. Return value is P[Pa].
      And yv[mol/mol], rhol[kg-mol/m3] and  rhov are set on (*c).v.y[1], 
      (*c).prop_l.rho, and (*c).v.rho.
	  Results with yb as yl indicate properties at bubble point.
	  */
  BWRM2   *m,
  double   t,   /* Temperature[K] */
  double   yl   /* molar fraction in liquid  [mol/mol] */
)
{
  BWRM2 mp, mm;
  double p, dp, omega, tc, pc, yv;
  double F, Tp, Tm, dP, dF, dTdP;
  int n = 0;
  
  if (t > m->c[1].Tc && t > m->c[2].Tc){
	  MESSAGE = BWRMERR_ARGUMENT_OUTOFRANGE;
	  return -1.0;
  }

  mp = mm = *m;
  MESSAGE = BWRMERR_NOPROBLEM;  
  omega = yl * (*m).c[1].omega  + (1-yl) * (*m).c[2].omega;
  pc    = yl * (*m).c[1].Pc     + (1-yl) * (*m).c[2].Pc;
  tc    = yl * (*m).c[1].Tc     + (1-yl) * (*m).c[2].Tc;


  /* define initial value of P */
  p = pc * pow(10.0, -7.0/3.0 * (1.0+omega) * (tc/t-1.0));

  /* define initial value of yv */
  if ((*m).c[1].Tc  < (*m).c[2].Tc)  yv = yl + (1.0 -yl) * 0.2;
  else                               yv = 0.9 * yl;
  mp.v.y[1] = mm.v.y[1] = yv;

  do{
	Tp = bwrmTyv_Pyl2(&mp, p+D, yl, t, mp.v.y[1]);
	Tm = bwrmTyv_Pyl2(&mp, p-D, yl, t, mm.v.y[1]);
	dTdP = (Tp - Tm)/(2.0*D);
	if (dP <= 0.0) dP = 1e+8;
	F = (Tp + Tm)/2 - t;
	dP = F / dTdP;
	p -= dP;
	/* printf("%e %e %f %f %f %f %f\n", p, dP, Tp, Tm, omega, t, yl);/**/
	if (n++ > 1000){
	  /* 
	  printf("bwrmPyv_Tyl warning: not converged(T=%lfK)\n", t);
	  /**/
	  MESSAGE = BWRMERR_CANT_CONVERGE ;
	  break;
	}
  }while(fabs(dP) > EPS );
  (*m).b.tmp.n_iteration = n;
  yv = (mp.v.y[1] + mm.v.y[1])/2;
  bwrmTyv_Pyl2(m, p, yl, t, yv);
  (*m).b.P = p;
  return((*m).b.P);
}



double
bwrmPyl_Tyv(/* 
      Calculates P[Pa]  & yv[mol/mol] from T[K] and yv[mol/mol] of binary 
      bibary mixtures at saturated condition. Return value is P[Pa].
      Also,  yl[mol/mol], rhol[kg-mol/m3] and  rhov are set on (*c).l.y[1], 
      (*m).l.rho, and (*m).v.rho.
	  Results with yb as yv indicate properties at dew point.
	  */
  BWRM2   *m,
  double   t,   /* Temperature[K] */
  double   yv   /* molar fraction in liquid  [mol/mol] */
)
{
  BWRM2 mp, mm;
  double p, dp, omega, tc, pc, yl;
  double F, Tp, Tm, dP, dF, dTdP;
  int n = 0;
  


  if (t > m->c[1].Tc && t > m->c[2].Tc){
	  MESSAGE = BWRMERR_ARGUMENT_OUTOFRANGE;
	  return -1.0;
  }


  mp = mm = *m;
  MESSAGE = BWRMERR_NOPROBLEM;  
  omega = yv * (*m).c[1].omega  + (1-yv) * (*m).c[2].omega;
  pc    = yv * (*m).c[1].Pc     + (1-yv) * (*m).c[2].Pc;
  tc    = yv * (*m).c[1].Tc     + (1-yv) * (*m).c[2].Tc;

  /* define initial value of P */
  p = pc * pow(10.0, -7.0/3.0 * (1.0+omega) * (tc/t-1.0));

  /* define initial value of yv */
  if ((*m).c[1].Tc  > (*m).c[2].Tc)  yl = yv + (1.0 -yv) * 0.2;
  else                               yl = 0.9 * yv;
  mp.l.y[1] = mm.l.y[1] = yl;

  do{
	Tp = bwrmTyl_Pyv2(&mp, p+D, yv, t, mp.l.y[1]);
	Tm = bwrmTyl_Pyv2(&mp, p-D, yv, t, mm.l.y[1]);
	dTdP = (Tp - Tm)/(2.0*D);
	if (dP <= 0.0) dP = 1e+8;
	F = (Tp + Tm)/2 - t;
	dP = F / dTdP;
	p -= dP;
	if (n++ > 1000){
	  /*
	  printf("bwrmPyl_Tyv warning: not converged(T=%lfK)\n", t);
	  /**/
	  MESSAGE = BWRMERR_CANT_CONVERGE ;
	  break;
	}
  }while(fabs(dP) > EPS );
  (*m).b.tmp.n_iteration = n;
  yl = (mp.l.y[1] + mm.l.y[1])/2;
  bwrmTyl_Pyv2(m, p, yv, t, yl);
  (*m).b.P= p;
  return(p);
}