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


int
bwrmSetHid200SI(/*
      Sets (*m).h_id, and determines (*m).h_adj 
      so that the result of bwrmH(t) for saturated liquid
      at 273.15[K] is 200,000 * (*m).M  [J/kg-mol] */
  BWRM *m,
  double (*f)(double t) /* function that calculates 'h'[J/kg-mole] at
                           ideal gas condition          ~~~~~~~~~~  */
){
  double rho;
  (*m).h_id = f;
  rho = bwrmRHOl(m, bwrmPs(m, T0bwrm), T0bwrm); 
  (*m).h_adj = -(bwrmdH(m, T0bwrm, rho) + (*m).h_id(T0bwrm));
  (*m).h_adj += 200e+3 * (*m).M;
}

double
bwrmH(/*
       Calculate 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 */
  double  t,      /* temperature [K]               */
  double  rho     /* molar density [kg-mol/m3]     */
){
  (*c).h = bwrmdH(c, t, rho) + (*c).h_id(t) + (*c).h_adj;
  return((*c).h );
}


double
bwrmdH(/*
       Calculates the deviation from ideal gas condition on
       molar enthalpy, 'h - h_id'[J/kg-mol]  */
  BWRM   *m,  /* structure of BWR coefficients */
  double  t,      /* temperature [K]               */
  double  rho     /* molar density [kg-mol/m3]     */
){
#define	C   (*m).cf
#define	MT (*m).tmp
  double h, gamma2rho4;

  RecalBWRM(m, t, rho);
  gamma2rho4 =  MT.gammarho2 * MT.gammarho2 ;

  h  = MT.h_T17 * (18.0 - (18.0 + 8.0*MT.gammarho2 - gamma2rho4) * MT.EXP1);
  h += MT.g_T8  * ( 9.0 - ( 9.0 + 3.5*MT.gammarho2 - gamma2rho4) * MT.EXP1);
  h += MT.c_T2  * ( 3.0 - ( 3.0 + 0.5*MT.gammarho2 - gamma2rho4) * MT.EXP1);
  h /= C.gamma;
  h += (1.2*MT.alphaa + 1.4*MT.alphad_T + 2.0*MT.alphae_T4 + 5.8*MT.alphaf_T23)
	    * MT.rho5;
  h += (MT.bRT - 1.5*C.a - 2.0*MT.d_T - 3.5*MT.e_T4 - 13.0*MT.f_T23) * MT.rho2;
  h += (MT.B0RT - 2.0*C.A0 - 4.0*MT.C0_T2 + 5.0*MT.D0_T3 - 6.0*MT.E0_T4) * rho;

  return(h);
#undef MT
#undef C
}


int
bwrmSetSid1SI(/*
      set (*m).s_id, and determines (*m).s_adj 
      so that the result of bwrmH(t) for saturated liquid
      at 273.15[K], 1[Pa] is 1 * (*m).M  [J/kg-mol] */
  BWRM *m,
  double (*f)(double t) /* function that calculates 's'[J/kg-mole] at
                           ideal gas condition          ~~~~~~~~~~  */
){
#define	P0 1.0
  double rho;
  (*m).s_id = f;
  rho = bwrmRHOl(m, bwrmPs(m, T0bwrm), T0bwrm); 
  (*m).s_adj = -( bwrmdS(m, T0bwrm, rho)
                + Rgas * log(P0/(rho*Rgas*T0bwrm)) + (*m).s_id(T0bwrm) 
				 );
  (*m).s_adj += 1.0e+3 * (*m).M;
}

double
bwrmdS(/*
      calculate the deviation from ideal gas condition at p=P0 on 
      molar entropy,  's-s0'[J/kg-molK]. */
  BWRM   *c,    /* structure of BWR coefficients */
  double  t,    /* temperature [K]               */
  double  rho   /* molar density [kg-mol/m3]     */
)
{
#define	C  (*c).cf
#define	MT (*c).tmp
  double 	s;
  RecalBWRM(c, t, rho);

  s  = (2.0*MT.c_T2 + 8.0*MT.g_T8 + 17.0*MT.h_T17) / C.gamma
	   * (1.0 - (1.0 + 0.5 * MT.gammarho2) * exp(-MT.gammarho2));
  s += (MT.d_T + 4.0*MT.e_T4 + 23.0*MT.f_T23) * C.alpha * MT.rho5 / 5.0;
  s += -0.5 * (MT.bRT + MT.d_T + 4.0*MT.e_T4 + 23.0*MT.f_T23) * MT.rho2;
  s += -(MT.B0RT + 2.0*MT.C0_T2 - 3.0*MT.D0_T3 + 4.0*MT.E0_T4) * rho;
  s /= t;

  return(s);
}
double
bwrmS(/*
       Calculates molar entropy  's'[J/kg-molK] using  BWR eq.
       ** This function require coef.s0() */
  BWRM   *c,  /* structure of BWR coefficients */
  double  t,      /* temperature [K]               */
  double  rho     /* molar density [kg-mol/m3]     */
){
#define	P0 1.0
  RecalBWRM(c, t, rho);
  (*c).s  = (*c).s_id(t) +  Rgas * log(P0/MT.rhoRT) + 
               bwrmdS(c, t, rho) + (*c).s_adj;
  return((*c).s );
}