CoolProp  6.6.1dev
An open-source fluid property and humid air property database
AbstractState.cpp
Go to the documentation of this file.
1 /*
2  * AbstractState.cpp
3  *
4  * Created on: 21 Dec 2013
5  * Author: jowr
6  */
7 
8 #ifndef _CRT_SECURE_NO_WARNINGS
9 #define _CRT_SECURE_NO_WARNINGS
10 #endif
11 
12 #include <stdlib.h>
13 #include "math.h"
14 #include "AbstractState.h"
15 #include "DataStructures.h"
21 
22 #if !defined(NO_TABULAR_BACKENDS)
25 #endif
26 
27 namespace CoolProp {
28 
31 
33 {
34  private:
35  std::map<backend_families, shared_ptr<AbstractStateGenerator>> backends;
36 
37  public:
38  void add_backend(const backend_families& bg, const shared_ptr<AbstractStateGenerator>& asg) {
39  backends[bg] = asg;
40  };
42  std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& generator,
43  std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& end) {
44  generator = backends.find(bg);
45  end = backends.end();
46  };
47  std::size_t size() {
48  return backends.size();
49  };
50 };
52  static BackendLibrary the_library;
53  return the_library;
54 }
55 
56 void register_backend(const backend_families& bf, shared_ptr<AbstractStateGenerator> gen) {
58 };
59 
61 {
62  public:
63  AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
64  if (fluid_names.size() == 1) { // Check that fluid_names[0] has only one component
65  std::string str = fluid_names[0]; // Check that the fluid name is an alias for "Water"
66  if ((upper(str) == "WATER") || (upper(str) == "H2O")) {
67  return new IF97Backend();
68  } else {
69  throw ValueError(format("The IF97 backend returns Water props only; fluid name [%s] not allowed", fluid_names[0].c_str()));
70  }
71  } else {
72  throw ValueError(format("The IF97 backend does not support mixtures, only Water"));
73  };
74  };
75 };
76 // This static initialization will cause the generator to register
77 static GeneratorInitializer<IF97BackendGenerator> if97_gen(IF97_BACKEND_FAMILY);
79 {
80  public:
81  AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
82  return new SRKBackend(fluid_names, get_config_double(R_U_CODATA));
83  };
84 };
85 static GeneratorInitializer<SRKGenerator> srk_gen(CoolProp::SRK_BACKEND_FAMILY);
87 {
88  public:
89  AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
90  return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA));
91  };
92 };
93 static GeneratorInitializer<PRGenerator> pr_gen(CoolProp::PR_BACKEND_FAMILY);
95 {
96  public:
97  AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
98  if (fluid_names.size() != 1) {
99  throw ValueError(format("For INCOMP backend, name vector must be one element long"));
100  }
101  return new IncompressibleBackend(fluid_names[0]);
102  };
103 };
104 // This static initialization will cause the generator to register
105 static GeneratorInitializer<IncompressibleBackendGenerator> incomp_gen(INCOMP_BACKEND_FAMILY);
107 {
108  public:
109  CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
110  return new CoolProp::VTPRBackend(fluid_names, CoolProp::get_config_double(R_U_CODATA));
111  };
112 };
113 // This static initialization will cause the generator to register
115 
117 {
118  public:
119  CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
120  return new CoolProp::PCSAFTBackend(fluid_names);
121  };
122 };
123 // This static initialization will cause the generator to register
125 
126 AbstractState* AbstractState::factory(const std::string& backend, const std::vector<std::string>& fluid_names) {
127  if (get_debug_level() > 0) {
128  std::cout << "AbstractState::factory(" << backend << "," << stringvec_to_string(fluid_names) << ")" << std::endl;
129  }
130 
131  backend_families f1;
132  std::string f2;
133  extract_backend_families_string(backend, f1, f2);
134 
135  std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator gen, end;
137 
138  if (get_debug_level() > 0) {
139  std::cout << "AbstractState::factory backend_library size: " << get_backend_library().size() << std::endl;
140  }
141 
142  if (gen != end) {
143  // One of the registered backends was able to match the given backend family
144  return gen->second->get_AbstractState(fluid_names);
145  }
146 #if !defined(NO_TABULAR_BACKENDS)
147  else if (f1 == TTSE_BACKEND_FAMILY) {
148  // Will throw if there is a problem with this backend
149  shared_ptr<AbstractState> AS(factory(f2, fluid_names));
150  return new TTSEBackend(AS);
151  } else if (f1 == BICUBIC_BACKEND_FAMILY) {
152  // Will throw if there is a problem with this backend
153  shared_ptr<AbstractState> AS(factory(f2, fluid_names));
154  return new BicubicBackend(AS);
155  }
156 #endif
157  else if (!backend.compare("?") || backend.empty()) {
158  std::size_t idel = fluid_names[0].find("::");
159  // Backend has not been specified, and we have to figure out what the backend is by parsing the string
160  if (idel == std::string::npos) // No '::' found, no backend specified, try HEOS, otherwise a failure
161  {
162  // Figure out what backend to use
163  return factory("HEOS", fluid_names);
164  } else {
165  // Split string at the '::' into two std::string, call again
166  return factory(std::string(fluid_names[0].begin(), fluid_names[0].begin() + idel),
167  std::string(fluid_names[0].begin() + (idel + 2), fluid_names[0].end()));
168  }
169  } else {
170  throw ValueError(format("Invalid backend name [%s] to factory function", backend.c_str()));
171  }
172 }
173 std::vector<std::string> AbstractState::fluid_names(void) {
174  return calc_fluid_names();
175 }
177  // Reset all instances of CachedElement and overwrite
178  // the internal double values with -_HUGE
179  this->_R = _HUGE;
180  this->_gas_constant.clear();
181  this->_molar_mass.clear();
182  this->_critical.fill(_HUGE);
183  this->_reducing.fill(_HUGE);
184  return true;
185 }
187  // Reset all instances of CachedElement and overwrite
188  // the internal double values with -_HUGE
189  this->_R = _HUGE;
190  this->_gas_constant.clear();
191  this->_molar_mass.clear();
192 
194  this->_rhoLanc.clear();
195  this->_rhoVanc.clear();
196  this->_pVanc.clear();
197  this->_pLanc.clear();
198  this->_TVanc.clear();
199  this->_TLanc.clear();
200 
201  this->_critical.fill(_HUGE);
202  this->_reducing.fill(_HUGE);
203 
205  this->_rhomolar = -_HUGE;
206  this->_T = -_HUGE;
207  this->_p = -_HUGE;
208  this->_Q = -_HUGE;
209  this->_tau.clear();
210  this->_delta.clear();
211 
212  this->_umolar.clear();
213  this->_cpmolar.clear();
214  this->_cp0molar.clear();
215  this->_cvmolar.clear();
216  this->_speed_sound.clear();
217  this->_hmolar.clear();
218  this->_smolar.clear();
219  this->_gibbsmolar.clear();
220  this->_helmholtzmolar.clear();
221  this->_logp.clear();
222  this->_logrhomolar.clear();
223 
224  this->_hmolar_excess.clear();
225  this->_smolar_excess.clear();
226  this->_gibbsmolar_excess.clear();
227  this->_volumemolar_excess.clear();
228  this->_umolar_excess.clear();
230 
231  this->_hmolar_residual.clear();
232  this->_smolar_residual.clear();
233  this->_gibbsmolar_residual.clear();
234 
236  this->_rho_spline.clear();
239 
241  this->_alpha0.clear();
242  this->_dalpha0_dTau.clear();
243  this->_dalpha0_dDelta.clear();
244  this->_d2alpha0_dTau2.clear();
246  this->_d2alpha0_dDelta2.clear();
247  this->_d3alpha0_dTau3.clear();
250  this->_d3alpha0_dDelta3.clear();
251  this->_alphar.clear();
252  this->_dalphar_dTau.clear();
253  this->_dalphar_dDelta.clear();
254  this->_d2alphar_dTau2.clear();
256  this->_d2alphar_dDelta2.clear();
257  this->_d3alphar_dTau3.clear();
260  this->_d3alphar_dDelta3.clear();
261 
262  this->_dalphar_dDelta_lim.clear();
266 
268  this->_rhoLmolar.clear();
269  this->_rhoVmolar.clear();
270 
272  this->_viscosity.clear();
273  this->_conductivity.clear();
274  this->_surface_tension.clear();
275 
276  return true;
277 }
279  // Check if a mass based input, convert it to molar units
280 
281  switch (input_pair) {
282  case DmassT_INPUTS:
283  //case HmassT_INPUTS: ///< Enthalpy in J/kg, Temperature in K (NOT CURRENTLY IMPLEMENTED)
284  case SmassT_INPUTS:
285  //case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)
286  case DmassP_INPUTS:
287  case DmassQ_INPUTS:
288  case HmassP_INPUTS:
289  case PSmass_INPUTS:
290  case PUmass_INPUTS:
291  case HmassSmass_INPUTS:
292  case SmassUmass_INPUTS:
293  case DmassHmass_INPUTS:
294  case DmassSmass_INPUTS:
295  case DmassUmass_INPUTS:
296  {
297  // Set the cache value for the molar mass if it hasn't been set yet
298  molar_mass();
299 
300  // Molar mass (just for compactness of the following switch)
301  CoolPropDbl mm = static_cast<CoolPropDbl>(_molar_mass);
302 
303  switch (input_pair) {
304  case DmassT_INPUTS:
305  input_pair = DmolarT_INPUTS;
306  value1 /= mm;
307  break;
308  //case HmassT_INPUTS: input_pair = HmolarT_INPUTS; value1 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
309  case SmassT_INPUTS:
310  input_pair = SmolarT_INPUTS;
311  value1 *= mm;
312  break;
313  //case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
314  case DmassP_INPUTS:
315  input_pair = DmolarP_INPUTS;
316  value1 /= mm;
317  break;
318  case DmassQ_INPUTS:
319  input_pair = DmolarQ_INPUTS;
320  value1 /= mm;
321  break;
322  case HmassP_INPUTS:
323  input_pair = HmolarP_INPUTS;
324  value1 *= mm;
325  break;
326  case PSmass_INPUTS:
327  input_pair = PSmolar_INPUTS;
328  value2 *= mm;
329  break;
330  case PUmass_INPUTS:
331  input_pair = PUmolar_INPUTS;
332  value2 *= mm;
333  break;
334  case HmassSmass_INPUTS:
335  input_pair = HmolarSmolar_INPUTS;
336  value1 *= mm;
337  value2 *= mm;
338  break;
339  case SmassUmass_INPUTS:
340  input_pair = SmolarUmolar_INPUTS;
341  value1 *= mm;
342  value2 *= mm;
343  break;
344  case DmassHmass_INPUTS:
345  input_pair = DmolarHmolar_INPUTS;
346  value1 /= mm;
347  value2 *= mm;
348  break;
349  case DmassSmass_INPUTS:
350  input_pair = DmolarSmolar_INPUTS;
351  value1 /= mm;
352  value2 *= mm;
353  break;
354  case DmassUmass_INPUTS:
355  input_pair = DmolarUmolar_INPUTS;
356  value1 /= mm;
357  value2 *= mm;
358  break;
359  default:
360  break;
361  }
362  break;
363  }
364  default:
365  return;
366  }
367 }
369  if (get_debug_level() >= 50)
370  std::cout << format("AbstractState: trivial_keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
371  switch (key) {
372  case imolar_mass:
373  return molar_mass();
374  case iacentric_factor:
375  return acentric_factor();
376  case igas_constant:
377  return gas_constant();
378  case iT_min:
379  return Tmin();
380  case iT_triple:
381  return Ttriple();
382  case iT_max:
383  return Tmax();
384  case iP_max:
385  return pmax();
386  case iP_min:
387  case iP_triple:
388  return this->p_triple();
389  case iT_reducing:
390  return calc_T_reducing();
391  case irhomolar_reducing:
392  return calc_rhomolar_reducing();
393  case iP_reducing:
394  return calc_p_reducing();
395  case iP_critical:
396  return this->p_critical();
397  case iT_critical:
398  return this->T_critical();
399  case irhomolar_critical:
400  return this->rhomolar_critical();
401  case irhomass_critical:
402  return this->rhomass_critical();
403  case iODP:
404  return this->calc_ODP();
405  case iGWP100:
406  return this->calc_GWP100();
407  case iGWP20:
408  return this->calc_GWP20();
409  case iGWP500:
410  return this->calc_GWP500();
411  case ifraction_min:
412  return this->calc_fraction_min();
413  case ifraction_max:
414  return this->calc_fraction_max();
415  case iT_freeze:
416  return this->calc_T_freeze();
417  case iFH:
418  return this->calc_flame_hazard();
419  case iHH:
420  return this->calc_health_hazard();
421  case iPH:
422  return this->calc_physical_hazard();
423  case idipole_moment:
424  return this->calc_dipole_moment();
425  default:
426  throw ValueError(
427  format("This input [%d: \"%s\"] is not valid for trivial_keyed_output", key, get_parameter_information(key, "short").c_str()));
428  }
429 }
431  if (get_debug_level() >= 50)
432  std::cout << format("AbstractState: keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
433  // Handle trivial inputs
434  if (is_trivial_parameter(key)) {
435  return trivial_keyed_output(key);
436  }
437  switch (key) {
438  case iQ:
439  return Q();
440  case iT:
441  return T();
442  case iP:
443  return p();
444  case iDmolar:
445  return rhomolar();
446  case iDmass:
447  return rhomass();
448  case iHmolar:
449  return hmolar();
450  case iHmolar_residual:
451  return hmolar_residual();
452  case iHmass:
453  return hmass();
454  case iSmolar:
455  return smolar();
456  case iSmolar_residual:
457  return smolar_residual();
458  case iSmass:
459  return smass();
460  case iUmolar:
461  return umolar();
462  case iUmass:
463  return umass();
464  case iGmolar:
465  return gibbsmolar();
466  case iGmolar_residual:
467  return gibbsmolar_residual();
468  case iGmass:
469  return gibbsmass();
470  case iHelmholtzmolar:
471  return helmholtzmolar();
472  case iHelmholtzmass:
473  return helmholtzmass();
474  case iCvmolar:
475  return cvmolar();
476  case iCvmass:
477  return cvmass();
478  case iCpmolar:
479  return cpmolar();
480  case iCp0molar:
481  return cp0molar();
482  case iCpmass:
483  return cpmass();
484  case iCp0mass:
485  return cp0mass();
486  case imolar_mass:
487  return molar_mass();
488  case iT_reducing:
489  return get_reducing_state().T;
490  case irhomolar_reducing:
491  return get_reducing_state().rhomolar;
492  case ispeed_sound:
493  return speed_sound();
494  case ialphar:
495  return alphar();
496  case ialpha0:
497  return alpha0();
499  return dalpha0_dDelta();
501  return d2alpha0_dDelta2();
503  return d3alpha0_dDelta3();
505  return dalpha0_dTau();
507  return dalphar_dDelta();
509  return dalphar_dTau();
510  case iBvirial:
511  return Bvirial();
512  case idBvirial_dT:
513  return dBvirial_dT();
514  case iCvirial:
515  return Cvirial();
516  case idCvirial_dT:
517  return dCvirial_dT();
524  case iviscosity:
525  return viscosity();
526  case iconductivity:
527  return conductivity();
528  case iPrandtl:
529  return Prandtl();
530  case isurface_tension:
531  return surface_tension();
532  case iPhase:
533  return phase();
534  case iZ:
535  return compressibility_factor();
536  case iPIP:
537  return PIP();
540  default:
541  throw ValueError(format("This input [%d: \"%s\"] is not valid for keyed_output", key, get_parameter_information(key, "short").c_str()));
542  }
543 }
544 
545 double AbstractState::tau(void) {
547  return _tau;
548 }
549 double AbstractState::delta(void) {
551  return _delta;
552 }
553 double AbstractState::Tmin(void) {
554  return calc_Tmin();
555 }
556 double AbstractState::Tmax(void) {
557  return calc_Tmax();
558 }
560  return calc_Ttriple();
561 }
562 double AbstractState::pmax(void) {
563  return calc_pmax();
564 }
566  return calc_T_critical();
567 }
569  if (!ValidNumber(_reducing.T)) {
571  }
572  return _reducing.T;
573 }
575  return calc_p_critical();
576 }
578  return calc_p_triple();
579 }
581  return calc_rhomolar_critical();
582 }
584  return calc_rhomolar_critical() * molar_mass();
585 }
589  }
590  return _reducing.rhomolar;
591 }
593  return rhomolar_reducing() * molar_mass();
594 }
595 double AbstractState::hmolar(void) {
596  if (!_hmolar) _hmolar = calc_hmolar();
597  return _hmolar;
598 }
601  return _hmolar_residual;
602 }
605  return _hmolar_excess;
606 }
607 double AbstractState::smolar(void) {
608  if (!_smolar) _smolar = calc_smolar();
609  return _smolar;
610 }
613  return _smolar_residual;
614 }
615 double AbstractState::neff(void) {
616  double tau = calc_T_reducing()/_T;
618  double Ar01 = delta*dalphar_dDelta();
619  double Ar11 = tau*delta*d2alphar_dDelta_dTau();
620  double Ar20 = tau*tau*d2alphar_dTau2();
621  return -3.0*(Ar01-Ar11)/Ar20;
622 }
625  return _smolar_excess;
626 }
627 double AbstractState::umolar(void) {
628  if (!_umolar) _umolar = calc_umolar();
629  return _umolar;
630 }
633  return _umolar_excess;
634 }
637  return _gibbsmolar;
638 }
641  return _gibbsmolar_residual;
642 }
645  return _gibbsmolar_excess;
646 }
649  return _helmholtzmolar;
650 }
653  return _helmholtzmolar_excess;
654 }
657  return _volumemolar_excess;
658 }
660  if (!_cpmolar) _cpmolar = calc_cpmolar();
661  return _cpmolar;
662 }
664  return calc_cpmolar_idealgas();
665 }
667  if (!_cvmolar) _cvmolar = calc_cvmolar();
668  return _cvmolar;
669 }
672  return _speed_sound;
673 }
676  return _viscosity;
677 }
680  return _conductivity;
681 }
682 double AbstractState::melting_line(int param, int given, double value) {
683  return calc_melting_line(param, given, value);
684 }
686  return calc_acentric_factor();
687 }
688 double AbstractState::saturation_ancillary(parameters param, int Q, parameters given, double value) {
689  return calc_saturation_ancillary(param, Q, given, value);
690 }
693  return _surface_tension;
694 }
697  return _molar_mass;
698 }
701  return _gas_constant;
702 }
703 double AbstractState::fugacity_coefficient(std::size_t i) {
704  // TODO: Cache the fug. coeff for each component
705  return calc_fugacity_coefficient(i);
706 }
708  // TODO: Cache the fug. coeff for each component
710 }
711 double AbstractState::fugacity(std::size_t i) {
712  // TODO: Cache the fug. coeff for each component
713  return calc_fugacity(i);
714 }
715 double AbstractState::chemical_potential(std::size_t i) {
716  // TODO: Cache the chemical potential for each component
717  return calc_chemical_potential(i);
718 }
719 void AbstractState::build_phase_envelope(const std::string& type) {
720  calc_phase_envelope(type);
721 }
723  return 1.0 / _rhomolar * first_partial_deriv(iDmolar, iP, iT);
724 }
726  return -1.0 / _rhomolar * first_partial_deriv(iDmolar, iT, iP);
727 }
730 }
732  return calc_Bvirial();
733 }
735  return calc_Cvirial();
736 }
738  return calc_dBvirial_dT();
739 }
741  return calc_dCvirial_dT();
742 }
745 }
746 
748  // See Colonna, FPE, 2010, Eq. 1
749  return 1 + this->second_partial_deriv(iP, iDmass, iSmolar, iDmass, iSmolar) * this->rhomass() / (2 * powInt(speed_sound(), 2));
750 };
751 
752 // Get the derivatives of the parameters in the partial derivative with respect to T and rho
754  CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), dT_dtau = -pow(T, 2) / Tr,
755  R = AS.gas_constant(), delta = rho / rhor, tau = Tr / T;
756 
757  switch (index) {
758  case iT:
759  dT = 1;
760  drho = 0;
761  break;
762  case iDmolar:
763  dT = 0;
764  drho = 1;
765  break;
766  case iDmass:
767  dT = 0;
768  drho = AS.molar_mass();
769  break;
770  case iP: {
771  // dp/drho|T
772  drho = R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
773  // dp/dT|rho
774  dT = rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau());
775  break;
776  }
777  case iHmass:
778  case iHmolar: {
779  // dh/dT|rho
780  dT = R
781  * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
782  + (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
783  // dh/drhomolar|T
784  drho = T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
785  if (index == iHmass) {
786  // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
787  drho /= AS.molar_mass();
788  dT /= AS.molar_mass();
789  }
790  break;
791  }
792  case iSmass:
793  case iSmolar: {
794  // ds/dT|rho
795  dT = R / T * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
796  // ds/drho|T
797  drho = R / rho * (-(1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
798  if (index == iSmass) {
799  // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
800  drho /= AS.molar_mass();
801  dT /= AS.molar_mass();
802  }
803  break;
804  }
805  case iUmass:
806  case iUmolar: {
807  // du/dT|rho
808  dT = R * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
809  // du/drho|T
810  drho = AS.T() * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau());
811  if (index == iUmass) {
812  // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
813  drho /= AS.molar_mass();
814  dT /= AS.molar_mass();
815  }
816  break;
817  }
818  case iGmass:
819  case iGmolar: {
820  // dg/dT|rho
821  double dTau_dT = 1 / dT_dtau;
822  dT = R * AS.T() * (AS.dalpha0_dTau() + AS.dalphar_dTau() + AS.delta() * AS.d2alphar_dDelta_dTau()) * dTau_dT
823  + R * (1 + AS.alpha0() + AS.alphar() + AS.delta() * AS.dalphar_dDelta());
824  // dg/drho|T
825  double dDelta_drho = 1 / rhor;
826  drho = AS.T() * R * (AS.dalpha0_dDelta() + AS.dalphar_dDelta() + AS.delta() * AS.d2alphar_dDelta2() + AS.dalphar_dDelta()) * dDelta_drho;
827  if (index == iGmass) {
828  // dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
829  drho /= AS.molar_mass();
830  dT /= AS.molar_mass();
831  }
832  break;
833  }
834  case iTau:
835  dT = 1 / dT_dtau;
836  drho = 0;
837  break;
838  case iDelta:
839  dT = 0;
840  drho = 1 / rhor;
841  break;
842  case iCvmolar:
843  case iCvmass: {
844  // use the second order derivative of internal energy
845  // make it cleaner by using the function get_dT_drho_second_derivatives directly?
846  // dcvdT|rho = d2u/dT2|rho
847  dT = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
848  // dcvdrho|T = d2u/dT/drho
849  drho = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
850  if (index == iCvmass) {
851  drho /= AS.molar_mass();
852  dT /= AS.molar_mass();
853  }
854  break;
855  }
856  case iCpmolar:
857  case iCpmass: {
858  // dcp/dT|rho = d2h/dT2 + dh/drho * dP/dT * d2P/drhodT / ( dp/drho )^2 - ( d2h/dTdrho * dP/dT + dh/drho * d2P/dT2 ) / ( dP/drho )
859  dT = R / T * pow(tau, 2)
860  * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
861  + delta * AS.d3alphar_dDelta_dTau2());
862  dT += (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
863  * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
864  * (R
865  * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
866  - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau()))
867  / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
868  dT -= ((R / rho * delta
869  * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
870  - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()))
871  * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
872  + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
873  * (rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2())))
874  / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
875  // dcpdrho|T = d2h/dTdrho + dh/drho * dP/dT * d2P/drho2 / ( dp/drho )^2 - ( d2h/drho2 * dP/dT + dh/drho * d2P/dTdrho ) / ( dP/drho )
876  drho = R / rho * delta
877  * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
878  - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()); //d2h/dTdrho
879  drho +=
880  (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
881  * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
882  * (T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3()))
883  / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
884  drho -= ((R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3()))
885  * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
886  + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
887  * (R
888  * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()
889  - 2 * delta * tau * AS.d2alphar_dDelta_dTau() - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau())))
890  / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
891  if (index == iCpmass) {
892  drho /= AS.molar_mass();
893  dT /= AS.molar_mass();
894  }
895  break;
896  }
897  case ispeed_sound: {
898  //dwdT
899  double aa = 1.0 + delta * AS.dalphar_dDelta() - delta * tau * AS.d2alphar_dDelta_dTau();
900  double bb = pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
901  double daa_dTau = -delta * tau * AS.d3alphar_dDelta_dTau2();
902  double dbb_dTau = pow(tau, 2) * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2.0 * tau * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
903  double w = AS.speed_sound();
904  dT = 1.0 / 2.0 / w / T
905  * (pow(w, 2)
906  - R * Tr / AS.molar_mass()
907  * (2.0 * delta * AS.d2alphar_dDelta_dTau() + pow(delta, 2) * AS.d3alphar_dDelta2_dTau()
908  - (2 * aa / bb * daa_dTau - pow(aa / bb, 2) * dbb_dTau)));
909  //dwdrho
910  double daa_dDelta =
911  AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2() - tau * (AS.d2alphar_dDelta_dTau() + delta * AS.d3alphar_dDelta2_dTau());
912  double dbb_dDelta = pow(tau, 2) * (AS.d3alpha0_dDelta_dTau2() + AS.d3alphar_dDelta_dTau2());
913  drho = R * T / 2.0 / AS.molar_mass() / w / rhor
914  * (2.0 * (AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2())
915  + (2.0 * delta * AS.d2alphar_dDelta2() + pow(delta, 2) * AS.d3alphar_dDelta3())
916  - (2 * aa / bb * daa_dDelta - pow(aa / bb, 2) * dbb_dDelta));
917  break;
918  }
919  default:
920  throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
921  }
922 }
924  CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), R = AS.gas_constant(), delta = rho / rhor,
925  tau = Tr / T;
926 
927  // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013,
928  // Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z
929 
930  switch (index) {
931  case iT:
932  case iDmass:
933  case iDmolar:
934  dT2 = 0; // d2rhomolar_dtau2
935  drho2 = 0;
936  drho_dT = 0;
937  break;
938  case iTau:
939  dT2 = 2 * Tr / pow(T, 3);
940  drho_dT = 0;
941  drho2 = 0;
942  break;
943  case iDelta:
944  dT2 = 0;
945  drho_dT = 0;
946  drho2 = 0;
947  break;
948  case iP: {
949  drho2 =
950  T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3());
951  dT2 = rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
952  drho_dT = R
953  * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
954  - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
955  break;
956  }
957  case iHmass:
958  case iHmolar: {
959  // d2h/drho2|T
960  drho2 = R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3());
961  // d2h/dT2|rho
962  dT2 = R / T * pow(tau, 2)
963  * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
964  + delta * AS.d3alphar_dDelta_dTau2());
965  // d2h/drho/dT
966  drho_dT = R / rho * delta
967  * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
968  - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau());
969  if (index == iHmass) {
970  drho2 /= AS.molar_mass();
971  drho_dT /= AS.molar_mass();
972  dT2 /= AS.molar_mass();
973  }
974  break;
975  }
976  case iSmass:
977  case iSmolar: {
978  // d2s/rho2|T
979  drho2 = R / pow(rho, 2) * (1 - pow(delta, 2) * AS.d2alphar_dDelta2() + tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
980  // d2s/dT2|rho
981  dT2 = R * pow(tau / T, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 3 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
982  // d2s/drho/dT
983  drho_dT = R / (T * rho) * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
984  if (index == iSmass) {
985  drho2 /= AS.molar_mass();
986  drho_dT /= AS.molar_mass();
987  dT2 /= AS.molar_mass();
988  }
989  break;
990  }
991  case iUmass:
992  case iUmolar: {
993  // d2u/rho2|T
994  drho2 = R * T * tau * pow(delta / rho, 2) * AS.d3alphar_dDelta2_dTau();
995  // d2u/dT2|rho
996  dT2 = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
997  // d2u/drho/dT
998  drho_dT = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
999  if (index == iUmass) {
1000  drho2 /= AS.molar_mass();
1001  drho_dT /= AS.molar_mass();
1002  dT2 /= AS.molar_mass();
1003  }
1004  break;
1005  }
1006  default:
1007  throw ValueError(format("input to get_dT_drho_second_derivatives[%s] is invalid", get_parameter_information(index, "short").c_str()));
1008  }
1009 }
1011  CoolPropDbl dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho;
1012 
1013  get_dT_drho(*this, Of, dOf_dT, dOf_drho);
1014  get_dT_drho(*this, Wrt, dWrt_dT, dWrt_drho);
1015  get_dT_drho(*this, Constant, dConstant_dT, dConstant_drho);
1016 
1017  return (dOf_dT * dConstant_drho - dOf_drho * dConstant_dT) / (dWrt_dT * dConstant_drho - dWrt_drho * dConstant_dT);
1018 }
1020  CoolPropDbl dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2, d2Wrt1_dT2,
1021  d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D,
1022  dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, dderiv1_drho, dderiv1_dT, second;
1023 
1024  // First and second partials needed for terms involved in first derivative
1025  get_dT_drho(*this, Of1, dOf1_dT, dOf1_drho);
1026  get_dT_drho(*this, Wrt1, dWrt1_dT, dWrt1_drho);
1027  get_dT_drho(*this, Constant1, dConstant1_dT, dConstant1_drho);
1028  get_dT_drho_second_derivatives(*this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2);
1029  get_dT_drho_second_derivatives(*this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2);
1030  get_dT_drho_second_derivatives(*this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2);
1031 
1032  // First derivatives of terms involved in the second derivative
1033  get_dT_drho(*this, Wrt2, dWrt2_dT, dWrt2_drho);
1034  get_dT_drho(*this, Constant2, dConstant2_dT, dConstant2_drho);
1035 
1036  // Numerator and denominator of first partial derivative term
1037  N = dOf1_dT * dConstant1_drho - dOf1_drho * dConstant1_dT;
1038  D = dWrt1_dT * dConstant1_drho - dWrt1_drho * dConstant1_dT;
1039 
1040  // Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant
1041  // They are of similar form, with Of1 and Wrt1 swapped
1042  dNdrho__T = dOf1_dT * d2Constant1_drho2 + d2Of1_drhodT * dConstant1_drho - dOf1_drho * d2Constant1_drhodT - d2Of1_drho2 * dConstant1_dT;
1043  dDdrho__T = dWrt1_dT * d2Constant1_drho2 + d2Wrt1_drhodT * dConstant1_drho - dWrt1_drho * d2Constant1_drhodT - d2Wrt1_drho2 * dConstant1_dT;
1044 
1045  // Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant
1046  // They are of similar form, with Of1 and Wrt1 swapped
1047  dNdT__rho = dOf1_dT * d2Constant1_drhodT + d2Of1_dT2 * dConstant1_drho - dOf1_drho * d2Constant1_dT2 - d2Of1_drhodT * dConstant1_dT;
1048  dDdT__rho = dWrt1_dT * d2Constant1_drhodT + d2Wrt1_dT2 * dConstant1_drho - dWrt1_drho * d2Constant1_dT2 - d2Wrt1_drhodT * dConstant1_dT;
1049 
1050  // First partial of first derivative term with respect to T
1051  dderiv1_drho = (D * dNdrho__T - N * dDdrho__T) / pow(D, 2);
1052 
1053  // First partial of first derivative term with respect to rho
1054  dderiv1_dT = (D * dNdT__rho - N * dDdT__rho) / pow(D, 2);
1055 
1056  // Complete second derivative
1057  second = (dderiv1_dT * dConstant2_drho - dderiv1_drho * dConstant2_dT) / (dWrt2_dT * dConstant2_drho - dWrt2_drho * dConstant2_dT);
1058 
1059  return second;
1060 }
1061 // // ----------------------------------------
1062 // // Smoothing functions for density
1063 // // ----------------------------------------
1064 // /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1065 // virtual double AbstractState::drhodh_constp_smoothed(double xend);
1066 // /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1067 // virtual double AbstractState::drhodp_consth_smoothed(double xend);
1068 // /// Density corresponding to the smoothed derivatives in the region of x=0 to x=xend
1069 // virtual void AbstractState::rho_smoothed(double xend, double *rho_spline, double *dsplinedh, double *dsplinedp);
1070 
1071 } /* namespace CoolProp */
1072 
1073 #ifdef ENABLE_CATCH
1074 
1075 #include <catch2/catch_all.hpp>
1076 
1077 TEST_CASE("Check AbstractState", "[AbstractState]") {
1078  SECTION("bad backend") {
1079  CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
1080  }
1081  SECTION("good backend - bad fluid") {
1082  CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "DEFINITELY_A_BAD_FLUID")));
1083  }
1084  SECTION("good backend - helmholtz") {
1085  CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Water")));
1086  }
1087  SECTION("good backend - incomp") {
1088  CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("INCOMP", "DEB")));
1089  }
1090  SECTION("good backend - REFPROP") {
1091  CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("REFPROP", "Water")));
1092  }
1093 }
1094 
1095 TEST_CASE("Check derivatives in first_partial_deriv", "[derivs_in_first_partial_deriv]") {
1096  shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("HEOS", "Water"));
1097  shared_ptr<CoolProp::AbstractState> WaterplusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1098  shared_ptr<CoolProp::AbstractState> WaterminusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1099  shared_ptr<CoolProp::AbstractState> Waterplusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1100  shared_ptr<CoolProp::AbstractState> Waterminusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1101 
1102  double dT = 1e-3, drho = 1;
1103  Water->update(CoolProp::PT_INPUTS, 101325, 300);
1104  WaterplusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 + dT);
1105  WaterminusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 - dT);
1106  Waterplusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() + drho, 300);
1107  Waterminusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() - drho, 300);
1108 
1109  // Numerical derivatives
1110  CoolPropDbl dP_dT_num = (WaterplusT->p() - WaterminusT->p()) / (2 * dT);
1111  CoolPropDbl dP_drho_num = (Waterplusrho->p() - Waterminusrho->p()) / (2 * drho);
1112 
1113  CoolPropDbl dHmolar_dT_num = (WaterplusT->hmolar() - WaterminusT->hmolar()) / (2 * dT);
1114  CoolPropDbl dHmolar_drho_num = (Waterplusrho->hmolar() - Waterminusrho->hmolar()) / (2 * drho);
1115  CoolPropDbl dHmass_dT_num = (WaterplusT->hmass() - WaterminusT->hmass()) / (2 * dT);
1116  CoolPropDbl dHmass_drho_num = (Waterplusrho->hmass() - Waterminusrho->hmass()) / (2 * drho);
1117 
1118  CoolPropDbl dSmolar_dT_num = (WaterplusT->smolar() - WaterminusT->smolar()) / (2 * dT);
1119  CoolPropDbl dSmolar_drho_num = (Waterplusrho->smolar() - Waterminusrho->smolar()) / (2 * drho);
1120  CoolPropDbl dSmass_dT_num = (WaterplusT->smass() - WaterminusT->smass()) / (2 * dT);
1121  CoolPropDbl dSmass_drho_num = (Waterplusrho->smass() - Waterminusrho->smass()) / (2 * drho);
1122 
1123  CoolPropDbl dUmolar_dT_num = (WaterplusT->umolar() - WaterminusT->umolar()) / (2 * dT);
1124  CoolPropDbl dUmolar_drho_num = (Waterplusrho->umolar() - Waterminusrho->umolar()) / (2 * drho);
1125  CoolPropDbl dUmass_dT_num = (WaterplusT->umass() - WaterminusT->umass()) / (2 * dT);
1126  CoolPropDbl dUmass_drho_num = (Waterplusrho->umass() - Waterminusrho->umass()) / (2 * drho);
1127 
1128  CoolPropDbl dGmolar_dT_num = (WaterplusT->gibbsmolar() - WaterminusT->gibbsmolar()) / (2 * dT);
1129  CoolPropDbl dGmolar_drho_num = (Waterplusrho->gibbsmolar() - Waterminusrho->gibbsmolar()) / (2 * drho);
1130  CoolPropDbl dGmass_dT_num = (WaterplusT->gibbsmass() - WaterminusT->gibbsmass()) / (2 * dT);
1131  CoolPropDbl dGmass_drho_num = (Waterplusrho->gibbsmass() - Waterminusrho->gibbsmass()) / (2 * drho);
1132 
1133  CoolPropDbl dCvmolar_dT_num = (WaterplusT->cvmolar() - WaterminusT->cvmolar()) / (2 * dT);
1134  CoolPropDbl dCvmolar_drho_num = (Waterplusrho->cvmolar() - Waterminusrho->cvmolar()) / (2 * drho);
1135  CoolPropDbl dCvmass_dT_num = (WaterplusT->cvmass() - WaterminusT->cvmass()) / (2 * dT);
1136  CoolPropDbl dCvmass_drho_num = (Waterplusrho->cvmass() - Waterminusrho->cvmass()) / (2 * drho);
1137 
1138  CoolPropDbl dCpmolar_dT_num = (WaterplusT->cpmolar() - WaterminusT->cpmolar()) / (2 * dT);
1139  CoolPropDbl dCpmolar_drho_num = (Waterplusrho->cpmolar() - Waterminusrho->cpmolar()) / (2 * drho);
1140  CoolPropDbl dCpmass_dT_num = (WaterplusT->cpmass() - WaterminusT->cpmass()) / (2 * dT);
1141  CoolPropDbl dCpmass_drho_num = (Waterplusrho->cpmass() - Waterminusrho->cpmass()) / (2 * drho);
1142 
1143  CoolPropDbl dspeed_sound_dT_num = (WaterplusT->speed_sound() - WaterminusT->speed_sound()) / (2 * dT);
1144  CoolPropDbl dspeed_sound_drho_num = (Waterplusrho->speed_sound() - Waterminusrho->speed_sound()) / (2 * drho);
1145 
1146  // Pressure
1147  CoolPropDbl dP_dT_analyt, dP_drho_analyt;
1148  CoolProp::get_dT_drho(*Water, CoolProp::iP, dP_dT_analyt, dP_drho_analyt);
1149  // Enthalpy
1150  CoolPropDbl dHmolar_dT_analyt, dHmolar_drho_analyt;
1151  CoolProp::get_dT_drho(*Water, CoolProp::iHmolar, dHmolar_dT_analyt, dHmolar_drho_analyt);
1152  CoolPropDbl dHmass_dT_analyt, dHmass_drho_analyt;
1153  CoolProp::get_dT_drho(*Water, CoolProp::iHmass, dHmass_dT_analyt, dHmass_drho_analyt);
1154  // Entropy
1155  CoolPropDbl dSmolar_dT_analyt, dSmolar_drho_analyt;
1156  CoolProp::get_dT_drho(*Water, CoolProp::iSmolar, dSmolar_dT_analyt, dSmolar_drho_analyt);
1157  CoolPropDbl dSmass_dT_analyt, dSmass_drho_analyt;
1158  CoolProp::get_dT_drho(*Water, CoolProp::iSmass, dSmass_dT_analyt, dSmass_drho_analyt);
1159  // Internal energy
1160  CoolPropDbl dUmolar_dT_analyt, dUmolar_drho_analyt;
1161  CoolProp::get_dT_drho(*Water, CoolProp::iUmolar, dUmolar_dT_analyt, dUmolar_drho_analyt);
1162  CoolPropDbl dUmass_dT_analyt, dUmass_drho_analyt;
1163  CoolProp::get_dT_drho(*Water, CoolProp::iUmass, dUmass_dT_analyt, dUmass_drho_analyt);
1164  // Gibbs
1165  CoolPropDbl dGmolar_dT_analyt, dGmolar_drho_analyt;
1166  CoolProp::get_dT_drho(*Water, CoolProp::iGmolar, dGmolar_dT_analyt, dGmolar_drho_analyt);
1167  CoolPropDbl dGmass_dT_analyt, dGmass_drho_analyt;
1168  CoolProp::get_dT_drho(*Water, CoolProp::iGmass, dGmass_dT_analyt, dGmass_drho_analyt);
1169  // Isochoric heat capacity
1170  CoolPropDbl dCvmolar_dT_analyt, dCvmolar_drho_analyt;
1171  CoolProp::get_dT_drho(*Water, CoolProp::iCvmolar, dCvmolar_dT_analyt, dCvmolar_drho_analyt);
1172  CoolPropDbl dCvmass_dT_analyt, dCvmass_drho_analyt;
1173  CoolProp::get_dT_drho(*Water, CoolProp::iCvmass, dCvmass_dT_analyt, dCvmass_drho_analyt);
1174  // Isobaric heat capacity
1175  CoolPropDbl dCpmolar_dT_analyt, dCpmolar_drho_analyt;
1176  CoolProp::get_dT_drho(*Water, CoolProp::iCpmolar, dCpmolar_dT_analyt, dCpmolar_drho_analyt);
1177  CoolPropDbl dCpmass_dT_analyt, dCpmass_drho_analyt;
1178  CoolProp::get_dT_drho(*Water, CoolProp::iCpmass, dCpmass_dT_analyt, dCpmass_drho_analyt);
1179  // Speed of sound
1180  CoolPropDbl dspeed_sound_dT_analyt, dspeed_sound_drho_analyt;
1181  CoolProp::get_dT_drho(*Water, CoolProp::ispeed_sound, dspeed_sound_dT_analyt, dspeed_sound_drho_analyt);
1182 
1183  double eps = 1e-3;
1184 
1185  CHECK(std::abs(dP_dT_analyt / dP_dT_num - 1) < eps);
1186  CHECK(std::abs(dP_drho_analyt / dP_drho_num - 1) < eps);
1187 
1188  CHECK(std::abs(dHmolar_dT_analyt / dHmolar_dT_num - 1) < eps);
1189  CHECK(std::abs(dHmolar_drho_analyt / dHmolar_drho_num - 1) < eps);
1190  CHECK(std::abs(dHmass_dT_analyt / dHmass_dT_num - 1) < eps);
1191  CHECK(std::abs(dHmass_drho_analyt / dHmass_drho_num - 1) < eps);
1192 
1193  CHECK(std::abs(dSmolar_dT_analyt / dSmolar_dT_num - 1) < eps);
1194  CHECK(std::abs(dSmolar_drho_analyt / dSmolar_drho_num - 1) < eps);
1195  CHECK(std::abs(dSmass_dT_analyt / dSmass_dT_num - 1) < eps);
1196  CHECK(std::abs(dSmass_drho_analyt / dSmass_drho_num - 1) < eps);
1197 
1198  CHECK(std::abs(dUmolar_dT_analyt / dUmolar_dT_num - 1) < eps);
1199  CHECK(std::abs(dUmolar_drho_analyt / dUmolar_drho_num - 1) < eps);
1200  CHECK(std::abs(dUmass_dT_analyt / dUmass_dT_num - 1) < eps);
1201  CHECK(std::abs(dUmass_drho_analyt / dUmass_drho_num - 1) < eps);
1202 
1203  CHECK(std::abs(dGmolar_dT_analyt / dGmolar_dT_num - 1) < eps);
1204  CHECK(std::abs(dGmolar_drho_analyt / dGmolar_drho_num - 1) < eps);
1205  CHECK(std::abs(dGmass_dT_analyt / dGmass_dT_num - 1) < eps);
1206  CHECK(std::abs(dGmass_drho_analyt / dGmass_drho_num - 1) < eps);
1207 
1208  CHECK(std::abs(dCvmolar_dT_analyt / dCvmolar_dT_num - 1) < eps);
1209  CHECK(std::abs(dCvmolar_drho_analyt / dCvmolar_drho_num - 1) < eps);
1210  CHECK(std::abs(dCvmass_dT_analyt / dCvmass_dT_num - 1) < eps);
1211  CHECK(std::abs(dCvmass_drho_analyt / dCvmass_drho_num - 1) < eps);
1212 
1213  CHECK(std::abs(dCpmolar_dT_analyt / dCpmolar_dT_num - 1) < eps);
1214  CHECK(std::abs(dCpmolar_drho_analyt / dCpmolar_drho_num - 1) < eps);
1215  CHECK(std::abs(dCpmass_dT_analyt / dCpmass_dT_num - 1) < eps);
1216  CHECK(std::abs(dCpmass_drho_analyt / dCpmass_drho_num - 1) < eps);
1217 
1218  CHECK(std::abs(dspeed_sound_dT_analyt / dspeed_sound_dT_num - 1) < eps);
1219  CHECK(std::abs(dspeed_sound_drho_analyt / dspeed_sound_drho_num - 1) < eps);
1220 }
1221 
1222 #endif