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

/* ------------------------------------------------------------------- 
                             DENSITY
   ------------------------------------------------------------------- */
#define EPS 1e-10   /* limitation of coverfence */
double 
bwrmRHO(/*
        Calculates molar density 'rho' in [kg-mole/m3] 
        using BWR eq. and Newton method */
  BWRM   *c,
  double  p,  /* pressure [Pa]   */
  double  t,  /* temperature [K] */
  double  rhostart
)
{
  int    i;
  double rho, delta;

  (*c).tmp.errorcode = BWRMERR_NOPROBLEM;
  rho = rhostart;
  if (t != (*c).tmp.Tprev)   bwrmCalTmpT(c, t); 

  for(i = 0 ; i < 100 ; i++)	{
    bwrmCalTmpRHO(c, rho);
    bwrmCalTmpOther(c, t, rho);

	/* Newton method ( Xn+1 = Xn - F(Xn)/F'(Xn) ) */
	delta = (bwrmP(c, t, rho) - p) / bwrmdPdrho(c, t, rho);
	rho  -= delta;
	if (rho < 0.0 ) {
	  if (rhostart > (*c).RHOc ) /* liquid */
		rho *= -1.0;
	  else
		rho = (rho+delta)/2.0;   /* gas */
	}
	if (rho > (*c).RHOc && rhostart < (*c).RHOc ) {
		rho = ((rho+delta) + (*c).RHOc)/2.0;   /* gas */
	}
	if (rho < (*c).RHOc && rhostart > (*c).RHOc ) {
		rho = ((rho+delta) + (*c).RHOc)/2.0;   /* liquid */
	}
	  
	if ( fabs(delta / rho) <= EPS ){ /* OK ? */
	  bwrmCalTmpRHO(c, rho);
	  bwrmCalTmpOther(c, t, rho);
	  (*c).tmp.errorcode = BWRMERR_NOPROBLEM;
	  goto RETURN;
	}
  }
  (*c).tmp.errorcode = BWRMERR_CANT_CONVERGE ;
/*   printf("bwrmRHO warning: not convergenced(p=%lfMPa, t=%f[K]\n",
            p/1.0e+6, t); /* Check from (*c).tmp.errorcode! */
 RETURN:
  (*c).rho = rho;
  (*c).tmp.n_iteration = i;
  return(rho);
}


double 
bwrmRHOl(
  BWRM    *c,
  double   p,  /* pressure [Pa]   */
  double   t   /* temperature [K] */
)
{
  double rhostart = 30.0 ;
  (*c).rhol = (*c).rho = bwrmRHO(c, p , t , rhostart ); 
  if ((*c).rhol < (*c).RHOc  ) (*c).tmp.errorcode = BWRMERR_FUNNY_RESULT;
  return( (*c).rhol );
}

double 
bwrmRHOv(
  BWRM *c,
  double    p,  /* pressure [Pa]   */
  double    t   /* temperature [K] */
)
{
  double rhostart = 0.000001 ;
  (*c).rhov = (*c).rho = bwrmRHO(c, p , t , rhostart );
  if ((*c).rhov > (*c).RHOc  ) (*c).tmp.errorcode = BWRMERR_FUNNY_RESULT;
  return( (*c).rhov );
}


#undef EPS



/* ------------------------------------------------------------------- 
                            TEMPERATURE
   ------------------------------------------------------------------- */
#define EPS 1e-7   /* limitation of coverfence */
double
bwrmT(/* Calculates temperature in [K] from P and rho */
  BWRM    *c,
  double   p,     /* pressure [Pa] */
  double   rho,   /* molaxr density [kg-mole/m3] */
  double   tstart /* initial value of temperature for iteration[K] */
)
{
  int 	i;
  double t,delta;
  double bwrmdPdT();  /*  dP/dT          */

  t = tstart;
  if (rho != (*c).tmp.rhoprev)    bwrmCalTmpRHO(c, t);

  (*c).tmp.errorcode = BWRMERR_NOPROBLEM;
  for (i = 0; i < 1000; i++ ){
	/* Newton method [ Xn+1 = Xn - F(Xn)/F'(Xn) ] */
    bwrmCalTmpT(c, t);
    bwrmCalTmpOther(c, t, rho);
	delta  = (bwrmP(c, t , rho) - p) / bwrmdPdT(c, t, rho);
	t -= delta;
	while ( t < 0.0 )
	  {
		delta /= 2.0;
		t += delta;
	  }
	if ( fabs(delta) <= EPS ){ /* OK ? */
	  bwrmCalTmpT(c, t);
	  bwrmCalTmpOther(c, t, rho);
	  goto RETURN;
	}
  }
  (*c).tmp.errorcode = BWRMERR_CANT_CONVERGE ;
  /*  printf("bwrmT warning: not convergenced(p=%lfMPa, rho=%f[kg-mole/m3]\n",
            p/1.0e+6, rho);/**/
 RETURN:
  (*c).tmp.n_iteration = i;
  (*c).T = t;
  return(t);
#undef EPS
}



#define EPS 1e-7   /* limitation of coverfence */

double 
bwrmRHO_L(
  BWRM   *c,
  double  p,  /* pressure [Pa]   */
  double  t,  /* temperature [K] */
  double  rhostart
)
{
  int    i;
  double rho, delta;
  double F, dF, Fpv, rhopv;

  (*c).tmp.errorcode = BWRMERR_NOPROBLEM;
  if (t != (*c).tmp.Tprev)   bwrmCalTmpT(c, t); 

  /* loop to set initial value of rho */
  rho = rhostart - 0.5;
  do{
	rho += 0.5;
	bwrmCalTmpRHO(c, rho); bwrmCalTmpOther(c, t, rho);
	F  = bwrmP(c, t, rho) - p;
  }while(F < 0.0); 
  dF = bwrmdPdrho(c, t, rho);

  /* main loop which calculates rho */
  i = 0;
  while(1){
    Fpv = F; rhopv = rho; 
	delta = F / dF;
	rho  -= delta;
	bwrmCalTmpRHO(c, rho); bwrmCalTmpOther(c, t, rho);
    if (fabs(delta / rho) <  EPS) break; /* -------- OK -------- */
	dF = bwrmdPdrho(c, t, rho);
	F  = bwrmP(c, t, rho) - p;

    /* loop to correct rho so that  F > 0.0 */
	while(F < 0.0){
	  rho = rhopv - (rhopv - rho) * Fpv / (Fpv - F);
	  bwrmCalTmpRHO(c, rho); bwrmCalTmpOther(c, t, rho);
	  F  = bwrmP(c, t, rho) - p;
	}/*---*/

	if ( i++ > 1000){/* --------- can't convergence. ---------- */
	  /* printf("bwrmRHO_L warning: not convergenced(P=%lfMPa, T=%f[C]\n",
			 p/1.0e+6, t-273.15);/* Check from (*c).tmp.errorcode! */
	  (*c).tmp.errorcode = BWRMERR_CANT_CONVERGE ;
	  goto RETURN;
	}
  }
  if (rho < 0.0 ) (*c).tmp.errorcode = BWRMERR_FUNNY_RESULT;

 RETURN:
  (*c).rho = rho;
  (*c).tmp.n_iteration = i;
  return(rho);
}
#ifdef AAA

#endif