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

double
bwrmT_Ph(/* 
      Determines temperature from pressure,P[Pa] and 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.  */
  BWRM   *c,    /* 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]     */
#define TMP     (*c).tmp
#define EPS     1e-4
{
#define LIQUID  (mode == 'l' || mode == 'L') 
#define VAPOR   (mode == 'v' || mode == 'V') 
  double t, dt;
  double rho;
  int    n = 0;
  double D = 0.001;

  if( LIQUID ){
	rho = 30.0;
	t = T0bwrm;
  }
  else if( VAPOR ){
	rho = 0.00001;
	t = (*c).Tc;
  }else{
	/* fprintf(stderr, "bwrmT_Ph FATAL ERROR: mode1 must be 'l' or 'v'\n");
	fprintf(stderr, "                    : Check your program!\n");
	exit(1);/**/
	TMP.errorcode = BWRMERR_CANT_CONVERGE ;
	return -1.0;
  }
  do{
	rho = bwrmRHO(c, p, t, rho);
	dt = (bwrmH(c, t+D, rho) - bwrmH(c, t-D, rho)) / (2.0*D);
	if (dt == 0.0) dt = 1e+10;
	dt = (bwrmH(c, t, rho) - h) / dt;
	t -= dt;
	if (n++ > 1000){
	  /* printf("bwrmT_Ph warning: not convergenced(h=%lf[kJ/kg])\n",
			 h/1.0e+3 / (*c).M); /**/
	  TMP.errorcode = BWRMERR_CANT_CONVERGE ;
	  goto RETURN;
	}
  }while(fabs(dt) > EPS );
  TMP.errorcode = BWRMERR_NOPROBLEM;
  goto RETURN;
  
 RETURN:
  (*c).rho = rho;
  if( LIQUID ) (*c).rhol = rho;
  else         (*c).rhov = rho;
  (*c).T = t;
  TMP.n_iteration = n;
  return((*c).T);
}