WaterPropsIAPWSphi.cpp

Go to the documentation of this file.

/**
 * @file WaterPropsIAPWSphi.cpp
 *  Definitions for Lowest level of the classes which support a real water model
 *  (see class \link Cantera::WaterPropsIAPWS WaterPropsIAPWS\endlink and class  #WaterPropsIAPWSphi).
 */
/*
 * Copywrite (2006) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
/*
 * $Id: WaterPropsIAPWSphi.cpp 311 2009-12-11 01:32:26Z hkmoffa $
 */

#include "WaterPropsIAPWSphi.h"

#include <cstdio>
#include <cmath>

using std::printf;
using std::sqrt;
using std::log;
using std::exp;
using std::pow;
using std::fabs;

/*
 * Critical Point values in mks units: Note, these aren't used in this 
 * routine, except for internal checks. All calculations here are done
 * in dimensionless units.
 */
static const doublereal  T_c = 647.096;  // Kelvin
static const doublereal  P_c = 22.064E6; // Pascals
static const doublereal  Rho_c = 322.;    // kg m-3
static const doublereal  M_water = 18.015268; // kg kmol-1

/*
 * The added constants were calculated so that u = s = 0
 * for liquid at the triple point. These where determined
 * by the program testPress. I'm not quite satisfied with
 * the result, but will let it stand for the moment.
 * H didn't turn out to be .611872 J/kg, but .611782 J/kg.
 * There may be a slight error here somehow.
 */
static const doublereal  ni0[9] = {
  0.0,
  -8.32044648201 - 0.000000001739715,
  6.6832105268   + 0.000000000793232,
  3.00632,
  0.012436,
  0.97315,
  1.27950,
  0.96956,
  0.24873
};

static const doublereal  gammi0[9] = {
  0.0,
  0.0,
  0.0,
  0.0,
  1.28728967,
  3.53734222,
  7.74073708,
  9.24437796,
  27.5075105
};

static const int ciR[56] = {
  0,  //  0
  0,  //  1
  0,
  0,
  0,
  0,  //  5
  0,
  0,
  1,
  1,
  1,  // 10
  1,
  1,
  1,
  1,
  1,  // 15
  1,
  1,
  1,
  1,
  1,  // 20
  1,   
  1,
  2,
  2,
  2,  // 25
  2,
  2,
  2,
  2,
  2,  // 30
  2,
  2,
  2,
  2,
  2,  // 35
  2,
  2,
  2,
  2,
  2,  // 40
  2,
  2,
  3,
  3,
  3,  // 45
  3,
  4,
  6,
  6,
  6,  // 50
  6,
  0,
  0,
  0,
  0   // 55
};

static const int diR[55] = {
  0,  //  0
  1,  //  1
  1,
  1,
  2,
  2,  //  5
  3,
  4,
  1,
  1,
  1,  // 10
  2,
  2,
  3,
  4,
  4,  // 15
  5,
  7,
  9,
  10,
  11, // 20
  13,   
  15,
  1,
  2,
  2,  // 25
  2,
  3,
  4,
  4,
  4,  // 30
  5,
  6,
  6,
  7,
  9,  // 35
  9,
  9,
  9,
  9,
  10, // 40
  10,
  12,
  3,
  4,
  4,  // 45
  5,
  14,
  3,
  6,
  6,  // 50
  6,
  3,
  3,
  3   // 54
};

static const int tiR[55] = {
  0,  //  0
  0,  //  1
  0,
  0,
  0,
  0,  //  5
  0,
  0,
  4,  //  8
  6,
  12, // 10
  1,
  5,
  4,
  2,
  13, // 15
  9,
  3,
  4,
  11,
  4,  // 20
  13,   
  1,
  7,
  1,
  9,  // 25
  10,
  10,
  3,
  7,
  10,  // 30
  10,
  6,
  10,
  10,
  1,  // 35
  2,
  3,
  4,
  8,
  6, // 40
  9,
  8,
  16,
  22,
  23,  // 45
  23,
  10,
  50,
  44,
  46,  // 50
  50,
  0,
  1,
  4   // 54
};

static const doublereal  ni[57] = {
  +0.0,
  +0.12533547935523E-1, //  1
  +0.78957634722828E1,  //  2
  -0.87803203303561E1,  //  3
  +0.31802509345418E0,  //  4
  -0.26145533859358E0,  //  5
  -0.78199751687981E-2, //  6
  +0.88089493102134E-2, //  7
  -0.66856572307965E0,  //  8
  +0.20433810950965,    //  9
  -0.66212605039687E-4, // 10
  -0.19232721156002E0,  // 11
  -0.25709043003438E0,  // 12
  +0.16074868486251E0,  // 13
  -0.40092828925807E-1, // 14
  +0.39343422603254E-6, // 15
  -0.75941377088144E-5, // 16
  +0.56250979351888E-3, // 17
  -0.15608652257135E-4, // 18
  +0.11537996422951E-8, // 19
  +0.36582165144204E-6, // 20
  -0.13251180074668E-11,// 21
  -0.62639586912454E-9, // 22
  -0.10793600908932E0,  // 23
  +0.17611491008752E-1, // 24
  +0.22132295167546E0,  // 25
  -0.40247669763528E0,  // 26
  +0.58083399985759E0,  // 27
  +0.49969146990806E-2, // 28
  -0.31358700712549E-1, // 29
  -0.74315929710341E0,  // 30
  +0.47807329915480E0,  // 31
  +0.20527940895948E-1, // 32
  -0.13636435110343E0,  // 33
  +0.14180634400617E-1, // 34
  +0.83326504880713E-2, // 35
  -0.29052336009585E-1, // 36
  +0.38615085574206E-1, // 37
  -0.20393486513704E-1, // 38
  -0.16554050063734E-2, // 39
  +0.19955571979541E-2, // 40
  +0.15870308324157E-3, // 41
  -0.16388568342530E-4, // 42
  +0.43613615723811E-1, // 43
  +0.34994005463765E-1, // 44
  -0.76788197844621E-1, // 45
  +0.22446277332006E-1, // 46
  -0.62689710414685E-4, // 47
  -0.55711118565645E-9, // 48
  -0.19905718354408E0,  // 49
  +0.31777497330738E0,  // 50
  -0.11841182425981E0,  // 51
  -0.31306260323435E2,  // 52
  +0.31546140237781E2,  // 53
  -0.25213154341695E4,  // 54
  -0.14874640856724E0,  // 55
  +0.31806110878444E0   // 56
};


static const doublereal  alphai[3] = {
  +20.,
  +20.,
  +20.
};

static const doublereal  betai[3] = {
  +150.,
  +150.,
  +250.
};

static const doublereal  gammai[3] = {
  +1.21,
  +1.21,
  +1.25
};

static const doublereal  epsi[3] = {
  +1.0,
  +1.0,
  +1.0
};

static const doublereal  ai[2] = {
  +3.5,
  +3.5
};

static const doublereal  bi[2] = {
  +0.85,
  +0.95
};

static const doublereal  Bi[2] = {
  +0.2,
  +0.2
};

static const doublereal  Ci[2] = {
  +28.0,
  +32.0
};

static const doublereal  Di[2] = {
  +700.,
  +800.
};

static const doublereal  Ai[2] = {
  +0.32,
  +0.32
};

static const doublereal  Bbetai[2] = {
  +0.3,
  +0.3
};

/**
 * Constructor for the object. 
 */
WaterPropsIAPWSphi::WaterPropsIAPWSphi() :
  TAUsave(-1.0),
  TAUsqrt(-1.0),
  DELTAsave(-1.0)
{
  for (int i = 0; i < 52; i++) {
    TAUp[i] = 1.0;
  }
  for (int i = 0; i < 16; i++) {
    DELTAp[i] = 1.0;
  }
}

/*
 * intCheck() calculates all of the functions at a one point and
 * prints out the result. It's used for conducting the internal
 * check.
 */
void WaterPropsIAPWSphi::intCheck(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau    = phi0();
  doublereal  res    = phiR();
  doublereal  res_d  = phiR_d();
  doublereal  nau_d  = phi0_d();
  doublereal  res_dd = phiR_dd();
  doublereal  nau_dd = phi0_dd();
  doublereal  res_t  = phiR_t();
  doublereal  nau_t  = phi0_t();
  doublereal  res_tt = phiR_tt();
  doublereal  nau_tt = phi0_tt();
  doublereal  res_dt = phiR_dt();
  doublereal  nau_dt = phi0_dt();

  std::printf("nau    = %20.12e\t\tres    = %20.12e\n", nau,    res);
  std::printf("nau_d  = %20.12e\t\tres_d  = %20.12e\n", nau_d,  res_d);
  printf("nau_dd = %20.12e\t\tres_dd = %20.12e\n", nau_dd, res_dd);
  printf("nau_t  = %20.12e\t\tres_t  = %20.12e\n", nau_t,  res_t);
  printf("nau_tt = %20.12e\t\tres_tt = %20.12e\n", nau_tt, res_tt);
  printf("nau_dt = %20.12e\t\tres_dt = %20.12e\n", nau_dt, res_dt);
}

void WaterPropsIAPWSphi::check1() {
  doublereal  T = 500.;
  doublereal  rho = 838.025;
  doublereal  tau = T_c/T;
  doublereal  delta = rho / Rho_c;
  printf(" T = 500 K, rho = 838.025 kg m-3\n");
  intCheck(tau, delta);
}

void WaterPropsIAPWSphi::check2() {
  doublereal  T = 647;
  doublereal  rho = 358.0;
  doublereal  tau = T_c/T;
  doublereal  delta = rho / Rho_c;
  printf(" T = 647 K, rho = 358.0 kg m-3\n");
  intCheck(tau, delta);
}

/*
 * Calculate the polynomials in tau and delta, and store them in static 
 * storage.
 */
void WaterPropsIAPWSphi::tdpolycalc(doublereal  tau, doublereal  delta) {
  if ((tau != TAUsave) || 1) {
    TAUsave = tau;
    TAUsqrt = sqrt(tau);
    TAUp[0] = 1.0;
    for (int i = 1; i < 51; i++) {
      TAUp[i] = TAUp[i-1] * tau;
    }
  }
  if ((delta != DELTAsave) || 1) {
    DELTAsave = delta;
    DELTAp[0] = 1.0;
    for (int i = 1; i <= 15; i++) {
      DELTAp[i] = DELTAp[i-1] * delta;
    }
  }
}

/*
 * Calculate Eqn. 6.5 for phi0, the ideal gas part of the
 * dimensionless Helmholtz free energy.
 */
doublereal  WaterPropsIAPWSphi::phi0() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  doublereal  retn = log(delta) + ni0[1] + ni0[2]*tau + ni0[3]*log(tau);

  retn += ni0[4] * log(1.0 - exp(-gammi0[4]*tau));
  retn += ni0[5] * log(1.0 - exp(-gammi0[5]*tau));
  retn += ni0[6] * log(1.0 - exp(-gammi0[6]*tau));
  retn += ni0[7] * log(1.0 - exp(-gammi0[7]*tau));
  retn += ni0[8] * log(1.0 - exp(-gammi0[8]*tau));
  return retn;
}

/*
 * Calculate Eqn. 6.6 for phiR, the residual part of the
 * dimensionless Helmholtz free energy.
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
    
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = (ni[1] * delta / TAUsqrt +
                ni[2] * delta * TAUsqrt * T375 +
                ni[3] * delta * tau +
                ni[4] * DELTAp[2] * TAUsqrt +
                ni[5] * DELTAp[2] * T375 * T375 +
                ni[6] * DELTAp[3] * T375 +
                ni[7] * DELTAp[4] * tau);
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * exp(-DELTAp[ciR[i]]));
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    val += (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * 
            exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    doublereal  atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
      
    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
    val += (ni[i] * triagtmp * delta * phi);
  }

  return val;
}

/*
 * Calculate the Phi function, which is basically the helmholtz free energy
 * Eqn. (6.4)
 */
doublereal  WaterPropsIAPWSphi::phi(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau = phi0();
  doublereal  res = phiR();
  doublereal  retn = nau + res;
  return retn;
}


/*
 * Calculate d_phiR_d(delta), the first derivative of phiR
 * wrt delta
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR_d() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
    
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = (ni[1] / TAUsqrt +
                ni[2] * TAUsqrt * T375 +
                ni[3] * tau +
                ni[4] * 2.0 * delta * TAUsqrt +
                ni[5] * 2.0 * delta * T375 * T375 +
                ni[6] * 3.0 * DELTAp[2] * T375 +
                ni[7] * 4.0 * DELTAp[3] * tau);
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    val += ( (ni[i] * exp(-DELTAp[ciR[i]]) *  DELTAp[diR[i] - 1] *
              TAUp[tiR[i]]) * (diR[i] - ciR[i]* DELTAp[ciR[i]]));
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    doublereal  tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * 
                  exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
    val += tmp * (diR[i]/delta - 2.0 * alphai[j] * dtmp);
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    doublereal  atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
    doublereal  triagtmpm1 = pow(triag, bi[j]-1.0);
    doublereal  atmpM1 = atmp - 1.0;
    doublereal  ptmp = pow(dtmp2,atmpM1); 
    doublereal  p2tmp = pow(dtmp2, ai[j]-1.0);
    doublereal  dtriagddelta = 
      deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
                2.0*Bi[j]*ai[j]*p2tmp);

    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
    doublereal  dphiddelta = -2.0*Ci[j]*deltam1*phi;
    doublereal  dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;

    doublereal  tmp = ni[i] * (triagtmp * (phi + delta*dphiddelta) +
                          dtriagtmpddelta * delta * phi);
    val += tmp;
  }

  return val;
}

/*
 * Calculate d_phi0_d(delta), the first derivative of phi0
 * wrt delta
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phi0_d() const {
  doublereal  delta = DELTAsave;
  return (1.0/delta);
}

/*
 * Calculate the dPhidDelta function, which is basically the derivative
 * of helmholtz free energy wrt delta
 * Eqn. (6.4)
 */
doublereal  WaterPropsIAPWSphi::phi_d(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau = phi0_d();
  doublereal  res = phiR_d();
  doublereal  retn = nau + res;
  return retn;
}

/*
 * Calculate the dimensionless pressure at tau and delta;
 *
 *       p/(rhoRT) = delta * phi_d()
 *
 * note: this is done so much, we have a seperate routine.
 */
doublereal  WaterPropsIAPWSphi::pressureM_rhoRT(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  res = phiR_d();
  doublereal  retn = 1.0 + delta * res; 
  return retn;
}

/*
 * Calculate d2_phiR_dd(delta), the second derivative of phiR
 * wrt delta
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR_dd() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
  doublereal  atmp;
    
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = (ni[4] * 2.0 * TAUsqrt +
                ni[5] * 2.0 * T375 * T375 +
                ni[6] * 6.0 * delta * T375 +
                ni[7] * 12.0 * DELTAp[2] * tau);
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    doublereal  dtmp = DELTAp[ciR[i]];
    doublereal  tmp = ni[i] * exp(-dtmp) * TAUp[tiR[i]];
    if (diR[i] == 1) {
      atmp = 1.0/delta;
    } else {
      atmp = DELTAp[diR[i] - 2];
    }
    tmp *= atmp *((diR[i] - ciR[i]*dtmp)*(diR[i]-1.0-ciR[i]*dtmp) -
                  ciR[i]*ciR[i]*dtmp);
    val += tmp;
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    doublereal  tmp = (ni[i] * TAUp[tiR[i]] * 
                  exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
    doublereal  deltmp =  DELTAp[diR[i]];
    doublereal  deltmpM1 = deltmp/delta;
    doublereal  deltmpM2 = deltmpM1 / delta;
    doublereal  d2tmp = dtmp * dtmp;
      
    val += tmp * (-2.0*alphai[j]*deltmp +
                  4.0 * alphai[j] * alphai[j] * deltmp * d2tmp -
                  4.0 * diR[i] * alphai[j] * deltmpM1 * dtmp +
                  diR[i] * (diR[i] - 1.0) * deltmpM2);
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
    doublereal  triagtmpm1 = pow(triag, bi[j]-1.0);
    doublereal  atmpM1 = atmp - 1.0;
    doublereal  ptmp = pow(dtmp2,atmpM1);
    doublereal  p2tmp = pow(dtmp2, ai[j]-1.0);
    doublereal  dtriagddelta = 
      deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
                2.0*Bi[j]*ai[j]*p2tmp);

    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
    doublereal  dphiddelta = -2.0*Ci[j]*deltam1*phi;
    doublereal  dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;


    doublereal  d2phiddelta2 = 2.0 * Ci[j] * phi * (2.0*Ci[j]*dtmp2 - 1.0);

    doublereal  pptmp = ptmp / dtmp2;
    doublereal  d2triagddelta2 = dtriagddelta / deltam1;
    d2triagddelta2 += 
      dtmp2 *(4.0*Bi[j]*ai[j]*(ai[j]-1.0)*pow(dtmp2,ai[j]-2.0) +
              2.0*Ai[j]*Ai[j]/(Bbetai[j]*Bbetai[j])*ptmp*ptmp +
              Ai[j]*theta*4.0/Bbetai[j]*(atmp-1.0)*pptmp);

    doublereal   d2triagtmpd2delta = 
      bi[j] * (triagtmpm1 * d2triagddelta2 +
               (bi[j]-1.0)*triagtmpm1/triag*dtriagddelta*dtriagddelta);

    doublereal  ctmp = (triagtmp * (2.0*dphiddelta + delta*d2phiddelta2) +
                   2.0*dtriagtmpddelta*(phi + delta * dphiddelta) +
                   d2triagtmpd2delta * delta * phi);
 
    val += ni[i] * ctmp;
  }

  return val;
}

/*
 * Calculate d2_phi0_dd(delta), the second derivative of phi0
 * wrt delta
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phi0_dd() const {
  doublereal  delta = DELTAsave;
  return (-1.0/(delta*delta));
}

/*
 * Calculate the d2_PhidDelta2 function, which is the second derivative
 * of helmholtz free energy wrt delta
 * Eqn. (6.4)
 */
doublereal  WaterPropsIAPWSphi::phi_dd(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau = phi0_dd();
  doublereal  res = phiR_dd();
  doublereal  retn = nau + res;
  return retn;
}

doublereal  WaterPropsIAPWSphi::dimdpdrho(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  res1 = phiR_d();
  doublereal  res2 = phiR_dd();
  doublereal  retn = 1.0 + delta * (2.0*res1 + delta*res2);
  return retn;
}

doublereal  WaterPropsIAPWSphi::dimdpdT(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  res1 = phiR_d();
  doublereal  res2 = phiR_dt();
  doublereal  retn = (1.0 + delta * res1) - tau * delta * (res2);
  return retn;
}

/*
 * Calculate d_phi0/d(tau)
 */
doublereal  WaterPropsIAPWSphi::phi0_t() const {
  doublereal  tau = TAUsave;
  doublereal  retn = ni0[2] + ni0[3]/tau;;
  retn += (ni0[4] * gammi0[4] * (1.0/(1.0 - exp(-gammi0[4]*tau)) - 1.0));
  retn += (ni0[5] * gammi0[5] * (1.0/(1.0 - exp(-gammi0[5]*tau)) - 1.0));
  retn += (ni0[6] * gammi0[6] * (1.0/(1.0 - exp(-gammi0[6]*tau)) - 1.0));
  retn += (ni0[7] * gammi0[7] * (1.0/(1.0 - exp(-gammi0[7]*tau)) - 1.0));
  retn += (ni0[8] * gammi0[8] * (1.0/(1.0 - exp(-gammi0[8]*tau)) - 1.0));
  return retn;
}

/*
 * Calculate Eqn. 6.6 for dphiRdtau, the derivative residual part of the
 * dimensionless Helmholtz free energy wrt temperature
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR_t() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
  doublereal  atmp, tmp;
    
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = ((-0.5) *ni[1] * delta / TAUsqrt / tau +
                ni[2] * delta * 0.875 / TAUsqrt * T375 +
                ni[3] * delta +
                ni[4] * DELTAp[2] * 0.5 / TAUsqrt +
                ni[5] * DELTAp[2] * 0.75 * T375 * T375 / tau +
                ni[6] * DELTAp[3] * 0.375 * T375 / tau +
                ni[7] * DELTAp[4]);
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-1] * exp(-DELTAp[ciR[i]]));
    val += tiR[i] * tmp;
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * 
           exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
    val += tmp *(tiR[i]/tau - 2.0 * betai[j]*ttmp);
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
      
    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
 

    doublereal  dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;

    doublereal  dphidtau = - 2.0 * Di[j] * ttmp * phi;
      
    val += ni[i] * delta * (dtriagtmpdtau * phi + triagtmp * dphidtau);
  }

  return val;
}

/*
 * Calculate the dPhidtau function, which is basically the derivative
 * of helmholtz free energy wrt tau
 * Eqn. (6.4)
 */
doublereal  WaterPropsIAPWSphi::phi_t(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau = phi0_t();
  doublereal  res = phiR_t();
  doublereal  retn = nau + res;
  return retn;
}

/*
 * Calculate d2_phi0/dtau2
 */
doublereal  WaterPropsIAPWSphi::phi0_tt() const {
  doublereal  tau = TAUsave;
  doublereal  tmp, itmp;
  doublereal  retn = - ni0[3]/(tau * tau);
  for (int i = 4; i <= 8; i++) {
    tmp = exp(-gammi0[i]*tau);
    itmp = 1.0 - tmp;
    retn -= (ni0[i] * gammi0[i] * gammi0[i] * tmp / (itmp * itmp));
  }
  return retn;
}

/*
 * Calculate Eqn. 6.6 for dphiRdtau, the second derivative residual part of the
 * dimensionless Helmholtz free energy wrt temperature
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR_tt() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
  doublereal  atmp, tmp;
    
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = ((-0.5) * (-1.5) * ni[1] * delta / (TAUsqrt * tau * tau) +
                ni[2] * delta * 0.875 * (-0.125) * T375 / (TAUsqrt * tau) +
                ni[4] * DELTAp[2] * 0.5 * (-0.5)/ (TAUsqrt * tau) +
                ni[5] * DELTAp[2] * 0.75 *(-0.25) * T375 * T375 / (tau * tau) +
                ni[6] * DELTAp[3] * 0.375 *(-0.625) * T375 / (tau * tau));
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    if (tiR[i] > 1) {
      tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]-2] * exp(-DELTAp[ciR[i]]));
      val += tiR[i] *  (tiR[i] - 1.0) * tmp;
    }
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * 
           exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
    atmp = tiR[i]/tau - 2.0 * betai[j]*ttmp;
    val += tmp *(atmp * atmp - tiR[i]/(tau*tau) - 2.0*betai[j]);
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
    doublereal  triagtmpM1 = triagtmp / triag;
      
    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
 

    doublereal  dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;

    doublereal  dphidtau = - 2.0 * Di[j] * ttmp * phi;

    doublereal  d2triagtmpdtau2 = 
      (2 * bi[j] * triagtmpM1 +
       4 * theta * theta  * bi[j] * (bi[j]-1.0) * triagtmpM1 / triag);

    doublereal  d2phidtau2 = 2.0*Di[j]*phi *(2.0*Di[j]*ttmp*ttmp - 1.0);

    tmp =  (d2triagtmpdtau2 * phi +
            2 * dtriagtmpdtau * dphidtau +
            triagtmp * d2phidtau2);
    val += ni[i] * delta * tmp;
  }

  return val;
}

/*
 * Calculate the d2Phidtau2 function, which is basically the second derivative
 * of helmholtz free energy wrt tau
 * Eqn. (6.4)
 */
doublereal  WaterPropsIAPWSphi::phi_tt(doublereal  tau, doublereal  delta) {
  tdpolycalc(tau, delta);
  doublereal  nau = phi0_tt();
  doublereal  res = phiR_tt();
  doublereal  retn = nau + res;
  return retn;
}

/**
 * Calculate d2_phi0/dtauddelta
 */
doublereal  WaterPropsIAPWSphi::phi0_dt() const {
  return 0.0;
}

/*
 * Calculate d2_phiR_d(delta)d(tau), the mixed derivative of phi
 * wrt delta and tau.
 *
 *  tau = dimensionless temperature
 *  delta = dimensionless pressure
 */
doublereal  WaterPropsIAPWSphi::phiR_dt() const {
  doublereal  tau = TAUsave;
  doublereal  delta = DELTAsave;
  int i, j;
  doublereal  tmp;
  /*
   * Write out the first seven polynomials in the expression
   */
  doublereal  T375 = pow(tau, 0.375);
  doublereal  val = (ni[1] * (-0.5) / (TAUsqrt * tau) +
                ni[2] * (0.875) * T375 / TAUsqrt +
                ni[3] +
                ni[4] * 2.0 * delta * (0.5) / TAUsqrt +
                ni[5] * 2.0 * delta * (0.75) * T375 * T375 / tau +
                ni[6] * 3.0 * DELTAp[2] * 0.375 * T375 / tau +
                ni[7] * 4.0 * DELTAp[3]);
  /*
   * Next, do polynomial contributions 8 to 51
   */
  for (i = 8; i <= 51; i++) {
    tmp = (ni[i] * tiR[i] * exp(-DELTAp[ciR[i]]) *  DELTAp[diR[i] - 1] *
           TAUp[tiR[i] - 1]);
    val += tmp * (diR[i] - ciR[i] * DELTAp[ciR[i]]);
  }

  /*
   * Next do contributions 52 to 54
   */
  for (j = 0; j < 3; j++) {
    i = 52 + j;
    doublereal  dtmp = delta - epsi[j];
    doublereal  ttmp = tau - gammai[j];
    tmp = (ni[i] * DELTAp[diR[i]] * TAUp[tiR[i]] * 
           exp(-alphai[j]*dtmp*dtmp - betai[j]*ttmp*ttmp));
    val += tmp * ((diR[i]/delta - 2.0 * alphai[j] * dtmp) * 
                  (tiR[i]/tau   - 2.0 * betai[j]  * ttmp));
  }

  /*
   * Next do contributions 55 and 56
   */
  for (j = 0; j < 2; j++) {
    i = 55 + j;
    doublereal  deltam1 = delta - 1.0;
    doublereal  dtmp2 = deltam1 * deltam1;
    doublereal  atmp = 0.5 / Bbetai[j];
    doublereal  theta = (1.0 - tau) + Ai[j] * pow(dtmp2, atmp);
    doublereal  triag = theta * theta + Bi[j] * pow(dtmp2, ai[j]);
    doublereal  ttmp = tau - 1.0;
   
    doublereal  triagtmp = pow(triag, bi[j]);
    doublereal  triagtmpm1 = pow(triag, bi[j]-1.0);
    doublereal  atmpM1 = atmp - 1.0;
    doublereal  ptmp = pow(dtmp2,atmpM1); 
    doublereal  p2tmp = pow(dtmp2, ai[j]-1.0);
    doublereal  dtriagddelta = 
      deltam1 *(Ai[j] * theta * 2.0 / Bbetai[j] * ptmp +
                2.0*Bi[j]*ai[j]*p2tmp);

    doublereal  phi = exp(-Ci[j]*dtmp2 - Di[j]*ttmp*ttmp);
    doublereal  dphiddelta = -2.0*Ci[j]*deltam1*phi;
    doublereal  dtriagtmpddelta = bi[j] * triagtmpm1 * dtriagddelta;


    doublereal  dtriagtmpdtau = -2.0*theta * bi[j] * triagtmp / triag;

    doublereal  dphidtau = - 2.0 * Di[j] * ttmp * phi;

    doublereal  d2phiddeltadtau = 4.0 * Ci[j] * Di[j] * deltam1 * ttmp * phi;

    doublereal   d2triagtmpddeltadtau =
      ( -Ai[j] * bi[j] * 2.0 / Bbetai[j] * triagtmpm1 * deltam1 * ptmp
        -2.0 * theta * bi[j] * (bi[j] - 1.0) * triagtmpm1 / triag * dtriagddelta);
              

    doublereal  tmp = ni[i] * (triagtmp * (dphidtau + delta*d2phiddeltadtau) +
                          delta * dtriagtmpddelta * dphidtau +
                          dtriagtmpdtau * (phi + delta * dphiddelta) +
                          d2triagtmpddeltadtau * delta * phi);
    val += tmp;
  }

  return val;
}

/**
 * This program computes the reduced density, given the reduced pressure
 * and the reduced temperature, tau. It takes an initial guess, deltaGuess.
 * DeltaGuess is important as this is a multivalued function below the
 * critical point.
 * 
 */
doublereal  WaterPropsIAPWSphi::dfind(doublereal  p_red, doublereal  tau, doublereal  deltaGuess) {
  doublereal  dd = deltaGuess;
  bool conv = false;
  doublereal  deldd = dd;
  doublereal  pcheck = 1.0E-30 + 1.0E-8 * p_red;
  for (int n = 0; n < 200; n++) {
    /*
     * Calculate the internal polynomials, and then calculate the
     * phi deriv functions needed by this routine.
     */
    tdpolycalc(tau, dd);
    doublereal  q1 = phiR_d();
    doublereal  q2 = phiR_dd();

    /*
     * Calculate the predicted reduced pressure, pred0, based on the
     * current tau and dd. 
     */
    doublereal  pred0 = dd + dd * dd * q1;
    /*
     * Calculate the derivative of the predicted reduced pressure
     * wrt the reduced density, dd, This is dpddelta
     */
    doublereal  dpddelta = 1.0 + 2.0 * dd * q1 + dd * dd * q2; 
    /*
     * If dpddelta is negative, then we are in the middle of the
     * 2 phase region, beyond the stability curve. We need to adjust
     * the initial guess outwards and start a new iteration.
     */
    if (dpddelta <= 0.0) {
      if (deltaGuess > 1.0) dd = dd * 1.05;
      if (deltaGuess < 1.0) dd = dd * 0.95;
      continue;
    }
    /*
     * Check for convergence
     */
    if (fabs(pred0-p_red) < pcheck) {
      conv = true;
      break;
    }

    /*
     * Dampen and crop the update
     */
    doublereal  dpdx = dpddelta;
    if (n < 10) {
      dpdx = dpddelta * 1.1;
    }
    if (dpdx < 0.001) dpdx = 0.001;
      
    /*
     * Formulate the update to reduced density using 
     * Newton's method. Then, crop it to a max value
     * of 0.02
     */
    deldd = - (pred0 - p_red) / dpdx;
    if (fabs(deldd) > 0.05) {
      deldd = deldd * 0.05 / fabs(deldd);
    }
    /*
     * updated the reduced density value
     */
    dd = dd + deldd;
    if (fabs(deldd/dd) < 1.0E-14) {
      conv = true;
      break;
    }
    /*
     * Check for negative densities
     */
    if (dd <= 0.0) {
      dd = 1.0E-24;
    }
  }
  /*
   * Check for convergence, and return 0.0 if it wasn't achieved.
   */
  if (! conv) {
    dd = 0.0;
  }
  return dd;
}

/**
 * Calculate the dimensionless gibbs free energy g/RT.
 */
doublereal  WaterPropsIAPWSphi::gibbs_RT() const {
  doublereal  delta = DELTAsave;
  doublereal  rd = phiR_d();
  doublereal  g = 1.0 + phi0() + phiR() + delta * rd;
  return g;
}

/**
 * Calculate the dimensionless enthalpy h/RT.
 */
doublereal  WaterPropsIAPWSphi::enthalpy_RT() const {
  doublereal  delta = DELTAsave;
  doublereal  tau   = TAUsave;
  doublereal  rd = phiR_d();
  doublereal  nt = phi0_t();
  doublereal  rt = phiR_t();
  doublereal  hRT = 1.0 + tau * (nt + rt) + delta * rd;
  return hRT;
}

/*
 * Calculate the dimensionless entropy s/R.
 */
doublereal  WaterPropsIAPWSphi::entropy_R() const {
  doublereal  tau   = TAUsave;
  doublereal  nt = phi0_t();
  doublereal  rt = phiR_t();
  doublereal  p0 = phi0();
  doublereal  pR = phiR();
  doublereal  sR = tau * (nt + rt) - p0 - pR;
  return sR;
}

/*
 * Calculate the dimensionless internal energy, u/RT.
 */
doublereal  WaterPropsIAPWSphi::intEnergy_RT() const {
  doublereal  tau   = TAUsave;
  doublereal  nt = phi0_t();
  doublereal  rt = phiR_t();
  doublereal  uR = tau * (nt + rt);
  return uR;
}

/*
 * Calculate the dimensionless constant volume Heat Capacity, Cv/R
 */
doublereal  WaterPropsIAPWSphi::cv_R() const {
  doublereal  tau   = TAUsave;
  doublereal  ntt = phi0_tt();
  doublereal  rtt = phiR_tt();
  doublereal  cvR = - tau * tau * (ntt + rtt);
  return cvR;
}

/*
 * Calculate the dimensionless constant pressure Heat Capacity, Cp/R
 */
doublereal  WaterPropsIAPWSphi::cp_R() const {
  doublereal  tau   = TAUsave;
  doublereal  delta = DELTAsave;
  doublereal  cvR = cv_R();
  //doublereal  nd = phi0_d();
  doublereal  rd = phiR_d();
  doublereal  rdd = phiR_dd();
  doublereal  rdt = phiR_dt();
  doublereal  num = (1.0 + delta * rd - delta * tau * rdt);
  doublereal  cpR = cvR + (num * num / 
                      (1.0 + 2.0 * delta * rd + delta * delta * rdd));
  return cpR;
}