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

double
bwrmP(/*
		 Calculate pressure 'P' in [Pa] by BWR eq.        */
  BWRM    *c,  /* structure of BWR cficients */
  double   t,      /* temperature [K]               */
  double   rho     /* molar density [kg-mol/m3]     */
){
#define	C (*c).cf
#define MT (*c).tmp
  double *p;

  p = &((*c).P);
  RecalBWRM(c, t, rho);

  *p =  (MT.c_T2 + MT.g_T8 + MT.h_T17)
		  * MT.rho3 * (1.0 + MT.gammarho2)  * MT.EXP1;
  *p += (MT.alphaa + MT.alphad_T + MT.alphae_T4 + MT.alphaf_T23) * MT.rho6;
  *p += (MT.bRT - C.a - MT.d_T - MT.e_T4 - MT.f_T23) * MT.rho3;
  *p += (MT.B0RT - C.A0 - MT.C0_T2 + MT.D0_T3 - MT.E0_T4) * MT.rho2;
  *p += MT.rhoRT;
  return(*p);
#undef C
#undef MT
}


double 
bwrmdPdrho(/*  calculate  dP/dRHO	*/
  BWRM   *c,
  double  t,      /* temperature [K]               */
  double  rho     /* molar density [kg-mol/m3]     */
)
{
#define	C (*c).cf
#define MT (*c).tmp
  double dp;

  RecalBWRM(c, t, rho);

  dp  = Rgas * t;
  dp += 2.0*( MT.B0RT - C.A0 - MT.C0_T2 + MT.D0_T3 - MT.E0_T4 )*rho;
  dp += 3.0*( MT.bRT - C.a - MT.d_T - MT.e_T4 - MT.f_T23 )*MT.rho2;
  dp += 6.0*( 
         MT.alphaa + MT.alphad_T + MT.alphae_T4 + MT.alphaf_T23 ) * MT.rho5;
  dp += ( MT.c_T2 + MT.g_T8 + MT.h_T17 )*MT.rho2*MT.EXP1
        * ( 3.0 + 3.0 * MT.gammarho2  - 2.0 * MT.gammarho2 * MT.gammarho2);
  return(dp);
#undef C
#undef MT
}


double bwrmdPdT( /* dP/dT */
  BWRM   *c,
  double  t,      /* temperature [K]               */
  double  rho     /* molar density [kg-mol/m3]     */
){
#define	C (*c)
#define	MT (*c).tmp
  double dp ;

  RecalBWRM(c, t, rho);

  dp  =  MT.rhoRT;
  dp += (MT.B0RT + 2.0*MT.C0_T2 - 3.0*MT.D0_T3 + 4.0*MT.E0_T4)*MT.rho2;
  dp += (MT.bRT + MT.d_T - 4.0*MT.e_T4 + 23.0*MT.f_T23)*MT.rho3;
  dp += (-MT.alphad_T -4.0*MT.alphae_T4 - 23.0*MT.alphaf_T23 )*MT.rho6;
  dp += (-2.0*MT.c_T2 - 8.0*MT.g_T8 - 17.0*MT.h_T17)*MT.rho3*(1.0 + MT.gammarho2)*MT.EXP1;
  dp /= t;
  return( dp );
#undef C
#undef	MT
}