{

private:

std::map<backend_families, shared_ptr<AbstractStateGenerator>> backends;

public:

void add_backend(const backend_families& bg, const shared_ptr<AbstractStateGenerator>& asg) {

backends[bg] = asg;

};

void get_generator_iterators(const backend_families& bg,

std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& generator,

std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& end) {

generator = backends.find(bg);

end = backends.end();

};

std::size_t size() {

return backends.size();

};

{

public:

AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

if (fluid_names.size() == 1) { // Check that fluid_names[0] has only one component

std::string str = fluid_names[0]; // Check that the fluid name is an alias for "Water"

if ((upper(str) == "WATER") || (upper(str) == "H2O")) {

return new IF97Backend();

} else {

throw ValueError(format("The IF97 backend returns Water props only; fluid name [%s] not allowed", fluid_names[0].c_str()));

}

} else {

throw ValueError(format("The IF97 backend does not support mixtures, only Water"));

};

};

{

public:

AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

return new SRKBackend(fluid_names, get_config_double(R_U_CODATA));

};

{

public:

AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA));

};

{

public:

AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

if (fluid_names.size() != 1) {

throw ValueError(format("For INCOMP backend, name vector must be one element long"));

}

return new IncompressibleBackend(fluid_names[0]);

};

{

public:

CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

return new CoolProp::VTPRBackend(fluid_names, CoolProp::get_config_double(R_U_CODATA));

};

{

public:

CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {

return new CoolProp::PCSAFTBackend(fluid_names);

};

// Reset all instances of CachedElement and overwrite

// the internal double values with -_HUGE

this->_R = _HUGE;

this->_gas_constant.clear();

this->_molar_mass.clear();

this->_critical.fill(_HUGE);

this->_reducing.fill(_HUGE);

return true;

// Reset all instances of CachedElement and overwrite

// the internal double values with -_HUGE

this->_R = _HUGE;

this->_gas_constant.clear();

this->_molar_mass.clear();

/// Ancillary curve values

this->_rhoLanc.clear();

this->_rhoVanc.clear();

this->_pVanc.clear();

this->_pLanc.clear();

this->_TVanc.clear();

this->_TLanc.clear();

this->_critical.fill(_HUGE);

this->_reducing.fill(_HUGE);

/// Bulk values

this->_rhomolar = -_HUGE;

this->_T = -_HUGE;

this->_p = -_HUGE;

this->_Q = -_HUGE;

this->_tau.clear();

this->_delta.clear();

this->_umolar.clear();

this->_cpmolar.clear();

this->_cp0molar.clear();

this->_cvmolar.clear();

this->_speed_sound.clear();

this->_hmolar.clear();

this->_smolar.clear();

this->_gibbsmolar.clear();

this->_helmholtzmolar.clear();

this->_logp.clear();

this->_logrhomolar.clear();

this->_hmolar_excess.clear();

this->_smolar_excess.clear();

this->_gibbsmolar_excess.clear();

this->_volumemolar_excess.clear();

this->_umolar_excess.clear();

this->_helmholtzmolar_excess.clear();

this->_hmolar_residual.clear();

this->_smolar_residual.clear();

this->_gibbsmolar_residual.clear();

/// Smoothing values

this->_rho_spline.clear();

this->_drho_spline_dh__constp.clear();

this->_drho_spline_dp__consth.clear();

/// Cached low-level elements for in-place calculation of other properties

this->_alpha0.clear();

this->_dalpha0_dTau.clear();

this->_dalpha0_dDelta.clear();

this->_d2alpha0_dTau2.clear();

this->_d2alpha0_dDelta_dTau.clear();

this->_d2alpha0_dDelta2.clear();

this->_d3alpha0_dTau3.clear();

this->_d3alpha0_dDelta_dTau2.clear();

this->_d3alpha0_dDelta2_dTau.clear();

this->_d3alpha0_dDelta3.clear();

this->_alphar.clear();

this->_dalphar_dTau.clear();

this->_dalphar_dDelta.clear();

this->_d2alphar_dTau2.clear();

this->_d2alphar_dDelta_dTau.clear();

this->_d2alphar_dDelta2.clear();

this->_d3alphar_dTau3.clear();

this->_d3alphar_dDelta_dTau2.clear();

this->_d3alphar_dDelta2_dTau.clear();

this->_d3alphar_dDelta3.clear();

this->_dalphar_dDelta_lim.clear();

this->_d2alphar_dDelta2_lim.clear();

this->_d2alphar_dDelta_dTau_lim.clear();

this->_d3alphar_dDelta2_dTau_lim.clear();

/// Two-Phase variables

this->_rhoLmolar.clear();

this->_rhoVmolar.clear();

/// Transport properties

this->_viscosity.clear();

this->_conductivity.clear();

this->_surface_tension.clear();

return true;

// Check if a mass based input, convert it to molar units

switch (input_pair) {

case DmassT_INPUTS: ///< Mass density in kg/m^3, Temperature in K

//case HmassT_INPUTS: ///< Enthalpy in J/kg, Temperature in K (NOT CURRENTLY IMPLEMENTED)

case SmassT_INPUTS: ///< Entropy in J/kg/K, Temperature in K

//case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)

case DmassP_INPUTS: ///< Mass density in kg/m^3, Pressure in Pa

case DmassQ_INPUTS: ///< Mass density in kg/m^3, molar quality

case HmassP_INPUTS: ///< Enthalpy in J/kg, Pressure in Pa

case PSmass_INPUTS: ///< Pressure in Pa, Entropy in J/kg/K

case PUmass_INPUTS: ///< Pressure in Pa, Internal energy in J/kg

case HmassSmass_INPUTS: ///< Enthalpy in J/kg, Entropy in J/kg/K

case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg

case DmassHmass_INPUTS: ///< Mass density in kg/m^3, Enthalpy in J/kg

case DmassSmass_INPUTS: ///< Mass density in kg/m^3, Entropy in J/kg/K

case DmassUmass_INPUTS: ///< Mass density in kg/m^3, Internal energy in J/kg

{

// Set the cache value for the molar mass if it hasn't been set yet

molar_mass();

// Molar mass (just for compactness of the following switch)

CoolPropDbl mm = static_cast<CoolPropDbl>(_molar_mass);

switch (input_pair) {

case DmassT_INPUTS:

input_pair = DmolarT_INPUTS;

value1 /= mm;

break;

//case HmassT_INPUTS: input_pair = HmolarT_INPUTS; value1 *= mm; break; (NOT CURRENTLY IMPLEMENTED)

case SmassT_INPUTS:

input_pair = SmolarT_INPUTS;

value1 *= mm;

break;

//case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm; break; (NOT CURRENTLY IMPLEMENTED)

case DmassP_INPUTS:

input_pair = DmolarP_INPUTS;

value1 /= mm;

break;

case DmassQ_INPUTS:

input_pair = DmolarQ_INPUTS;

value1 /= mm;

break;

case HmassP_INPUTS:

input_pair = HmolarP_INPUTS;

value1 *= mm;

break;

case PSmass_INPUTS:

input_pair = PSmolar_INPUTS;

value2 *= mm;

break;

case PUmass_INPUTS:

input_pair = PUmolar_INPUTS;

value2 *= mm;

break;

case HmassSmass_INPUTS:

input_pair = HmolarSmolar_INPUTS;

value1 *= mm;

value2 *= mm;

break;

case SmassUmass_INPUTS:

input_pair = SmolarUmolar_INPUTS;

value1 *= mm;

value2 *= mm;

break;

case DmassHmass_INPUTS:

input_pair = DmolarHmolar_INPUTS;

value1 /= mm;

value2 *= mm;

break;

case DmassSmass_INPUTS:

input_pair = DmolarSmolar_INPUTS;

value1 /= mm;

value2 *= mm;

break;

case DmassUmass_INPUTS:

input_pair = DmolarUmolar_INPUTS;

value1 /= mm;

value2 *= mm;

break;

default:

break;

}

break;

}

default:

return;

}

if (get_debug_level() >= 50)

std::cout << format("AbstractState: trivial_keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;

switch (key) {

case imolar_mass:

return molar_mass();

case iacentric_factor:

return acentric_factor();

case igas_constant:

return gas_constant();

case iT_min:

return Tmin();

case iT_triple:

return Ttriple();

case iT_max:

return Tmax();

case iP_max:

return pmax();

case iP_min:

case iP_triple:

return this->p_triple();

case iT_reducing:

return calc_T_reducing();

case irhomolar_reducing:

return calc_rhomolar_reducing();

case iP_reducing:

return calc_p_reducing();

case iP_critical:

return this->p_critical();

case iT_critical:

return this->T_critical();

case irhomolar_critical:

return this->rhomolar_critical();

case irhomass_critical:

return this->rhomass_critical();

case iODP:

return this->calc_ODP();

case iGWP100:

return this->calc_GWP100();

case iGWP20:

return this->calc_GWP20();

case iGWP500:

return this->calc_GWP500();

case ifraction_min:

return this->calc_fraction_min();

case ifraction_max:

return this->calc_fraction_max();

case iT_freeze:

return this->calc_T_freeze();

case iFH:

return this->calc_flame_hazard();

case iHH:

return this->calc_health_hazard();

case iPH:

return this->calc_physical_hazard();

case idipole_moment:

return this->calc_dipole_moment();

default:

throw ValueError(

format("This input [%d: \"%s\"] is not valid for trivial_keyed_output", key, get_parameter_information(key, "short").c_str()));

}

if (get_debug_level() >= 50)

std::cout << format("AbstractState: keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;

// Handle trivial inputs

if (is_trivial_parameter(key)) {

return trivial_keyed_output(key);

}

switch (key) {

case iQ:

return Q();

case iT:

return T();

case iP:

return p();

case iDmolar:

return rhomolar();

case iDmass:

return rhomass();

case iHmolar:

return hmolar();

case iHmolar_residual:

return hmolar_residual();

case iHmass:

return hmass();

case iSmolar:

return smolar();

case iSmolar_residual:

return smolar_residual();

case iSmass:

return smass();

case iUmolar:

return umolar();

case iUmass:

return umass();

case iGmolar:

return gibbsmolar();

case iGmolar_residual:

return gibbsmolar_residual();

case iGmass:

return gibbsmass();

case iHelmholtzmolar:

return helmholtzmolar();

case iHelmholtzmass:

return helmholtzmass();

case iCvmolar:

return cvmolar();

case iCvmass:

return cvmass();

case iCpmolar:

return cpmolar();

case iCp0molar:

return cp0molar();

case iCpmass:

return cpmass();

case iCp0mass:

return cp0mass();

case imolar_mass:

return molar_mass();

case iT_reducing:

return get_reducing_state().T;

case irhomolar_reducing:

return get_reducing_state().rhomolar;

case ispeed_sound:

return speed_sound();

case ialphar:

return alphar();

case ialpha0:

return alpha0();

case idalpha0_ddelta_consttau:

return dalpha0_dDelta();

case id2alpha0_ddelta2_consttau:

return d2alpha0_dDelta2();

case id3alpha0_ddelta3_consttau:

return d3alpha0_dDelta3();

case idalpha0_dtau_constdelta:

return dalpha0_dTau();

case idalphar_ddelta_consttau:

return dalphar_dDelta();

case idalphar_dtau_constdelta:

return dalphar_dTau();

case iBvirial:

return Bvirial();

case idBvirial_dT:

return dBvirial_dT();

case iCvirial:

return Cvirial();

case idCvirial_dT:

return dCvirial_dT();

case iisothermal_compressibility:

return isothermal_compressibility();

case iisobaric_expansion_coefficient:

return isobaric_expansion_coefficient();

case iisentropic_expansion_coefficient:

return isentropic_expansion_coefficient();

case iviscosity:

return viscosity();

case iconductivity:

return conductivity();

case iPrandtl:

return Prandtl();

case isurface_tension:

return surface_tension();

case iPhase:

return phase();

case iZ:

return compressibility_factor();

case iPIP:

return PIP();

case ifundamental_derivative_of_gas_dynamics:

return fundamental_derivative_of_gas_dynamics();

default:

throw ValueError(format("This input [%d: \"%s\"] is not valid for keyed_output", key, get_parameter_information(key, "short").c_str()));

}

if (!ValidNumber(_reducing.T)) {

calc_reducing_state();

}

return _reducing.T;

if (!ValidNumber(_reducing.rhomolar)) {

calc_reducing_state();

}

return _reducing.rhomolar;

double tau = calc_T_reducing()/_T;

double delta = _rhomolar/calc_rhomolar_reducing();

double Ar01 = delta*dalphar_dDelta();

double Ar11 = tau*delta*d2alphar_dDelta_dTau();

double Ar20 = tau*tau*d2alphar_dTau2();

return -3.0*(Ar01-Ar11)/Ar20;

CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), dT_dtau = -pow(T, 2) / Tr,

R = AS.gas_constant(), delta = rho / rhor, tau = Tr / T;

switch (index) {

case iT:

dT = 1;

drho = 0;

break;

case iDmolar:

dT = 0;

drho = 1;

break;

case iDmass:

dT = 0;

drho = AS.molar_mass();

break;

case iP: {

// dp/drho|T

drho = R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());

// dp/dT|rho

dT = rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau());

break;

}

case iHmass:

case iHmolar: {

// dh/dT|rho

dT = R

* (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())

+ (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));

// dh/drhomolar|T

drho = T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());

if (index == iHmass) {

// dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case iSmass:

case iSmolar: {

// ds/dT|rho

dT = R / T * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));

// ds/drho|T

drho = R / rho * (-(1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));

if (index == iSmass) {

// ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case iUmass:

case iUmolar: {

// du/dT|rho

dT = R * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));

// du/drho|T

drho = AS.T() * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau());

if (index == iUmass) {

// du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case iGmass:

case iGmolar: {

// dg/dT|rho

double dTau_dT = 1 / dT_dtau;

dT = R * AS.T() * (AS.dalpha0_dTau() + AS.dalphar_dTau() + AS.delta() * AS.d2alphar_dDelta_dTau()) * dTau_dT

+ R * (1 + AS.alpha0() + AS.alphar() + AS.delta() * AS.dalphar_dDelta());

// dg/drho|T

double dDelta_drho = 1 / rhor;

drho = AS.T() * R * (AS.dalpha0_dDelta() + AS.dalphar_dDelta() + AS.delta() * AS.d2alphar_dDelta2() + AS.dalphar_dDelta()) * dDelta_drho;

if (index == iGmass) {

// dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case iTau:

dT = 1 / dT_dtau;

drho = 0;

break;

case iDelta:

dT = 0;

drho = 1 / rhor;

break;

case iCvmolar:

case iCvmass: {

// use the second order derivative of internal energy

// make it cleaner by using the function get_dT_drho_second_derivatives directly?

// dcvdT|rho = d2u/dT2|rho

dT = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));

// dcvdrho|T = d2u/dT/drho

drho = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());

if (index == iCvmass) {

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case iCpmolar:

case iCpmass: {

// dcp/dT|rho = d2h/dT2 + dh/drho * dP/dT * d2P/drhodT / ( dp/drho )^2 - ( d2h/dTdrho * dP/dT + dh/drho * d2P/dT2 ) / ( dP/drho )

dT = R / T * pow(tau, 2)

* (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())

+ delta * AS.d3alphar_dDelta_dTau2());

dT += (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))

* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))

* (R

* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()

- tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau()))

/ pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);

dT -= ((R / rho * delta

* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()

- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()))

* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))

+ (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))

* (rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2())))

/ (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));

// dcpdrho|T = d2h/dTdrho + dh/drho * dP/dT * d2P/drho2 / ( dp/drho )^2 - ( d2h/drho2 * dP/dT + dh/drho * d2P/dTdrho ) / ( dP/drho )

drho = R / rho * delta

* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()

- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()); //d2h/dTdrho

drho +=

(T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))

* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))

* (T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3()))

/ pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);

drho -= ((R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3()))

* (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))

+ (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))

* (R

* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()

- 2 * delta * tau * AS.d2alphar_dDelta_dTau() - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau())))

/ (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));

if (index == iCpmass) {

drho /= AS.molar_mass();

dT /= AS.molar_mass();

}

break;

}

case ispeed_sound: {

//dwdT

double aa = 1.0 + delta * AS.dalphar_dDelta() - delta * tau * AS.d2alphar_dDelta_dTau();

double bb = pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());

double daa_dTau = -delta * tau * AS.d3alphar_dDelta_dTau2();

double dbb_dTau = pow(tau, 2) * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2.0 * tau * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());

double w = AS.speed_sound();

dT = 1.0 / 2.0 / w / T

* (pow(w, 2)

- R * Tr / AS.molar_mass()

* (2.0 * delta * AS.d2alphar_dDelta_dTau() + pow(delta, 2) * AS.d3alphar_dDelta2_dTau()

- (2 * aa / bb * daa_dTau - pow(aa / bb, 2) * dbb_dTau)));

//dwdrho

double daa_dDelta =

AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2() - tau * (AS.d2alphar_dDelta_dTau() + delta * AS.d3alphar_dDelta2_dTau());

double dbb_dDelta = pow(tau, 2) * (AS.d3alpha0_dDelta_dTau2() + AS.d3alphar_dDelta_dTau2());

drho = R * T / 2.0 / AS.molar_mass() / w / rhor

* (2.0 * (AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2())

+ (2.0 * delta * AS.d2alphar_dDelta2() + pow(delta, 2) * AS.d3alphar_dDelta3())

- (2 * aa / bb * daa_dDelta - pow(aa / bb, 2) * dbb_dDelta));

break;

}

default:

throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));

}

CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), R = AS.gas_constant(), delta = rho / rhor,

tau = Tr / T;

// Here we use T and rho as independent variables since derivations are already done by Thorade, 2013,

// Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z

switch (index) {

case iT:

case iDmass:

case iDmolar:

dT2 = 0; // d2rhomolar_dtau2

drho2 = 0;

drho_dT = 0;

break;

case iTau:

dT2 = 2 * Tr / pow(T, 3);

drho_dT = 0;

drho2 = 0;

break;

case iDelta:

dT2 = 0;

drho_dT = 0;

drho2 = 0;

break;

case iP: {

drho2 =

T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3());

dT2 = rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());

drho_dT = R

* (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()

- tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());

break;

}

case iHmass:

case iHmolar: {

// d2h/drho2|T

drho2 = R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3());

// d2h/dT2|rho

dT2 = R / T * pow(tau, 2)

* (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())

+ delta * AS.d3alphar_dDelta_dTau2());

// d2h/drho/dT

drho_dT = R / rho * delta

* (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()

- tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau());

if (index == iHmass) {

drho2 /= AS.molar_mass();

drho_dT /= AS.molar_mass();

dT2 /= AS.molar_mass();

}

break;

}

case iSmass:

case iSmolar: {

// d2s/rho2|T

drho2 = R / pow(rho, 2) * (1 - pow(delta, 2) * AS.d2alphar_dDelta2() + tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());

// d2s/dT2|rho

dT2 = R * pow(tau / T, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 3 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));

// d2s/drho/dT

drho_dT = R / (T * rho) * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());

if (index == iSmass) {

drho2 /= AS.molar_mass();

drho_dT /= AS.molar_mass();

dT2 /= AS.molar_mass();

}

break;

}

case iUmass:

case iUmolar: {

// d2u/rho2|T

drho2 = R * T * tau * pow(delta / rho, 2) * AS.d3alphar_dDelta2_dTau();

// d2u/dT2|rho

dT2 = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));

// d2u/drho/dT

drho_dT = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());

if (index == iUmass) {

drho2 /= AS.molar_mass();