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

double
bwrmCp(/*
       Calculate molar specific heat with constant pressure,
       'cp' [J/kg-mol K].
       ** This function require coef.cp0() **  */
  BWRM    *c,   /* structure of BWR coefficients */
  double   t,   /* temperature [K]               */
  double   rho  /* molar density [kg-mol/m3]     */
){
  (*c).cp = bwrmdCp(c, t, rho) + (*c).cp_id(t);
  return((*c).cp);
}
double
bwrmdCp(/*
       Calculates the deviation from the ideal gas condition on 
       molar specific heat with constant pressure. [J/kg-mol K] */
  BWRM    *coef,
  double   t,      /* temperature [K]               */
  double   rho     /* molar density [kg-mol/m3]     */
)
{
  double  dCp;
  double  dPdT2();
  int mode = 0;
  double bwrmdPdT(),bwrmdPdrho();/* in bwrm3.c (dP/dRHO) */
  double bwrmdHdT();
  double bwrmdHdRHO();

  RecalBWRM(coef, t, rho);
  dCp  = bwrmdHdT(coef, t , rho  );
  dCp -= bwrmdHdRHO(coef, t , rho) * bwrmdPdT(coef, t , rho)
             /  bwrmdPdrho(coef ,t , rho);
  return( dCp );
}


	
double bwrmdHdT( /* dh/dT */
  BWRM    *coef,
  double   t,      /* temperature [K]               */
  double   rho     /* molar density [kg-mol/m3]     */
)
{
#define	MT (*coef).tmp
#define	C (*coef).cf
  double dh , g2d4 ;
  g2d4 = MT.gammarho2 * MT.gammarho2;
  
  dh  = (MT.B0RT + 8.0*MT.C0_T2 - 15.0*MT.D0_T3 + 24.0*MT.E0_T4)*rho;
  dh += (MT.bRT + 2.0*MT.d_T + 14.0*MT.e_T4 + 299.0*MT.f_T23 )*MT.rho2;
  dh += (-1.4*MT.alphad_T - 8.0*MT.alphae_T4 - 133.4*MT.alphaf_T23 )*MT.rho5;
  dh += -2.0*MT.c_T2/C.gamma*(3.0 - (3.0 + MT.gammarho2/2.0 - g2d4)*MT.EXP1);
  dh += -8.0*MT.g_T8/C.gamma*(9.0 - (9.0 + 3.5*MT.gammarho2 - g2d4)*MT.EXP1);
  dh += -17.0*MT.h_T17/C.gamma*(18.0 - (18.0 + 8.0*MT.gammarho2 - g2d4)*MT.EXP1);
  dh /= t;
  
  return( dh );
#undef	MT
#undef	C
}



double bwrmdHdRHO( /* dh/drho */
  BWRM    *coef,
  double   t,  /* temperature [K]            */
  double   rho /* molar density [kg-mol/m3]  */
)
{
#define	MT (*coef).tmp
#define	C (*coef).cf
  double dh , g2d4 ;

  g2d4 = MT.gammarho2 * MT.gammarho2;

  dh  = (MT.B0RT - 2.0*C.A0 - 4.0*MT.C0_T2 + 5.0*MT.D0_T3 - 6.0*MT.E0_T4);
  dh += (MT.bRT - 1.5*C.a - 2.0*MT.d_T - 3.5*MT.e_T4 - 13.0*MT.f_T23 )*2.0*rho;
  dh += (6.0*MT.alphaa + 7.0*MT.alphad_T + 10.0*MT.alphae_T4 + 29.0*MT.alphaf_T23 )
          * MT.rho4;
  dh += 2.0*MT.c_T2*(2.5 + 2.5*MT.gammarho2 - g2d4)*MT.EXP1*rho;
  dh += 2.0*MT.g_T8*(5.5 + 5.5*MT.gammarho2 - g2d4)*MT.EXP1*rho;
  dh += 2.0*MT.h_T17*(10.0 + 10.0*MT.gammarho2 - g2d4)*MT.EXP1*rho;
  
  return( dh );
#undef	MT
#undef	C
}