CoolProp  4.2.5
An open-source fluid property and humid air property database
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TTSE.cpp
Go to the documentation of this file.
1 #if defined(_MSC_VER)
2 #define _CRT_SECURE_NO_WARNINGS
3 #endif
4 
5 #include "CoolPropTools.h"
6 
7 #if defined(__ISWINDOWS__)
8 #include <windows.h> // for the CreateDirectory function
9 #else
10  #if !defined(__powerpc__)
11  #include <pwd.h>
12  #endif
13 #include <unistd.h>
14 #include <sys/types.h>
15 #include <sys/stat.h>
16 #endif
17 
18 
19 
20 #include "CoolProp.h"
21 #include "CPState.h"
22 #include "TTSE.h"
23 #include <vector>
24 #include <stdlib.h>
25 #include <stdio.h>
26 #include <iostream>
27 #include <fstream>
28 #include <streambuf>
29 #include <time.h>
30 #include <sys/types.h>
31 #include <sys/stat.h>
32 #include <errno.h>
33 #include <cerrno>
34 #include <float.h>
35 
36 // An ugly hack to disable the timing function on PPC since the sysClkRate function not found
37 #if defined(__powerpc__)
38 #define CLOCKS_PER_SEC 1000
39 #endif
40 
41 // The revision of the TTSE tables, only use tables with the same revision. Increment this macro if any non-forward compatible changes are made
42 #define TTSEREV 6
43 
44 std::string get_home_dir(void)
45 {
46 
47  // See http://stackoverflow.com/questions/2552416/how-can-i-find-the-users-home-dir-in-a-cross-platform-manner-using-c
48  #if defined(__ISLINUX__)
49  char *home = NULL;
50  home = getenv("HOME");
51  return std::string(home);
52  #elif defined(__ISAPPLE__)
53  char *home = NULL;
54  home = getenv("HOME");
55  if (home==NULL) {
56  struct passwd* pwd = getpwuid(getuid());
57  if (pwd) {
58  home = pwd->pw_dir;
59  }
60  }
61  if (home==NULL) {
62  throw NotImplementedError("Could not detect home directory.");
63  }
64 // throw NotImplementedError("This function is not defined for your platform.");
65  return std::string(home);
66  #elif defined(__ISWINDOWS__)
67  char * pUSERPROFILE = getenv("USERPROFILE");
68  if (pUSERPROFILE != NULL) {
69  return std::string(pUSERPROFILE);
70  } else {
71  char * pHOMEDRIVE = getenv("HOMEDRIVE");
72  char * pHOMEPATH = getenv("HOMEPATH");
73  if (pHOMEDRIVE != NULL && pHOMEPATH != NULL) {
74  return std::string(pHOMEDRIVE) + std::string(pHOMEPATH);
75  } else {
76  return std::string("");
77  }
78  }
79  #else
80  throw NotImplementedError("This function is not defined for your platform.");
81  #endif
82 }
83 
85 const static double Ainv[16][16] = {
86  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
87  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
88  {-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
89  { 2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
90  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
91  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
92  { 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0},
93  { 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0},
94  {-3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0},
95  { 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0},
96  { 9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1},
97  {-6, 6, 6, -6, -3, -3, 3, 3, -4, 4, -2, 2, -2, -2, -1, -1},
98  { 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0},
99  { 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0},
100  {-6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1},
101  { 4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1}
102  };
103 
104 double round(double r) {
105  return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
106 }
107 
108 void matrix_vector_product(std::vector<double> *x, std::vector<double> *y)
109 {
110  double sum;
111  for (unsigned int i = 0; i < 16; i++)
112  {
113  sum = 0;
114  for (unsigned int j = 0; j < 16; j++)
115  {
116  sum += Ainv[i][j]*(*x)[j];
117  }
118  (*y)[i] = sum;
119  }
120 }
121 
122 TTSESinglePhaseTableClass::TTSESinglePhaseTableClass(){
123  this->enable_writing_tables_to_files = true;
124  this->enable_transport = true;
125  SatL = NULL;
126  SatV = NULL;
127 
128  // Default to use the "normal" TTSE evaluation
129  mode = TTSE_MODE_TTSE;
130 }
131 
132 TTSESinglePhaseTableClass::TTSESinglePhaseTableClass(Fluid *pFluid) {
133  this->pFluid = pFluid;
134  // The default data location for the LUT
135  #if defined(__ISWINDOWS__)
136  this->root_path = get_home_dir()+std::string("\\.CoolProp-TTSEData\\")+pFluid->get_name();
137  #else
138  this->root_path = get_home_dir()+std::string("/.CoolProp-TTSEData/")+pFluid->get_name();
139  #endif
140 
141  // Seed the generator for random number generation
142  srand((unsigned int)time(NULL));
143  this->enable_writing_tables_to_files = true;
144  this->enable_transport = true;
145  SatL = NULL;
146  SatV = NULL;
147 
148  // Instantiate the storage vectors for bicubic interpolation; instantiated once since instantation is quite slow
149  alpha_bicubic = std::vector<double>(16);
150  z_bicubic = std::vector<double>(16);
151 
152  // Default to use the "normal" TTSE evaluation
153  mode = TTSE_MODE_TTSE;
154 }
155 
156 int TTSESinglePhaseTableClass::get_mode() {
157  return this->mode;
158 }
159 
160 int TTSESinglePhaseTableClass::set_mode(int mode) {
161  if (mode == TTSE_MODE_TTSE || mode == TTSE_MODE_BICUBIC) {
162  this->mode = mode;
163  return true;
164  } else {
165  return false;
166  }
167 }
168 
169 void TTSESinglePhaseTableClass::set_size_ph(unsigned int Np, unsigned int Nh) {
170  this->Nh = Nh;
171  this->Np = Np;
172 
173  h.resize(Nh);
174  p.resize(Np);
175 
176  s.resize(Nh, std::vector<double>(Np, _HUGE));
177  dsdh.resize(Nh, std::vector<double>(Np, _HUGE));
178  dsdp.resize(Nh, std::vector<double>(Np, _HUGE));
179  d2sdh2.resize(Nh, std::vector<double>(Np, _HUGE));
180  d2sdp2.resize(Nh, std::vector<double>(Np, _HUGE));
181  d2sdhdp.resize(Nh, std::vector<double>(Np, _HUGE));
182 
183  T.resize(Nh, std::vector<double>(Np, _HUGE));
184  dTdh.resize(Nh, std::vector<double>(Np, _HUGE));
185  dTdp.resize(Nh, std::vector<double>(Np, _HUGE));
186  d2Tdh2.resize(Nh, std::vector<double>(Np, _HUGE));
187  d2Tdp2.resize(Nh, std::vector<double>(Np, _HUGE));
188  d2Tdhdp.resize(Nh, std::vector<double>(Np, _HUGE));
189 
190  rho.resize(Nh, std::vector<double>(Np, _HUGE));
191  drhodh.resize(Nh, std::vector<double>(Np, _HUGE));
192  drhodp.resize(Nh, std::vector<double>(Np, _HUGE));
193  d2rhodh2.resize(Nh, std::vector<double>(Np, _HUGE));
194  d2rhodp2.resize(Nh, std::vector<double>(Np, _HUGE));
195  d2rhodhdp.resize(Nh, std::vector<double>(Np, _HUGE));
196 
197  IL.resize(Np);
198  IV.resize(Np);
199  TL.resize(Np);
200  TV.resize(Np);
201  SL.resize(Np);
202  SV.resize(Np);
203  DL.resize(Np);
204  DV.resize(Np);
205 
206  // Instantiate the cell matrices for the bicubic interpolation
207  bicubic_cells = BiCubicCellsContainerClass();
208 
209  for(unsigned int i = 0; i < Nh; i++){
210  bicubic_cells.cells.push_back(std::vector<BiCubicCellClass>(Np));
211  }
212 }
213 
214 void TTSESinglePhaseTableClass::set_size_Trho(unsigned int NT, unsigned int Nrho) {
215  this->NT = NT;
216  this->Nrho = Nrho;
217 
218  T_Trho.resize(NT);
219  rho_Trho.resize(Nrho);
220 
221  s_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
222  dsdT_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
223  dsdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
224  d2sdT2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
225  d2sdrho2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
226  d2sdTdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
227 
228  p_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
229  dpdT_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
230  dpdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
231  d2pdT2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
232  d2pdrho2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
233  d2pdTdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
234 
235  h_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
236  dhdT_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
237  dhdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
238  d2hdT2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
239  d2hdrho2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
240  d2hdTdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
241 
242  k_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
243  dkdT_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
244  dkdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
245  d2kdT2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
246  d2kdrho2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
247  d2kdTdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
248 
249  mu_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
250  dmudT_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
251  dmudrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
252  d2mudT2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
253  d2mudrho2_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
254  d2mudTdrho_Trho.resize(NT, std::vector<double>(Nrho, _HUGE));
255 }
256 
257 void TTSESinglePhaseTableClass::nearest_good_neighbor(int *i, int *j) {
258  // Left
259  if (*i>0 && ValidNumber(rho[*i-1][*j]) && ValidNumber(T[*i-1][*j])){
260  *i -= 1;
261  return;
262  }
263  // Right
264  else if (*i<(int)Nh-1 && ValidNumber(rho[*i+1][*j]) && ValidNumber(T[*i+1][*j])){
265  *i += 1;
266  return;
267  }
268  // Down
269  else if (*j>0 && ValidNumber(rho[*i][*j-1]) && ValidNumber(T[*i][*j-1])){
270  *j -= 1;
271  return;
272  }
273  // Up
274  else if (*j<(int)Np-1 && ValidNumber(rho[*i][*j+1]) && ValidNumber(T[*i][*j+1])){
275  *j += 1;
276  return;
277  }
278  // Two Left
279  if (*i>1 && ValidNumber(rho[*i-2][*j]) && ValidNumber(T[*i-2][*j])){
280  *i -= 2;
281  return;
282  }
283  // Two Right
284  else if (*i<(int)Nh-2 && ValidNumber(rho[*i+2][*j]) && ValidNumber(T[*i+2][*j])){
285  *i += 2;
286  return;
287  }
288  // Two Down
289  else if (*j>1 && ValidNumber(rho[*i][*j-2]) && ValidNumber(T[*i][*j-2])){
290  *j -= 2;
291  return;
292  }
293  // Two Up
294  else if (*j<(int)Np-2 && ValidNumber(rho[*i][*j+2]) && ValidNumber(T[*i][*j+2])){
295  *j += 2;
296  return;
297  }
298 
299  else
300  {
301  throw ValueError(format("No neighbors found for %d,%d",i,j));
302  return;
303  }
304 }
305 
306 void TTSESinglePhaseTableClass::nearest_good_neighbor_Trho_interpolate(int *i, int *j)
307 {
308  if (bicubic_cells.cells[*i][*j].valid_Trho){
309  return;
310  }
311  else if (bicubic_cells.cells[*i+1][*j].valid_Trho)
312  {
313  *i += 1; return;
314  }
315  else if (bicubic_cells.cells[*i-1][*j].valid_Trho)
316  {
317  *i -= 1; return;
318  }
319  else if (bicubic_cells.cells[*i][*j-1].valid_Trho)
320  {
321  *j -= 1; return;
322  }
323  else if (bicubic_cells.cells[*i][*j+1].valid_Trho)
324  {
325  *j += 1; return;
326  }
327 
328  else if (bicubic_cells.cells[*i+1][*j+1].valid_Trho)
329  {
330  *i += 1; *j += 1; return;
331  }
332  else if (bicubic_cells.cells[*i-1][*j-1].valid_Trho)
333  {
334  *i -= 1; *j -= 1; return;
335  }
336  else if (bicubic_cells.cells[*i+1][*j-1].valid_Trho)
337  {
338  *j -= 1; *i += 1; return;
339  }
340  else if (bicubic_cells.cells[*i-1][*j+1].valid_Trho)
341  {
342  *j += 1; *i -= 1; return;
343  }
344  else
345  {
346  throw ValueError(format("No neighbors found for %d,%d",*i,*j));
347  return;
348  }
349 }
350 
351 void TTSESinglePhaseTableClass::nearest_good_neighbor_ph_interpolate(int *i, int *j)
352 {
353  if (bicubic_cells.cells[*i][*j].valid_hp){
354  return;
355  }
356  else if (bicubic_cells.cells[*i+1][*j].valid_hp)
357  {
358  *i += 1; return;
359  }
360  else if (bicubic_cells.cells[*i-1][*j].valid_hp)
361  {
362  *i -= 1; return;
363  }
364  else if (bicubic_cells.cells[*i][*j-1].valid_hp)
365  {
366  *j -= 1; return;
367  }
368  else if (bicubic_cells.cells[*i][*j+1].valid_hp)
369  {
370  *j += 1; return;
371  }
372 
373  else if (bicubic_cells.cells[*i+1][*j+1].valid_hp)
374  {
375  *i += 1; *j += 1; return;
376  }
377  else if (bicubic_cells.cells[*i-1][*j-1].valid_hp)
378  {
379  *i -= 1; *j -= 1; return;
380  }
381  else if (bicubic_cells.cells[*i+1][*j-1].valid_hp)
382  {
383  *j -= 1; *i += 1; return;
384  }
385  else if (bicubic_cells.cells[*i-1][*j+1].valid_hp)
386  {
387  *j += 1; *i -= 1; return;
388  }
389  else
390  {
391  throw ValueError(format("No neighbors found for %d,%d",*i,*j));
392  return;
393  }
394 }
395 
396 void TTSESinglePhaseTableClass::nearest_good_neighbor_Trho(int *i, int *j)
397 {
398  // Left
399  if (*i>0 && ValidNumber(h_Trho[*i-1][*j]) && ValidNumber(p_Trho[*i-1][*j])){
400  *i -= 1;
401  return;
402  }
403  // Right
404  else if (*i<(int)Nh-1 && ValidNumber(h_Trho[*i+1][*j]) && ValidNumber(p_Trho[*i+1][*j])){
405  *i += 1;
406  return;
407  }
408  // Down
409  else if (*j>0 && ValidNumber(h_Trho[*i][*j-1]) && ValidNumber(p_Trho[*i][*j-1])){
410  *j -= 1;
411  return;
412  }
413  // Up
414  else if (*j<(int)Np-1 && ValidNumber(h_Trho[*i][*j+1]) && ValidNumber(p_Trho[*i][*j+1])){
415  *j += 1;
416  return;
417  }
418  else
419  {
420  throw ValueError(format("No neighbors found for %d,%d",i,j));
421  return;
422  }
423 }
424 
425 void TTSESinglePhaseTableClass::nearest_neighbor_ph(int i, int j, double *T0, double *rho0)
426 {
427  // Left
428  if (i>0 && ValidNumber(rho[i-1][j]) && ValidNumber(T[i-1][j])){
429  *T0 = T[i-1][j];
430  *rho0 = rho[i-1][j];
431  return;
432  }
433  // Right
434  else if (i<(int)Nh-1 && ValidNumber(rho[i+1][j]) && ValidNumber(T[i+1][j])){
435  *T0 = T[i+1][j];
436  *rho0 = rho[i+1][j];
437  return;
438  }
439  // Down
440  else if (j>0 && ValidNumber(rho[i][j-1]) && ValidNumber(T[i][j-1])){
441  *T0 = T[i][j-1];
442  *rho0 = rho[i][j-1];
443  return;
444  }
445  // Up
446  else if (j<(int)Np-1 && ValidNumber(rho[i][j+1]) && ValidNumber(T[i][j+1])){
447  *T0 = T[i][j+1];
448  *rho0 = rho[i][j+1];
449  return;
450  }
451  else
452  {
453  *T0 = -1;
454  *rho0 = -1;
455  return;
456  }
457 }
458 std::string join(std::vector<std::string> strings, char delim)
459 {
460  std::string output = strings[0];
461  for (unsigned int i = 1; i < strings.size(); i++)
462  {
463  output += format("%c%s",delim,strings[i].c_str());
464  }
465  return output;
466 }
467 
468 void TTSESinglePhaseTableClass::matrix_to_file(std::string fName, std::vector< std::vector<double> > *A)
469 {
470  FILE *pFile;
471  pFile = fopen(fName.c_str(),"wb");
472  if (pFile != NULL)
473  {
474  for (unsigned int i = 0; i < Nh; i++)
475  {
476  fwrite(&(((*A)[i])[0]), sizeof(double),Np,pFile);
477  }
478  fclose(pFile);
479  }
480 }
481 
482 void TTSESinglePhaseTableClass::matrix_from_file(std::string fName, std::vector<std::vector<double> > *A)
483 {
484  FILE *pFile;
485  pFile = fopen(fName.c_str(),"rb");
486  if (pFile != NULL)
487  {
488  for (unsigned int i = 0; i < Nh; i++)
489  {
490  /*std::vector<double> row(Nh,0);
491  fread(&(row[0]), sizeof(double), Np, pFile);
492  double tgregt = 1;*/
493  fread(&(((*A)[i])[0]), sizeof(double),Np,pFile);
494  }
495  fclose(pFile);
496  }
497 }
498 
499 void TTSESinglePhaseTableClass::vector_to_file(std::string fName, std::vector<double>*A)
500 {
501  FILE *pFile;
502  pFile = fopen(fName.c_str(),"wb");
503  if (pFile != NULL)
504  {
505  fwrite((const char *)&(*A).front(), sizeof(double),(*A).size(),pFile);
506  fclose(pFile);
507  }
508 }
509 void TTSESinglePhaseTableClass::vector_from_file(std::string fName, int N, std::vector<double> *vec)
510 {
511  FILE *pFile;
512  pFile = fopen(fName.c_str(),"rb");
513  if (pFile != NULL)
514  {
515  fread(&((*vec)[0]), sizeof(double), N, pFile);
516  fclose(pFile);
517  }
518 }
519 bool pathExists(std::string path)
520 {
521  #if defined(__ISWINDOWS__) // Defined for 32-bit and 64-bit windows
522  struct _stat buf;
523  // Get data associated with path using the windows libraries,
524  // and if you can (result == 0), the path exists
525  if ( _stat( path.c_str(), &buf) == 0)
526  return true;
527  else
528  return false;
529  #elif defined(__ISLINUX__) || defined(__ISAPPLE__)
530  struct stat st;
531 // if(stat(path.c_str(),&st) == 0) {
532 // if(st.st_mode & S_IFDIR != 0) return true;
533 // if(st.st_mode & S_IFREG != 0) return true;
534 // return false;
535 // } else {
536 // return false;
537 // }
538  if(lstat(path.c_str(),&st) == 0) {
539  if(S_ISDIR(st.st_mode)) return true;
540  if(S_ISREG(st.st_mode)) return true;
541  return false;
542  } else {
543  return false;
544  }
545 
546  #else
547  throw NotImplementedError("This function is not defined for your platform.");
548  #endif
549 }
550 //bool fileExists(const char *fileName)
551 //{
552 // std::ifstream infile(fileName);
553 // return infile.good();
554 //}
555 void make_dirs(std::string file_path)
556 {
557  std::vector<std::string> pathsplit = strsplit(file_path,'/');
558  std::string path = pathsplit[0];
559  if (pathsplit.size()>0)
560  {
561  for (unsigned int i = 0; i < pathsplit.size(); i++)
562  {
563  if (!pathExists(path))
564  {
565  #ifdef _WIN32
566  #if defined(_UNICODE)
567  CreateDirectoryA((LPCSTR)path.c_str(),NULL);
568  #else
569  CreateDirectory((LPCSTR)path.c_str(),NULL);
570  #endif
571  #else
572  #if defined(__powerpc__)
573  #else
574  mkdir(path.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
575  #endif
576  #endif
577  }
578  if (i < pathsplit.size()-1)
579  path += format("/%s",pathsplit[i+1].c_str());
580  }
581  }
582  else
583  {
584  throw ValueError(format("Could not make path [%s]",file_path.c_str()));
585  }
586 }
587 
588 bool TTSESinglePhaseTableClass::read_all_from_file(std::string root_path)
589 {
590  std::string Fluid;
591  double hmin,hmax,pmin,pmax;
592  int Np, Nh, TTSERev;
593 
594  // Replace any '\' with '/' in the path
595  for (unsigned int i = 0; i<root_path.length(); i++)
596  {
597  if (root_path[i] == '\\') root_path[i] = '/';
598  }
599  // If it ends with a '/', remove it for now
600  if (root_path[root_path.size()-1] == '/')
601  root_path = std::string(root_path,0,root_path.size()-1);
602 
603  if (!pathExists(root_path)) return false;
604 
605  // Append a '/'
606  root_path += format("/");
607 
608  if (!pathExists(root_path+format("Info_ph.txt"))) return false;
609 
610  // Parse the info file to find the dimensions and boundaries and check that they are the same
611  std::string info = get_file_contents((root_path+format("Info_ph.txt")).c_str());
612 
613  // Split into lines
614 #if defined(__ISWINDOWS__)
615  std::vector<std::string> lines = strsplit(info,'\r');
616 #else
617  std::vector<std::string> lines = strsplit(info,'\n');
618 #endif
619 
620  for (unsigned int i = 0; i< lines.size(); i++) {
621 
622  // Split at ':'
623  std::vector<std::string> line = strsplit(lines[i],':');
624  if (line.size() == 1) // No : found
625  continue;
626 
627  if (line[0].find("Fluid")!= std::string::npos){
628  Fluid = line[1];}
629  else if (line[0].find("hmin")!= std::string::npos) {
630  hmin = strtod(line[1].c_str(),NULL);}
631  else if (line[0].find("hmax")!= std::string::npos) {
632  hmax = strtod(line[1].c_str(),NULL);}
633  else if (line[0].find("pmin")!= std::string::npos) {
634  pmin = strtod(line[1].c_str(),NULL);}
635  else if (line[0].find("pmax")!= std::string::npos) {
636  pmax = strtod(line[1].c_str(),NULL);}
637  else if (line[0].find("Np")!= std::string::npos) {
638  Np = (int)strtol(line[1].c_str(),NULL,0);}
639  else if (line[0].find("Nh")!= std::string::npos) {
640  Nh = (int)strtol(line[1].c_str(),NULL,0);}
641  else if (line[0].find("Tmin")!= std::string::npos) {
642  Tmin = strtod(line[1].c_str(),NULL);}
643  else if (line[0].find("Tmax")!= std::string::npos) {
644  Tmax = strtod(line[1].c_str(),NULL);}
645  else if (line[0].find("rhomin")!= std::string::npos) {
646  rhomin = strtod(line[1].c_str(),NULL);}
647  else if (line[0].find("rhomax")!= std::string::npos) {
648  rhomax = strtod(line[1].c_str(),NULL);}
649  else if (line[0].find("NT")!= std::string::npos) {
650  NT = (int)strtol(line[1].c_str(),NULL,0);}
651  else if (line[0].find("Nrho")!= std::string::npos) {
652  Nrho = (int)strtol(line[1].c_str(),NULL,0);}
653  else if (line[0].find("TTSERev")!= std::string::npos) {
654  TTSERev = (int)strtol(line[1].c_str(),NULL,0);}
655  }
656 
657  // Didn't work since at least one of the parameters was different
658  // so we need to build the tables again
659  if (!(!Fluid.compare(pFluid->get_name())
660  && Nh == (int)this->Nh
661  && Np == (int)this->Np
662  && fabs(pmin - this->pmin)<10*DBL_EPSILON
663  && fabs(pmax - this->pmax)<10*DBL_EPSILON
664  && fabs(hmin - this->hmin)<10*DBL_EPSILON
665  && fabs(hmax - this->hmax)<10*DBL_EPSILON
666  && NT == this->NT
667  && Nrho == this->Nrho
668  && fabs(Tmin - this->Tmin)<10*DBL_EPSILON
669  && fabs(Tmax - this->Tmax)<10*DBL_EPSILON
670  && fabs(rhomin - this->rhomin)<10*DBL_EPSILON
671  && fabs(rhomax - this->rhomax)<10*DBL_EPSILON
672  && TTSERev == TTSEREV
673  )) return false;
674 
675  // Read all the data from the binary files
676  vector_from_file(root_path + std::string("p_ph.ttse"),Np,&p);
677  vector_from_file(root_path + std::string("h_ph.ttse"),Nh,&h);
678  matrix_from_file(root_path + std::string("T_ph.ttse"),&T);
679  matrix_from_file(root_path + std::string("dTdh_ph.ttse"),&dTdh);
680  matrix_from_file(root_path + std::string("dTdp_ph.ttse"),&dTdp);
681  matrix_from_file(root_path + std::string("d2Tdh2_ph.ttse"),&d2Tdh2);
682  matrix_from_file(root_path + std::string("d2Tdp2_ph.ttse"),&d2Tdp2);
683  matrix_from_file(root_path + std::string("d2Tdhdp_ph.ttse"),&d2Tdhdp);
684  matrix_from_file(root_path + std::string("s_ph.ttse"),&s);
685  matrix_from_file(root_path + std::string("dsdh_ph.ttse"),&dsdh);
686  matrix_from_file(root_path + std::string("dsdp_ph.ttse"),&dsdp);
687  matrix_from_file(root_path + std::string("d2sdh2_ph.ttse"),&d2sdh2);
688  matrix_from_file(root_path + std::string("d2sdp2_ph.ttse"),&d2sdp2);
689  matrix_from_file(root_path + std::string("d2sdhdp_ph.ttse"),&d2sdhdp);
690  matrix_from_file(root_path + std::string("rho_ph.ttse"),&rho);
691  matrix_from_file(root_path + std::string("drhodh_ph.ttse"),&drhodh);
692  matrix_from_file(root_path + std::string("drhodp_ph.ttse"),&drhodp);
693  matrix_from_file(root_path + std::string("d2rhodh2_ph.ttse"),&d2rhodh2);
694  matrix_from_file(root_path + std::string("d2rhodp2_ph.ttse"),&d2rhodp2);
695  matrix_from_file(root_path + std::string("d2rhoTdhdp_ph.ttse"),&d2rhodhdp);
696 
697  vector_from_file(root_path + std::string("T_Trho.ttse"),NT,&T_Trho);
698  vector_from_file(root_path + std::string("rho_Trho.ttse"),Nrho,&rho_Trho);
699  matrix_from_file(root_path + std::string("p_Trho.ttse"),&p_Trho);
700  matrix_from_file(root_path + std::string("dpdT_Trho.ttse"),&dpdT_Trho);
701  matrix_from_file(root_path + std::string("dpdrho_Trho.ttse"),&dpdrho_Trho);
702  matrix_from_file(root_path + std::string("d2pdT2_Trho.ttse"),&d2pdT2_Trho);
703  matrix_from_file(root_path + std::string("d2pdrho2_Trho.ttse"),&d2pdrho2_Trho);
704  matrix_from_file(root_path + std::string("d2pdTdrho_Trho.ttse"),&d2pdTdrho_Trho);
705  matrix_from_file(root_path + std::string("s_Trho.ttse"),&s_Trho);
706  matrix_from_file(root_path + std::string("dsdT_Trho.ttse"),&dsdT_Trho);
707  matrix_from_file(root_path + std::string("dsdrho_Trho.ttse"),&dsdrho_Trho);
708  matrix_from_file(root_path + std::string("d2sdT2_Trho.ttse"),&d2sdT2_Trho);
709  matrix_from_file(root_path + std::string("d2sdrho2_Trho.ttse"),&d2sdrho2_Trho);
710  matrix_from_file(root_path + std::string("d2sdTdrho_Trho.ttse"),&d2sdTdrho_Trho);
711  matrix_from_file(root_path + std::string("h_Trho.ttse"),&h_Trho);
712  matrix_from_file(root_path + std::string("dhdT_Trho.ttse"),&dhdT_Trho);
713  matrix_from_file(root_path + std::string("dhdrho_Trho.ttse"),&dhdrho_Trho);
714  matrix_from_file(root_path + std::string("d2hdT2_Trho.ttse"),&d2hdT2_Trho);
715  matrix_from_file(root_path + std::string("d2hdrho2_Trho.ttse"),&d2hdrho2_Trho);
716  matrix_from_file(root_path + std::string("d2hdTdrho_Trho.ttse"),&d2hdTdrho_Trho);
717  matrix_from_file(root_path + std::string("mu_Trho.ttse"),&mu_Trho);
718  matrix_from_file(root_path + std::string("dmudT_Trho.ttse"),&dmudT_Trho);
719  matrix_from_file(root_path + std::string("dmudrho_Trho.ttse"),&dmudrho_Trho);
720  matrix_from_file(root_path + std::string("d2mudT2_Trho.ttse"),&d2mudT2_Trho);
721  matrix_from_file(root_path + std::string("d2mudrho2_Trho.ttse"),&d2mudrho2_Trho);
722  matrix_from_file(root_path + std::string("d2mudTdrho_Trho.ttse"),&d2mudTdrho_Trho);
723  matrix_from_file(root_path + std::string("k_Trho.ttse"),&k_Trho);
724  matrix_from_file(root_path + std::string("dkdT_Trho.ttse"),&dkdT_Trho);
725  matrix_from_file(root_path + std::string("dkdrho_Trho.ttse"),&dkdrho_Trho);
726  matrix_from_file(root_path + std::string("d2kdT2_Trho.ttse"),&d2kdT2_Trho);
727  matrix_from_file(root_path + std::string("d2kdrho2_Trho.ttse"),&d2kdrho2_Trho);
728  matrix_from_file(root_path + std::string("d2kdTdrho_Trho.ttse"),&d2kdTdrho_Trho);
729 
730  this->pratio = pow(pmax/pmin,1/((double)Np-1));
731  this->logpratio = log(pratio); // For speed since log() is a slow function
732  this->logpmin = log(pmin);
733  this->rhoratio = pow(rhomax/rhomin,1/((double)Nrho-1));
734  this->logrhoratio = log(rhoratio); // For speed since log() is a slow function
735  this->logrhomin = log(rhomin);
736  this->jpcrit_floor = (int)floor((log(pFluid->reduce.p.Pa)-logpmin)/logpratio);
737  this->jpcrit_ceil = (int)ceil((log(pFluid->reduce.p.Pa)-logpmin)/logpratio);
738 
739  update_saturation_boundary_indices();
740 
741  update_cell_validity();
742 
743  return true;
744 }
745 void TTSESinglePhaseTableClass::write_all_to_file(std::string root_path)
746 {
747  // Replace any '\' with '/' in the path
748  for (unsigned int i = 0; i<root_path.length(); i++)
749  {
750  if (root_path[i] == '\\') root_path[i] = '/';
751  }
752  // If it ends with a '/', remove it for now
753  if (root_path[root_path.size()-1] == '/')
754  root_path = std::string(root_path,0,root_path.size()-1);
755 
756  if (!pathExists(root_path))
757  make_dirs(root_path);
758 
759  // Append a '/'
760  root_path += format("/");
761 
762  std::string header = std::string("Data for the TTSE method\nDO NOT CHANGE ANY OF THESE PARAMETERS FOR ANY REASON!\n\n");
763 
764  header += format("TTSERev:%d\nFluid:%s\npmin:%23.19g\npmax:%23.19g\nNp:%25d\nhmin:%23.19g\nhmax:%23.19g\nNh:%25d\nTmin:%23.19g\nTmax:%23.19g\nNT:%25d\nrhomin:%23.19g\nrhomax:%23.19g\nNrho:%25d\n",TTSEREV,pFluid->get_name().c_str(),pmin,pmax,Np,hmin,hmax,Nh,Tmin,Tmax,NT,rhomin,rhomax,Nrho);
765 
766  clock_t t1,t2;
767  t1 = clock();
768 
769  // Write the header information to a text file
770  FILE *fp;
771  fp = fopen((root_path+std::string("Info_ph.txt")).c_str(),"w");
772  fprintf(fp,"%s",header.c_str());
773  fclose(fp);
774 
775  // Write each of these files in binary mode
776  vector_to_file(root_path + std::string("p_ph.ttse"),&p);
777  vector_to_file(root_path + std::string("h_ph.ttse"),&h);
778  matrix_to_file(root_path + std::string("T_ph.ttse"),&T);
779  matrix_to_file(root_path + std::string("dTdh_ph.ttse"),&dTdh);
780  matrix_to_file(root_path + std::string("dTdp_ph.ttse"),&dTdp);
781  matrix_to_file(root_path + std::string("d2Tdh2_ph.ttse"),&d2Tdh2);
782  matrix_to_file(root_path + std::string("d2Tdp2_ph.ttse"),&d2Tdp2);
783  matrix_to_file(root_path + std::string("d2Tdhdp_ph.ttse"),&d2Tdhdp);
784  matrix_to_file(root_path + std::string("s_ph.ttse"),&s);
785  matrix_to_file(root_path + std::string("dsdh_ph.ttse"),&dsdh);
786  matrix_to_file(root_path + std::string("dsdp_ph.ttse"),&dsdp);
787  matrix_to_file(root_path + std::string("d2sdh2_ph.ttse"),&d2sdh2);
788  matrix_to_file(root_path + std::string("d2sdp2_ph.ttse"),&d2sdp2);
789  matrix_to_file(root_path + std::string("d2sdhdp_ph.ttse"),&d2sdhdp);
790  matrix_to_file(root_path + std::string("rho_ph.ttse"),&rho);
791  matrix_to_file(root_path + std::string("drhodh_ph.ttse"),&drhodh);
792  matrix_to_file(root_path + std::string("drhodp_ph.ttse"),&drhodp);
793  matrix_to_file(root_path + std::string("d2rhodh2_ph.ttse"),&d2rhodh2);
794  matrix_to_file(root_path + std::string("d2rhodp2_ph.ttse"),&d2rhodp2);
795  matrix_to_file(root_path + std::string("d2rhoTdhdp_ph.ttse"),&d2rhodhdp);
796 
797  vector_to_file(root_path + std::string("T_Trho.ttse"),&T_Trho);
798  vector_to_file(root_path + std::string("rho_Trho.ttse"),&rho_Trho);
799  matrix_to_file(root_path + std::string("p_Trho.ttse"),&p_Trho);
800  matrix_to_file(root_path + std::string("dpdT_Trho.ttse"),&dpdT_Trho);
801  matrix_to_file(root_path + std::string("dpdrho_Trho.ttse"),&dpdrho_Trho);
802  matrix_to_file(root_path + std::string("d2pdT2_Trho.ttse"),&d2pdT2_Trho);
803  matrix_to_file(root_path + std::string("d2pdrho2_Trho.ttse"),&d2pdrho2_Trho);
804  matrix_to_file(root_path + std::string("d2pdTdrho_Trho.ttse"),&d2pdTdrho_Trho);
805  matrix_to_file(root_path + std::string("s_Trho.ttse"),&s_Trho);
806  matrix_to_file(root_path + std::string("dsdT_Trho.ttse"),&dsdT_Trho);
807  matrix_to_file(root_path + std::string("dsdrho_Trho.ttse"),&dsdrho_Trho);
808  matrix_to_file(root_path + std::string("d2sdT2_Trho.ttse"),&d2sdT2_Trho);
809  matrix_to_file(root_path + std::string("d2sdrho2_Trho.ttse"),&d2sdrho2_Trho);
810  matrix_to_file(root_path + std::string("d2sdTdrho_Trho.ttse"),&d2sdTdrho_Trho);
811  matrix_to_file(root_path + std::string("h_Trho.ttse"),&h_Trho);
812  matrix_to_file(root_path + std::string("dhdT_Trho.ttse"),&dhdT_Trho);
813  matrix_to_file(root_path + std::string("dhdrho_Trho.ttse"),&dhdrho_Trho);
814  matrix_to_file(root_path + std::string("d2hdT2_Trho.ttse"),&d2hdT2_Trho);
815  matrix_to_file(root_path + std::string("d2hdrho2_Trho.ttse"),&d2hdrho2_Trho);
816  matrix_to_file(root_path + std::string("d2hdTdrho_Trho.ttse"),&d2hdTdrho_Trho);
817  matrix_to_file(root_path + std::string("k_Trho.ttse"),&k_Trho);
818  matrix_to_file(root_path + std::string("dkdT_Trho.ttse"),&dkdT_Trho);
819  matrix_to_file(root_path + std::string("dkdrho_Trho.ttse"),&dkdrho_Trho);
820  matrix_to_file(root_path + std::string("d2kdT2_Trho.ttse"),&d2kdT2_Trho);
821  matrix_to_file(root_path + std::string("d2kdrho2_Trho.ttse"),&d2kdrho2_Trho);
822  matrix_to_file(root_path + std::string("d2kdTdrho_Trho.ttse"),&d2kdTdrho_Trho);
823  matrix_to_file(root_path + std::string("mu_Trho.ttse"),&mu_Trho);
824  matrix_to_file(root_path + std::string("dmudT_Trho.ttse"),&dmudT_Trho);
825  matrix_to_file(root_path + std::string("dmudrho_Trho.ttse"),&dmudrho_Trho);
826  matrix_to_file(root_path + std::string("d2mudT2_Trho.ttse"),&d2mudT2_Trho);
827  matrix_to_file(root_path + std::string("d2mudrho2_Trho.ttse"),&d2mudrho2_Trho);
828  matrix_to_file(root_path + std::string("d2mudTdrho_Trho.ttse"),&d2mudTdrho_Trho);
829 
830  t2 = clock();
831  std::cout << "write time: " << (double)(t2-t1)/CLOCKS_PER_SEC << std::endl;
832 }
833 
834 double TTSESinglePhaseTableClass::build_ph(double hmin, double hmax, double pmin, double pmax, TTSETwoPhaseTableClass *SatL, TTSETwoPhaseTableClass *SatV)
835 {
836  bool SinglePhase = false;
837 
838  this->hmin = hmin;
839  this->hmax = hmax;
840  this->pmin = pmin;
841  this->pmax = pmax;
842  this->logpmin = log(pmin);
843 
844  CoolPropStateClassSI CPS(pFluid);
845 
846  long iFluid = get_Fluid_index(CPS.pFluid->get_name());
847 
848  double dh = (hmax - hmin)/(Nh - 1);
849  pratio = pow(pmax/pmin,1/((double)Np-1));
850  logpratio = log(pratio);
851 
852  clock_t t1,t2;
853  t1 = clock();
854  for (unsigned int i = 0; i<Nh; i++)
855  {
856  double hval = hmin + i*dh;
857  h[i] = hval;
858  for (unsigned int j = 0; j<Np; j++)
859  {
860  double pval = pmin*pow(pratio,(int)j);
861  p[j] = pval;
862 
863  // Check whether the point is single phase
864  // If pressure between ptriple point and pcrit, might be two-phase or single phase, otherwise definitely single phase
865  if (pval <= pFluid->reduce.p.Pa && pval >= pFluid->params.ptriple)
866  {
867  if (SatL == NULL || SatV == NULL){
868  // Not using TTSE method, use saturation (slow...)
869  CPS.update(iP,pval,iQ,0.5);
870  SinglePhase = (hval < CPS.hL() || hval > CPS.hV());
871  }
872  else{
873  // Using the TTSE method, nice and fast
874  SinglePhase = (hval < SatL->evaluate(iH,pval) || hval > SatV->evaluate(iH,pval));
875  }
876  }
877  else
878  {
879  SinglePhase = true;
880  }
881 
882  // If enthalpy is outside the saturation region, continue and do the calculation as a function of p,h
883  if (SinglePhase)
884  {
885  try
886  {
887  double T0=-1,rho0=-1,T,rho,rhoL,rhoV,TsatL,TsatV;
888 
889  // Find a good point around this point that is single-phase if any of its neighbors have been calculated
890  nearest_neighbor_ph(i,j,&T0,&rho0);
891 
892  // If good T,rho was returned, use it as a guess value to calculate T,rho from p,h more quickly
893  if (T0 > 0 && rho0 > 0)
894  {
895  // Get T,rho from p,h using our guess value
896  CPS.pFluid->temperature_ph(pval,hval,T,rho,rhoL,rhoV,TsatL,TsatV,T0,rho0);
897  CPS.flag_SinglePhase = true;
898  CPS.update(iT, T, iD, rho);
899  }
900  else
901  {
902  // Probably here you are close to the saturation boundary since TTSE says it is single phase, but no neighbors to lean on
903  CPS.flag_SinglePhase = true;
904 
905  double hsatLTTSE,hsatVTTSE;
906  if (SatL == NULL || SatV == NULL){
907  // Not using TTSE method, use saturation (slow...)
908  CPS.update(iP,pval,iQ,0.5);
909  hsatLTTSE = CPS.hL();
910  hsatVTTSE = CPS.hV();
911  }
912  else
913  {
914  hsatLTTSE = SatL->evaluate(iH,pval);
915  hsatVTTSE = SatV->evaluate(iH,pval);
916  }
917 
918  if (fabs(hval-hsatLTTSE)<10)
919  {
920  // Close to the saturated liquid boundary
921  T0 = SatL->evaluate(iT,pval);
922  rho0 = SatL->evaluate(iD,pval);
923  CPS.pFluid->temperature_ph(pval,hval,T,rho,rhoL,rhoV,TsatL,TsatV,T0,rho0);
924  CPS.flag_SinglePhase = true;
925  CPS.update(iT, T, iD, rho);
926  }
927  else if (fabs(hval-hsatVTTSE)<10)
928  {
929  // Close to the saturated vapor boundary
930  T0 = SatV->evaluate(iT,pval);
931  rho0 = SatV->evaluate(iD,pval);
932  CPS.pFluid->temperature_ph(pval,hval,T,rho,rhoL,rhoV,TsatL,TsatV,T0,rho0);
933  CPS.flag_SinglePhase = true;
934  CPS.update(iT, T, iD, rho);
935  }
936  else
937  {
938  CPS.update(iP, pval, iH, hval);
939  }
940  //std::cout << format("%d %d %g %g\n",i,j,pval,hval);
941  }
942 
943  T = CPS.T();
944  rho = CPS.rho();
945  double cp = CPS.cp();
946 
947  double A = CPS.dpdT_constrho()*CPS.dhdrho_constT()-CPS.dpdrho_constT()*CPS.dhdT_constrho();
948 
949  this->T[i][j] = T;
950  dTdh[i][j] = 1/cp;
951  dTdp[i][j] = 1/A*CPS.dhdrho_constT();
952  this->rho[i][j] = rho;
953  drhodh[i][j] = 1/A*CPS.dpdT_constrho();
954  drhodp[i][j] = -1/A*CPS.dhdT_constrho();
955  s[i][j] = CPS.s();
956  dsdh[i][j] = 1/T;
957  dsdp[i][j] = -1/(T*rho);
958 
959  // Matrices for second derivatives of entropy as a function of pressure and enthalpy
960  d2sdh2[i][j] = -1/(T*T)*dTdh[i][j];
961  d2sdhdp[i][j] = -1/(T*T)*dTdp[i][j];
962  d2sdp2[i][j] = 1/(T*T*rho)*dTdp[i][j]+1/(T*rho*rho)*drhodp[i][j];
963 
964  // These are common terms needed for a range of terms for T(h,p) as well as rho(h,p)
965  double dAdT_constrho = CPS.d2pdT2_constrho()*CPS.dhdrho_constT()+CPS.dpdT_constrho()*CPS.d2hdrhodT()-CPS.d2pdrhodT()*CPS.dhdT_constrho()-CPS.dpdrho_constT()*CPS.d2hdT2_constrho();
966  double dAdrho_constT = CPS.d2pdrhodT()*CPS.dhdrho_constT()+CPS.dpdT_constrho()*CPS.d2hdrho2_constT()-CPS.d2pdrho2_constT()*CPS.dhdT_constrho()-CPS.dpdrho_constT()*CPS.d2hdrhodT();
967 
968  //Matrices for temperature as a function of pressure and enthalpy
969  double ddT_dTdp_h_constrho = 1/A*CPS.d2hdrhodT()-1/(A*A)*dAdT_constrho*CPS.dhdrho_constT(); //[check]
970  double ddrho_dTdp_h_constT = 1/A*CPS.d2hdrho2_constT()-1/(A*A)*dAdrho_constT*CPS.dhdrho_constT(); //[check
971  double ddT_dTdh_p_constrho = -1/(cp*cp)*(CPS.d2hdT2_constrho()-CPS.dhdp_constT()*CPS.d2pdT2_constrho()+CPS.d2hdrhodT()*CPS.drhodT_constp()-CPS.dhdp_constT()*CPS.drhodT_constp()*CPS.d2pdrhodT());
972  double ddrho_dTdh_p_constT = -1/(cp*cp)*(CPS.d2hdrhodT()-CPS.dhdp_constT()*CPS.d2pdrhodT()+CPS.d2hdrho2_constT()*CPS.drhodT_constp()-CPS.dhdp_constT()*CPS.drhodT_constp()*CPS.d2pdrho2_constT());
973 
974  d2Tdh2[i][j] = ddT_dTdh_p_constrho/CPS.dhdT_constp()+ddrho_dTdh_p_constT/CPS.dhdrho_constp();
975  d2Tdhdp[i][j] = ddT_dTdp_h_constrho/CPS.dhdT_constp()+ddrho_dTdp_h_constT/CPS.dhdrho_constp();
976  d2Tdp2[i][j] = ddT_dTdp_h_constrho/CPS.dpdT_consth()+ddrho_dTdp_h_constT/CPS.dpdrho_consth();
977 
979  double ddT_drhodp_h_constrho = -1/A*CPS.d2hdT2_constrho()+1/(A*A)*dAdT_constrho*CPS.dhdT_constrho();
980  double ddrho_drhodp_h_constT = -1/A*CPS.d2hdrhodT()+1/(A*A)*dAdrho_constT*CPS.dhdT_constrho();
981  double ddT_drhodh_p_constrho = 1/A*CPS.d2pdT2_constrho()-1/(A*A)*dAdT_constrho*CPS.dpdT_constrho();
982  double ddrho_drhodh_p_constT = 1/A*CPS.d2pdrhodT()-1/(A*A)*dAdrho_constT*CPS.dpdT_constrho();
983 
984  d2rhodh2[i][j] = ddT_drhodh_p_constrho/CPS.dhdT_constp()+ddrho_drhodh_p_constT/CPS.dhdrho_constp();
985  d2rhodhdp[i][j] = ddT_drhodp_h_constrho/CPS.dhdT_constp()+ddrho_drhodp_h_constT/CPS.dhdrho_constp();
986  d2rhodp2[i][j] = ddT_drhodp_h_constrho/CPS.dpdT_consth()+ddrho_drhodp_h_constT/CPS.dpdrho_consth();
987  }
988  catch(std::exception &)
989  {
990  s[i][j] = _HUGE;
991  dsdh[i][j] = _HUGE;
992  dsdp[i][j] = _HUGE;
993  d2sdh2[i][j] = _HUGE;
994  d2sdhdp[i][j] = _HUGE;
995  d2sdp2[i][j] = _HUGE;
996 
997  this->T[i][j] = _HUGE;
998  dTdh[i][j] = _HUGE;
999  dTdp[i][j] = _HUGE;
1000  d2Tdh2[i][j] = _HUGE;
1001  d2Tdhdp[i][j] = _HUGE;
1002  d2Tdp2[i][j] = _HUGE;
1003 
1004  this->rho[i][j] = _HUGE;
1005  drhodh[i][j] = _HUGE;
1006  drhodp[i][j] = _HUGE;
1007  d2rhodh2[i][j] = _HUGE;
1008  d2rhodhdp[i][j] = _HUGE;
1009  d2rhodp2[i][j] = _HUGE;
1010  }
1011  }
1012  else
1013  {
1014  s[i][j] = _HUGE;
1015  dsdh[i][j] = _HUGE;
1016  dsdp[i][j] = _HUGE;
1017  d2sdh2[i][j] = _HUGE;
1018  d2sdhdp[i][j] = _HUGE;
1019  d2sdp2[i][j] = _HUGE;
1020 
1021  this->T[i][j] = _HUGE;
1022  dTdh[i][j] = _HUGE;
1023  dTdp[i][j] = _HUGE;
1024  d2Tdh2[i][j] = _HUGE;
1025  d2Tdhdp[i][j] = _HUGE;
1026  d2Tdp2[i][j] = _HUGE;
1027 
1028  this->rho[i][j] = _HUGE;
1029  drhodh[i][j] = _HUGE;
1030  drhodp[i][j] = _HUGE;
1031  d2rhodh2[i][j] = _HUGE;
1032  d2rhodhdp[i][j] = _HUGE;
1033  d2rhodp2[i][j] = _HUGE;
1034 
1035 
1036  }
1037  }
1038  }
1039  t2 = clock();
1040  double elap = (double)(t2-t1)/CLOCKS_PER_SEC;
1041  std::cout << elap << " to build single phase table with p,h" << std::endl;
1042 
1043  // Update the boundaries of the points within the single-phase regions
1044  update_saturation_boundary_indices();
1045 
1046  // Update the cell validity for all cells
1047  update_cell_validity();
1048 
1049  return elap;
1050 }
1051 double TTSESinglePhaseTableClass::build_Trho(double Tmin, double Tmax, double rhomin, double rhomax, TTSETwoPhaseTableClass *SatL, TTSETwoPhaseTableClass *SatV)
1052 {
1053  bool SinglePhase = false;
1054 
1055  if (Tmin < 0 && Tmax < 0 && rhomin < 0 && rhomax < 0)
1056  {
1057  rhomin = 9e9;
1058  rhomax = 0;
1059  Tmin = 9e9;
1060  Tmax = 0;
1061  // Use single-phase table to figure out the range for T,rho
1062  for (unsigned int i = 0; i<Nh; i++)
1063  {
1064  for (unsigned int j = 0; j<Np; j++)
1065  {
1066  if (ValidNumber(rho[i][j]) && rho[i][j] > rhomax){
1067  rhomax = rho[i][j];
1068  }
1069  if (ValidNumber(rho[i][j]) && rho[i][j] < rhomin){
1070  rhomin = rho[i][j];
1071  }
1072  if (ValidNumber(T[i][j]) && T[i][j] > Tmax){
1073  Tmax = T[i][j];
1074  }
1075  if (ValidNumber(T[i][j]) && T[i][j] < Tmin){
1076  Tmin = T[i][j];
1077  }
1078  }
1079  }
1080  }
1081  if (Tmin < pFluid->limits.Tmin){ Tmin = pFluid->limits.Tmin; }
1082  this->Tmin = Tmin;
1083  this->Tmax = Tmax;
1084  this->rhomin = rhomin;
1085  this->rhomax = rhomax;
1086  this->logrhomin = log(rhomin);
1087 
1088  rhoratio = pow(rhomax/rhomin,1/((double)Nrho-1));
1089  logrhoratio = log(rhoratio);
1090 
1091  CoolPropStateClassSI CPS(pFluid);
1092 
1093  long iFluid = get_Fluid_index(pFluid->get_name());
1094 
1095  double dT = (Tmax - Tmin)/((double)NT - 1);
1096 
1097  clock_t t1,t2;
1098  t1 = clock();
1099 
1100  for (unsigned int i = 0; i<NT; i++)
1101  {
1102  double Tval = Tmin + i*dT;
1103  T_Trho[i] = Tval;
1104  for (unsigned int j = 0; j<Nrho; j++)
1105  {
1106  double rhoval = rhomin*pow(rhoratio,(int)j);
1107  rho_Trho[j] = rhoval;
1108 
1109  // Check whether the point is single phase
1110  // If pressure between Ttriple point and Tcrit, might be two-phase or single phase, otherwise definitely single phase
1111  if (Tval <= pFluid->crit.T && Tval >= pFluid->params.Ttriple)
1112  {
1113  if (SatL == NULL || SatV == NULL){
1114  // Not using TTSE method, use saturation (slow...)
1115  CPS.update(iT,Tval,iQ,0.5);
1116  SinglePhase = (rhoval < CPS.rhoV() || rhoval > CPS.rhoL());
1117  }
1118  else{
1119  // Using the TTSE method, nice and fast
1120  double psatV = SatV->evaluate_T(Tval);
1121  double psatL = SatL->evaluate_T(Tval);
1122  SinglePhase = (rhoval < SatV->evaluate(iD,psatV) || rhoval > SatL->evaluate(iD,psatL));
1123  }
1124  }
1125  else
1126  {
1127  SinglePhase = true;
1128  }
1129 
1130  // If enthalpy is outside the saturation region, continue and do the calculation as a function of T,rho
1131  if (SinglePhase)
1132  {
1133  CPS.update(iT,Tval,iD,rhoval);
1134 
1135  s_Trho[i][j] = CPS.s();
1136  dsdT_Trho[i][j] = CPS.dsdT_constrho();
1137  dsdrho_Trho[i][j] = CPS.dsdrho_constT();
1138  d2sdT2_Trho[i][j] = CPS.d2sdT2_constrho();
1139  d2sdTdrho_Trho[i][j] = CPS.d2sdrhodT();
1140  d2sdrho2_Trho[i][j] = CPS.d2sdrho2_constT();
1141 
1142  h_Trho[i][j] = CPS.h();
1143  dhdT_Trho[i][j] = CPS.dhdT_constrho();
1144  dhdrho_Trho[i][j] = CPS.dhdrho_constT();
1145  d2hdT2_Trho[i][j] = CPS.d2hdT2_constrho();
1146  d2hdTdrho_Trho[i][j] = CPS.d2hdrhodT();
1147  d2hdrho2_Trho[i][j] = CPS.d2hdrho2_constT();
1148 
1149  p_Trho[i][j] = CPS.p();
1150  dpdT_Trho[i][j] = CPS.dpdT_constrho();
1151  dpdrho_Trho[i][j] = CPS.dpdrho_constT();
1152  d2pdT2_Trho[i][j] = CPS.d2pdT2_constrho();
1153  d2pdTdrho_Trho[i][j] = CPS.d2pdrhodT();
1154  d2pdrho2_Trho[i][j] = CPS.d2pdrho2_constT();
1155 
1156 
1159  // Using second-order centered finite differences to calculate the transport property derivatives
1160  double deltaT = 1e-3, deltarho = 1e-4;
1161 
1162  if (deltarho > rhoval)
1163  deltarho = rhoval/100;
1164 
1165  // The viscosity values
1166  double muval = IPropsSI(iV,iT,Tval, iD,rhoval, iFluid);
1167  double muplusrho = IPropsSI(iV,iT,Tval, iD,rhoval+deltarho,iFluid);
1168  double muminusrho = IPropsSI(iV,iT,Tval, iD,rhoval-deltarho,iFluid);
1169  double muplusT = IPropsSI(iV,iT,Tval+deltaT,iD,rhoval, iFluid);
1170  double muminusT = IPropsSI(iV,iT,Tval-deltaT,iD,rhoval, iFluid);
1171  double muplusT_plusrho = IPropsSI(iV,iT,Tval+deltaT,iD,rhoval+deltarho,iFluid);
1172  double muplusT_minusrho = IPropsSI(iV,iT,Tval+deltaT,iD,rhoval-deltarho,iFluid);
1173  double muminusT_plusrho = IPropsSI(iV,iT,Tval-deltaT,iD,rhoval+deltarho,iFluid);
1174  double muminusT_minusrho = IPropsSI(iV,iT,Tval-deltaT,iD,rhoval-deltarho,iFluid);
1175 
1176  mu_Trho[i][j] = muval;
1177  dmudT_Trho[i][j] = (-muminusT + muplusT)/(2*deltaT);
1178  dmudrho_Trho[i][j] = (-muminusrho + muplusrho)/(2*deltarho);
1179  d2mudT2_Trho[i][j] = (muminusT - 2*muval + muplusT)/(deltaT*deltaT);
1180  d2mudrho2_Trho[i][j] = (muminusrho - 2*muval + muplusrho)/(deltarho*deltarho);
1181  d2mudTdrho_Trho[i][j] = (muplusT_plusrho - muplusT_minusrho - muminusT_plusrho + muminusT_minusrho)/(2*deltaT*deltarho);
1182 
1183  // The thermal conductivity values
1184  double kval = IPropsSI(iL,iT,Tval, iD,rhoval, iFluid);
1185  double kplusrho = IPropsSI(iL,iT,Tval, iD,rhoval+deltarho,iFluid);
1186  double kminusrho = IPropsSI(iL,iT,Tval, iD,rhoval-deltarho,iFluid);
1187  double kplusT = IPropsSI(iL,iT,Tval+deltaT,iD,rhoval, iFluid);
1188  double kminusT = IPropsSI(iL,iT,Tval-deltaT,iD,rhoval, iFluid);
1189  double kplusT_plusrho = IPropsSI(iL,iT,Tval+deltaT,iD,rhoval+deltarho,iFluid);
1190  double kplusT_minusrho = IPropsSI(iL,iT,Tval+deltaT,iD,rhoval-deltarho,iFluid);
1191  double kminusT_plusrho = IPropsSI(iL,iT,Tval-deltaT,iD,rhoval+deltarho,iFluid);
1192  double kminusT_minusrho = IPropsSI(iL,iT,Tval-deltaT,iD,rhoval-deltarho,iFluid);
1193 
1194  k_Trho[i][j] = kval;
1195  dkdT_Trho[i][j] = (-kminusT + kplusT)/(2*deltaT);
1196  dkdrho_Trho[i][j] = (-kminusrho + kplusrho)/(2*deltarho);
1197  d2kdT2_Trho[i][j] = (kminusT - 2*kval + kplusT)/(deltaT*deltaT);
1198  d2kdrho2_Trho[i][j] = (kminusrho - 2*kval + kplusrho)/(deltarho*deltarho);
1199  d2kdTdrho_Trho[i][j] = (kplusT_plusrho - kplusT_minusrho - kminusT_plusrho + kminusT_minusrho)/(2*deltaT*deltarho);
1200  }
1201  else
1202  {
1203  s_Trho[i][j] = _HUGE;
1204  dsdT_Trho[i][j] = _HUGE;
1205  dsdrho_Trho[i][j] = _HUGE;
1206  d2sdT2_Trho[i][j] = _HUGE;
1207  d2sdTdrho_Trho[i][j] = _HUGE;
1208  d2sdrho2_Trho[i][j] = _HUGE;
1209 
1210  h_Trho[i][j] = _HUGE;
1211  dhdT_Trho[i][j] = _HUGE;
1212  dhdrho_Trho[i][j] = _HUGE;
1213  d2hdT2_Trho[i][j] = _HUGE;
1214  d2hdTdrho_Trho[i][j] = _HUGE;
1215  d2hdrho2_Trho[i][j] = _HUGE;
1216 
1217  p_Trho[i][j] = _HUGE;
1218  dpdT_Trho[i][j] = _HUGE;
1219  dpdrho_Trho[i][j] = _HUGE;
1220  d2pdT2_Trho[i][j] = _HUGE;
1221  d2pdTdrho_Trho[i][j] = _HUGE;
1222  d2pdrho2_Trho[i][j] = _HUGE;
1223 
1224  mu_Trho[i][j] = _HUGE;
1225  dmudT_Trho[i][j] = _HUGE;
1226  dmudrho_Trho[i][j] = _HUGE;
1227  d2mudT2_Trho[i][j] = _HUGE;
1228  d2mudTdrho_Trho[i][j] = _HUGE;
1229  d2mudrho2_Trho[i][j] = _HUGE;
1230 
1231  k_Trho[i][j] = _HUGE;
1232  dkdT_Trho[i][j] = _HUGE;
1233  dkdrho_Trho[i][j] = _HUGE;
1234  d2kdT2_Trho[i][j] = _HUGE;
1235  d2kdTdrho_Trho[i][j] = _HUGE;
1236  d2kdrho2_Trho[i][j] = _HUGE;
1237  }
1238  //std::cout << format("%d %d\n",i,j);
1239  }
1240  }
1241  t2 = clock();
1242  double elap = (double)(t2-t1)/CLOCKS_PER_SEC;
1243  std::cout << elap << " to build single phase table for T,rho" << std::endl;
1244 
1245  // Update the boundaries of the points within the single-phase regions
1246  update_saturation_boundary_indices();
1247 
1248  update_cell_validity();
1249 
1250  return elap;
1251 }
1252 
1253 void TTSESinglePhaseTableClass::update_cell_validity()
1254 {
1255  for (unsigned int i = 0; i < NT; i++)
1256  {
1257  for (unsigned int j = 0; j < Nrho; j++)
1258  {
1259  if (i < NT-1 && j < Nrho-1)
1260  {
1261  if(ValidNumber(s_Trho[i][j]) && ValidNumber(s_Trho[i][j+1]) && ValidNumber(s_Trho[i+1][j]) && ValidNumber(s_Trho[i+1][j+1]))
1262  {
1263  bicubic_cells.cells[i][j].valid_Trho = true;
1264  }
1265  else
1266  {
1267  bicubic_cells.cells[i][j].valid_Trho = false;
1268  }
1269  }
1270  }
1271  }
1272 
1273  for (unsigned int i = 0; i < Nh; i++)
1274  {
1275  for (unsigned int j = 0; j < Np; j++)
1276  {
1277  if (i < Nh-1 && j < Np-1)
1278  {
1279  if(ValidNumber(s[i][j]) && ValidNumber(s[i][j+1]) && ValidNumber(s[i+1][j]) && ValidNumber(s[i+1][j+1]))
1280  {
1281  bicubic_cells.cells[i][j].valid_hp = true;
1282  }
1283  else
1284  {
1285  bicubic_cells.cells[i][j].valid_hp = false;
1286  }
1287  }
1288  }
1289  }
1290 }
1291 //void TTSESinglePhaseTableClass::update_Trho_map()
1292 //{
1293 // int ii,jj;
1294 // double Tmin,Tmax,rhomin,rhomax,rhoL,rhoV,TsatL,TsatV,dummy;
1295 // // Get the bounding values
1296 // Tmin = T[0][0];
1297 // rhomax = rho[0][0];
1298 // Tmax = T[Nh-1][Np-1];
1299 // rhomin = rho[Nh-1][0];
1300 // // Resize the arrays, using the same sizes as the base matrices
1301 // T_Trho.resize(Nh);
1302 // rho_Trho.resize(Np);
1303 // i_Trho.resize(Nh, std::vector<int>(Np, -1));
1304 // j_Trho.resize(Nh, std::vector<int>(Np, -1));
1305 //
1306 // for (unsigned int i = 0; i < Nh; i++)
1307 // {
1308 // double T = (Tmax-Tmin)/(Nh-1)*i+Tmin;
1309 // T_Trho[i] = T;
1310 //
1311 // for (unsigned int j = 0; j < Np; j++)
1312 // {
1313 // double rho = (rhomax-rhomin)/(Np-1)*j+rhomin;
1314 // rho_Trho[j] = rho;
1315 //
1316 // // T,rho --> p,h
1317 // double p = pFluid->pressure_Trho(T,rho);
1318 // double h = pFluid->enthalpy_Trho(T,rho);
1319 // double rhooV = pFluid->rhosatV(T);
1320 // double rhooL = pFluid->rhosatL(T);
1321 // double pV = pFluid->psatV_anc(T);
1322 // double pp = pFluid->pressure_Trho(T,rhooV);
1323 //
1324 // // Find i,j from p,h
1325 // ii = (int)round(((h-hmin)/(hmax-hmin)*(Nh-1)));
1326 // jj = (int)round((log(p)-logpmin)/logpratio);
1327 //
1328 // // Only keep values that are within the range for the table
1329 // if ( ii>=0 && ii < (int)Nh && jj>=0 && jj< (int)Np)
1330 // {
1331 // i_Trho[i][j] = ii;
1332 // j_Trho[i][j] = jj;
1333 // }
1334 // else
1335 // {
1336 // i_Trho[i][j] = -1;
1337 // j_Trho[i][j] = -1;
1338 // }
1339 // }
1340 // }
1341 //}
1342 void TTSESinglePhaseTableClass::update_saturation_boundary_indices()
1343 {
1344  // Store some information about the phase boundaries so that we
1345  // can use other inputs than p,h more easily
1346 
1347  for (unsigned int j = 0; j < Np; j++)
1348  {
1349  if (p[j] < pFluid->reduce.p.Pa)
1350  {
1351  IL[j] = -1;
1352  // Sweep left to right to find a phase boundary, use the first one that fails in the saturation region
1353  for (unsigned int i = 0; i < Nh; i++)
1354  {
1355  if (!ValidNumber(T[i][j]))
1356  {
1357  IL[j] = i;
1358  break;
1359  }
1360  }
1361  IV[j] = -1;
1362  // Sweep right to left to find a phase boundary, use the first one that fails in the saturation region
1363  for (int i = Nh-1; i > 0; i--)
1364  {
1365  if (!ValidNumber(T[i][j]))
1366  {
1367  IV[j] = i;
1368  break;
1369  }
1370  }
1371  if (SatL != NULL && SatV != NULL)
1372  {
1373  TL[j] = SatL->evaluate(iT,p[j]);
1374  SL[j] = SatL->evaluate(iS,p[j]);
1375  DL[j] = SatL->evaluate(iD,p[j]);
1376  TV[j] = SatV->evaluate(iT,p[j]);
1377  SV[j] = SatV->evaluate(iS,p[j]);
1378  DV[j] = SatV->evaluate(iD,p[j]);
1379  }
1380  else
1381  {
1382  throw ValueError("SatL and SatV must be provided");
1383  }
1384  }
1385  else
1386  {
1387  IL[j] = -1;
1388  IV[j] = -1;
1389  TL[j] = _HUGE;
1390  SL[j] = _HUGE;
1391  DL[j] = _HUGE;
1392  TV[j] = _HUGE;
1393  SV[j] = _HUGE;
1394  DV[j] = _HUGE;
1395  }
1396  }
1397 }
1398 void TTSESinglePhaseTableClass::write_dotdrawing_tofile(char fName[])
1399 {
1400  FILE *fp;
1401  fp = fopen(fName,"w");
1402  for (int j = Np-1; j>=0; j--)
1403  {
1404  for (unsigned int i = 0; i<Nh; i++)
1405  {
1406  if (ValidNumber(rho[i][j]))
1407  {
1408  fprintf(fp,".");
1409  }
1410  else
1411  {
1412  fprintf(fp,"X");
1413  }
1414  }
1415  fprintf(fp,"\n");
1416  }
1417  fclose(fp);
1418 }
1419 
1420 double TTSESinglePhaseTableClass::check_randomly(long iParam, unsigned int N, std::vector<double> *h, std::vector<double> *p, std::vector<double> *EOSv, std::vector<double> *TTSE)
1421 {
1422  double val=0;
1423  h->resize(N);
1424  p->resize(N);
1425  EOSv->resize(N);
1426  TTSE->resize(N);
1427 
1428  CoolPropStateClassSI CPS(pFluid);
1429 
1430  for (unsigned int i = 0; i < N; i++)
1431  {
1432  double p1 = ((double)rand()/(double)RAND_MAX)*(pmax-pmin)+pmin;
1433  double h1 = ((double)rand()/(double)RAND_MAX)*(hmax-hmin)+hmin;
1434 
1435  CPS.update(iH,h1,iP,p1);
1436  double sEOS = CPS.s();
1437  double cpEOS = CPS.cp();
1438  double TEOS = CPS.T();
1439  double rhoEOS = CPS.rho();
1440 
1441  // Store the inputs
1442  (*h)[i] = h1;
1443  (*p)[i] = p1;
1444 
1445  // Get the value from TTSE
1446  (*TTSE)[i] = evaluate(iParam,p1,CPS._logp,h1);
1447 
1448  // Get the value from EOS
1449  switch (iParam)
1450  {
1451  case iS:
1452  (*EOSv)[i] = sEOS; break;
1453  case iT:
1454  (*EOSv)[i] = TEOS; break;
1455  case iC:
1456  (*EOSv)[i] = cpEOS; break;
1457  case iD:
1458  (*EOSv)[i] = rhoEOS; break;
1459  default:
1460  throw ValueError();
1461  }
1462 
1463  std::cout << format("%g %g %g %g %g (h,p,EOS,TTSE,diff)\n",h1,p1,(*EOSv)[i],(*TTSE)[i],(*EOSv)[i]-(*TTSE)[i]).c_str();
1464  }
1465  return val;
1466 }
1467 
1468 double TTSESinglePhaseTableClass::evaluate_randomly(long iParam, unsigned int N)
1469 {
1470  clock_t t1,t2;
1471  t1 = clock();
1472  for (unsigned int i = 0; i < N; i++)
1473  {
1474  double p1 = ((double)rand()/(double)RAND_MAX)*(pmax-pmin)+pmin;
1475  double h1 = ((double)rand()/(double)RAND_MAX)*(hmax-hmin)+hmin;
1476 
1477  if (p1 > pFluid->TTSESatL.pmax || h1 > pFluid->TTSESatV.evaluate(iH,p1) || h1 < pFluid->TTSESatL.evaluate(iH,p1))
1478  {
1479  // Get the value from TTSE
1480  evaluate(iParam,p1,log(p1),h1);
1481  }
1482  }
1483  t2 = clock();
1484  return (double)(t2-t1)/CLOCKS_PER_SEC/(double)N*1e6;
1485 }
1486 
1487 double TTSESinglePhaseTableClass::evaluate(long iParam, double p, double logp, double h)
1488 {
1489  if (!this->within_range(iP,p,iH,h))
1490  {
1491  throw ValueError(format("Input to TTSE [p = %0.16g, h = %0.16g] is out of range",p,h));
1492  }
1493 
1494  // Use Bicubic interpolation if requested
1495  if (mode == TTSE_MODE_BICUBIC){
1496  return interpolate_bicubic_ph(iParam,p,logp,h);
1497  }
1498 
1499  int i = (int)round(((h-hmin)/(hmax-hmin)*(Nh-1)));
1500  int j = (int)round((logp-logpmin)/logpratio);
1501 
1502  // If the value at i,j is too close to the saturation boundary, the nearest point i,j
1503  // might be in the two-phase region which is not defined for single-phase table.
1504  // Therefore, search around its neighbors for a better choice
1505  if (!ValidNumber(T[i][j])){
1506  nearest_good_neighbor(&i,&j);
1507  }
1508 
1509  // Distances from the node
1510  double deltap = p-this->p[j];
1511  double deltah = h-this->h[i];
1512 
1513  switch (iParam)
1514  {
1515  case iS:
1516  return s[i][j]+deltah*dsdh[i][j]+deltap*dsdp[i][j]+0.5*deltah*deltah*d2sdh2[i][j]+0.5*deltap*deltap*d2sdp2[i][j]+deltap*deltah*d2sdhdp[i][j]; break;
1517  case iT:
1518  return T[i][j]+deltah*dTdh[i][j]+deltap*dTdp[i][j]+0.5*deltah*deltah*d2Tdh2[i][j]+0.5*deltap*deltap*d2Tdp2[i][j]+deltap*deltah*d2Tdhdp[i][j]; break;
1519  case iD:
1520  return rho[i][j]+deltah*drhodh[i][j]+deltap*drhodp[i][j]+0.5*deltah*deltah*d2rhodh2[i][j]+0.5*deltap*deltap*d2rhodp2[i][j]+deltap*deltah*d2rhodhdp[i][j]; break;
1521  default:
1522  throw ValueError(format("Output key value [%d] to evaluate is invalid",iParam));
1523  }
1524 }
1525 
1526 void TTSESinglePhaseTableClass::bicubic_cell_coordinates_Trho(double Tval, double rhoval, double logrhoval, int *i, int *j)
1527 {
1528  if (!this->within_range(iT,Tval,iD,rhoval))
1529  {
1530  throw ValueError(format("Input to TTSE [T = %0.16g, logrho = %0.16g] is out of range",Tval,logrhoval));
1531  }
1532 
1533  *i = (int)round((Tval-Tmin)/(Tmax-Tmin)*(NT-1));
1534  *j = (int)round((logrhoval-logrhomin)/logrhoratio);
1535 
1536  if (Tval < this->T_Trho[*i])
1537  {
1538  *i -= 1;
1539  }
1540  if (rhoval < this->rho_Trho[*j])
1541  {
1542  *j -= 1;
1543  }
1544 
1545  nearest_good_neighbor_Trho_interpolate(i, j);
1546 }
1547 
1548 void TTSESinglePhaseTableClass::bicubic_cell_coordinates_ph(double hval, double pval, double logpval, int *i, int *j)
1549 {
1550  if (!this->within_range(iP,pval,iH,hval))
1551  {
1552  throw ValueError(format("Input to TTSE [p = %0.16g, h = %0.16g] is out of range",pval,hval));
1553  }
1554 
1555  *i = (int)round(((hval-hmin)/(hmax-hmin)*(Nh-1)));
1556  *j = (int)round((logpval-logpmin)/logpratio);
1557 
1558  if (hval < this->h[*i])
1559  {
1560  *i -= 1;
1561  }
1562  if (pval < p[*j])
1563  {
1564  *j -= 1;
1565  }
1566 
1567  nearest_good_neighbor_ph_interpolate(i, j);
1568 }
1569 
1570 std::vector<double> * TTSESinglePhaseTableClass::bicubic_cell_coeffs_ph(long iParam, int i, int j)
1571 {
1572  std::vector<double> *alpha = NULL;
1573  std::vector< std::vector<double> > *f = NULL, *dfdp = NULL, *dfdh = NULL, *d2fdhdp = NULL;
1574 
1575  switch(iParam)
1576  {
1577  case iS:
1578  if (!bicubic_cells.cells[i][j].alpha_s_hp.empty())
1579  {
1580  // Make alpha point to the pre-calculated values
1581  alpha = &bicubic_cells.cells[i][j].alpha_s_hp;
1582  break;
1583  }
1584  else
1585  {
1586  f = &s; dfdh = &dsdh; dfdp = &dsdp; d2fdhdp = &d2sdhdp; break;
1587  }
1588  case iT:
1589  if (!bicubic_cells.cells[i][j].alpha_T_hp.empty())
1590  {
1591  // Make alpha point to the pre-calculated values
1592  alpha = &bicubic_cells.cells[i][j].alpha_T_hp;
1593  break;
1594  }
1595  else
1596  {
1597  f = &T; dfdh = &dTdh; dfdp = &dTdp; d2fdhdp = &d2Tdhdp; break;
1598  }
1599  case iD:
1600  if (!bicubic_cells.cells[i][j].alpha_rho_hp.empty())
1601  {
1602  // Make alpha point to the pre-calculated values
1603  alpha = &bicubic_cells.cells[i][j].alpha_rho_hp;
1604  break;
1605  }
1606  else
1607  {
1608  f = &rho; dfdh = &drhodh; dfdp = &drhodp; d2fdhdp = &d2rhodhdp; break;
1609  }
1610  default:
1611  throw ValueError("iParam is invalid");
1612  }
1613 
1614  if (alpha == NULL)
1615  {
1616  double dhdx = (h[i+1]-h[i]);
1617  double dpdy = (p[j+1]-p[j]);
1618 
1619  z_bicubic[0] = (*f)[i][j];
1620  z_bicubic[1] = (*f)[i+1][j];
1621  z_bicubic[2] = (*f)[i][j+1];
1622  z_bicubic[3] = (*f)[i+1][j+1];
1623 
1624  z_bicubic[4] = (*dfdh)[i][j]*dhdx;
1625  z_bicubic[5] = (*dfdh)[i+1][j]*dhdx;
1626  z_bicubic[6] = (*dfdh)[i][j+1]*dhdx;
1627  z_bicubic[7] = (*dfdh)[i+1][j+1]*dhdx;
1628 
1629  z_bicubic[8] = (*dfdp)[i][j]*dpdy;
1630  z_bicubic[9] = (*dfdp)[i+1][j]*dpdy;
1631  z_bicubic[10] = (*dfdp)[i][j+1]*dpdy;
1632  z_bicubic[11] = (*dfdp)[i+1][j+1]*dpdy;
1633 
1634  z_bicubic[12] = (*d2fdhdp)[i][j]*dpdy*dhdx;
1635  z_bicubic[13] = (*d2fdhdp)[i+1][j]*dpdy*dhdx;
1636  z_bicubic[14] = (*d2fdhdp)[i][j+1]*dpdy*dhdx;
1637  z_bicubic[15] = (*d2fdhdp)[i+1][j+1]*dpdy*dhdx;
1638 
1639  // Find the alpha values by doing the matrix operation alpha = Ainv*z
1640  matrix_vector_product(&z_bicubic, &alpha_bicubic);
1641  // Set the pointer in this function
1642  alpha = &alpha_bicubic;
1643 
1644  // Store this array of bicubic coefficients
1645  switch (iParam)
1646  {
1647  case iS:
1648  bicubic_cells.cells[i][j].alpha_s_hp = alpha_bicubic; break;
1649  case iT:
1650  bicubic_cells.cells[i][j].alpha_T_hp = alpha_bicubic; break;
1651  case iD:
1652  bicubic_cells.cells[i][j].alpha_rho_hp = alpha_bicubic; break;
1653  default:
1654  throw ValueError("iParam is invalid");
1655  }
1656  }
1657  return alpha;
1658 }
1659 
1660 std::vector<double> * TTSESinglePhaseTableClass::bicubic_cell_coeffs_Trho(long iParam, int i, int j)
1661 {
1662 
1663  std::vector<double> *alpha = NULL;
1664  std::vector< std::vector<double> > *f = NULL, *dfdT = NULL, *dfdrho = NULL, *d2fdTdrho = NULL;
1665 
1666  // Find the coefficients for the bicubic function
1667  switch(iParam)
1668  {
1669  case iS:
1670  if (!bicubic_cells.cells[i][j].alpha_s_Trho.empty())
1671  {
1672  // Make alpha point to the pre-calculated values
1673  alpha = &bicubic_cells.cells[i][j].alpha_s_Trho;
1674  break;
1675  }
1676  else
1677  {
1678  f = &s_Trho; dfdT = &dsdT_Trho; dfdrho = &dsdrho_Trho; d2fdTdrho = &d2sdTdrho_Trho; break;
1679  }
1680  case iH:
1681  if (!bicubic_cells.cells[i][j].alpha_h_Trho.empty())
1682  {
1683  // Make alpha point to the pre-calculated values
1684  alpha = &bicubic_cells.cells[i][j].alpha_h_Trho;
1685  break;
1686  }
1687  else
1688  {
1689  f = &h_Trho; dfdT = &dhdT_Trho; dfdrho = &dhdrho_Trho; d2fdTdrho = &d2hdTdrho_Trho; break;
1690  }
1691  case iP:
1692  if (!bicubic_cells.cells[i][j].alpha_p_Trho.empty())
1693  {
1694  // Make alpha point to the pre-calculated values
1695  alpha = &bicubic_cells.cells[i][j].alpha_p_Trho;
1696  break;
1697  }
1698  else
1699  {
1700  f = &p_Trho; dfdT = &dpdT_Trho; dfdrho = &dpdrho_Trho; d2fdTdrho = &d2pdTdrho_Trho; break;
1701  }
1702  case iV:
1703  if (!bicubic_cells.cells[i][j].alpha_mu_Trho.empty())
1704  {
1705  // Make alpha point to the pre-calculated values
1706  alpha = &bicubic_cells.cells[i][j].alpha_mu_Trho;
1707  break;
1708  }
1709  else
1710  {
1711  f = &mu_Trho; dfdT = &dmudT_Trho; dfdrho = &dmudrho_Trho; d2fdTdrho = &d2mudTdrho_Trho; break;
1712  }
1713  case iL:
1714  if (!bicubic_cells.cells[i][j].alpha_k_Trho.empty())
1715  {
1716  // Make alpha point to the pre-calculated values
1717  alpha = &bicubic_cells.cells[i][j].alpha_k_Trho;
1718  break;
1719  }
1720  else
1721  {
1722  f = &k_Trho; dfdT = &dkdT_Trho; dfdrho = &dkdrho_Trho; d2fdTdrho = &d2kdTdrho_Trho; break;
1723  }
1724  default:
1725  throw ValueError(format("iParam [%d] is invalid",iParam).c_str());
1726  }
1727 
1728  if (alpha == NULL)
1729  {
1730  double dTdx = (T_Trho[i+1]-T_Trho[i]);
1731  double drhody = (rho_Trho[j+1]-rho_Trho[j]);
1732 
1733  z_bicubic[0] = (*f)[i][j];
1734  z_bicubic[1] = (*f)[i+1][j];
1735  z_bicubic[2] = (*f)[i][j+1];
1736  z_bicubic[3] = (*f)[i+1][j+1];
1737 
1738  z_bicubic[4] = (*dfdT)[i][j]*dTdx;
1739  z_bicubic[5] = (*dfdT)[i+1][j]*dTdx;
1740  z_bicubic[6] = (*dfdT)[i][j+1]*dTdx;
1741  z_bicubic[7] = (*dfdT)[i+1][j+1]*dTdx;
1742 
1743  z_bicubic[8] = (*dfdrho)[i][j]*drhody;
1744  z_bicubic[9] = (*dfdrho)[i+1][j]*drhody;
1745  z_bicubic[10] = (*dfdrho)[i][j+1]*drhody;
1746  z_bicubic[11] = (*dfdrho)[i+1][j+1]*drhody;
1747 
1748  z_bicubic[12] = (*d2fdTdrho)[i][j]*drhody*dTdx;
1749  z_bicubic[13] = (*d2fdTdrho)[i+1][j]*drhody*dTdx;
1750  z_bicubic[14] = (*d2fdTdrho)[i][j+1]*drhody*dTdx;
1751  z_bicubic[15] = (*d2fdTdrho)[i+1][j+1]*drhody*dTdx;
1752 
1753  // Find the alpha values by doing the matrix operation alpha = Ainv*z
1754  matrix_vector_product(&z_bicubic, &alpha_bicubic);
1755  // Set the pointer in this function
1756  alpha = &alpha_bicubic;
1757 
1758  // Store this array of bicubic coefficients
1759  switch (iParam)
1760  {
1761  case iS:
1762  bicubic_cells.cells[i][j].alpha_s_Trho = alpha_bicubic; break;
1763  case iH:
1764  bicubic_cells.cells[i][j].alpha_h_Trho = alpha_bicubic; break;
1765  case iP:
1766  bicubic_cells.cells[i][j].alpha_p_Trho = alpha_bicubic; break;
1767  case iV:
1768  bicubic_cells.cells[i][j].alpha_mu_Trho = alpha_bicubic; break;
1769  case iL:
1770  bicubic_cells.cells[i][j].alpha_k_Trho = alpha_bicubic; break;
1771  default:
1772  throw ValueError(format("iParam [%d] is invalid",iParam).c_str());
1773  }
1774  }
1775 
1776  return alpha;
1777 }
1778 
1825 double TTSESinglePhaseTableClass::interpolate_bicubic_ph(long iParam, double pval, double logp, double hval)
1826 {
1827  std::vector<double> *alpha = NULL;
1828  int i, j;
1829 
1830  // Get the i,j coordinates of the cell
1831  this->bicubic_cell_coordinates_ph(hval, pval,logp, &i, &j);
1832 
1833  double dhdx = (h[i+1]-h[i]);
1834  double dpdy = (p[j+1]-p[j]);
1835  double x = (hval-h[i])/dhdx;
1836  double y = (pval-p[j])/dpdy;
1837 
1838  // Find the coefficients for the cubic
1839  alpha = this->bicubic_cell_coeffs_ph(iParam, i, j);
1840 
1841  double x_0 = 1, x_1 = x, x_2 = x*x, x_3 = x*x*x;
1842  double y_0 = 1, y_1 = y, y_2 = y*y, y_3 = y*y*y;
1843 
1844  double summer;
1845  // Old version was "summer += (*alpha)[ii+jj*4]*x^i*y^j" for i in [0,3] and j in [0,3]
1846  summer = (*alpha)[0+0*4]*x_0*y_0+(*alpha)[0+1*4]*x_0*y_1+(*alpha)[0+2*4]*x_0*y_2+(*alpha)[0+3*4]*x_0*y_3
1847  +(*alpha)[1+0*4]*x_1*y_0+(*alpha)[1+1*4]*x_1*y_1+(*alpha)[1+2*4]*x_1*y_2+(*alpha)[1+3*4]*x_1*y_3
1848  +(*alpha)[2+0*4]*x_2*y_0+(*alpha)[2+1*4]*x_2*y_1+(*alpha)[2+2*4]*x_2*y_2+(*alpha)[2+3*4]*x_2*y_3
1849  +(*alpha)[3+0*4]*x_3*y_0+(*alpha)[3+1*4]*x_3*y_1+(*alpha)[3+2*4]*x_3*y_2+(*alpha)[3+3*4]*x_3*y_3;
1850 
1851  return summer;
1852 }
1853 double TTSESinglePhaseTableClass::bicubic_evaluate_first_derivative_ph(long iOF, long iWRT, long iCONSTANT, double p, double logp, double h)
1854 {
1855  std::vector<double> *alpha = NULL;
1856  int i, j;
1857 
1858  // Get the i,j coordinates of the cell
1859  this->bicubic_cell_coordinates_ph(h, p, logp, &i, &j);
1860 
1861  double dhdx = (this->h[i+1]-this->h[i]);
1862  double dpdy = (this->p[j+1]-this->p[j]);
1863  double x = (h-this->h[i])/dhdx;
1864  double y = (p-this->p[j])/dpdy;
1865 
1866  // Find the coefficients for the cubic
1867  alpha = this->bicubic_cell_coeffs_ph(iOF, i, j);
1868 
1869  double x_0 = 1, x_1 = x, x_2 = x*x, x_3 = x*x*x;
1870  double y_0 = 1, y_1 = y, y_2 = y*y, y_3 = y*y*y;
1871 
1872  double summer;
1873  if (iWRT == iH)
1874  {
1875  // Old version was "summer += (*alpha)[ii+jj*4]*i*x^(i-1)*y^j" for i in [0,3] and j in [0,3]
1876  double dsummer_dx = +(*alpha)[1+0*4]*1*x_0*y_0+(*alpha)[1+1*4]*1*x_0*y_1+(*alpha)[1+2*4]*1*x_0*y_2+(*alpha)[1+3*4]*1*x_0*y_3
1877  +(*alpha)[2+0*4]*2*x_1*y_0+(*alpha)[2+1*4]*2*x_1*y_1+(*alpha)[2+2*4]*2*x_1*y_2+(*alpha)[2+3*4]*2*x_1*y_3
1878  +(*alpha)[3+0*4]*3*x_2*y_0+(*alpha)[3+1*4]*3*x_2*y_1+(*alpha)[3+2*4]*3*x_2*y_2+(*alpha)[3+3*4]*3*x_2*y_3;
1879  return dsummer_dx/dhdx;
1880  }
1881  else if (iWRT == iP)
1882  {
1883  // Old version was "summer += (*alpha)[ii+jj*4]*j*x^i*y^(j-1)" for i in [0,3] and j in [0,3]
1884  double dsummer_dy = (*alpha)[0+1*4]*1*x_0*y_0+(*alpha)[0+2*4]*2*x_0*y_1+(*alpha)[0+3*4]*3*x_0*y_2
1885  +(*alpha)[1+1*4]*1*x_1*y_0+(*alpha)[1+2*4]*2*x_1*y_1+(*alpha)[1+3*4]*3*x_1*y_2
1886  +(*alpha)[2+1*4]*1*x_2*y_0+(*alpha)[2+2*4]*2*x_2*y_1+(*alpha)[2+3*4]*3*x_2*y_2
1887  +(*alpha)[3+1*4]*1*x_3*y_0+(*alpha)[3+2*4]*2*x_3*y_1+(*alpha)[3+3*4]*3*x_3*y_2;
1888  return dsummer_dy/dpdy;
1889  }
1890  else
1891  {
1892  throw ValueError(format("iWRT [%d] is invalid", iWRT).c_str());
1893  }
1894 
1895  return summer;
1896 }
1897 
1898 double TTSESinglePhaseTableClass::interpolate_bicubic_Trho(long iParam, double Tval, double rhoval, double logrho)
1899 {
1900  // A pointer to the values of the coefficients for the bicubic interpolation
1901  std::vector<double> *alpha = NULL;
1902  int i,j;
1903  //std::vector< std::vector<double> > *f = NULL, *dfdT = NULL, *dfdrho = NULL, *d2fdTdrho = NULL;
1904 
1905  bicubic_cell_coordinates_Trho(Tval, rhoval, logrho, &i, &j);
1906 
1907  double dTdx = (T_Trho[i+1]-T_Trho[i]);
1908  double drhody = (rho_Trho[j+1]-rho_Trho[j]);
1909  double x = (Tval-T_Trho[i])/dTdx;
1910  double y = (rhoval-rho_Trho[j])/drhody;
1911 
1912  // Find the coefficients for the cubic
1913  alpha = this->bicubic_cell_coeffs_Trho(iParam, i, j);
1914 
1915  double x_0 = 1, x_1 = x, x_2 = x*x, x_3 = x*x*x;
1916  double y_0 = 1, y_1 = y, y_2 = y*y, y_3 = y*y*y;
1917 
1918  double summer;
1919  // Old version was "summer += (*alpha)[ii+jj*4]*x^i*y^j" for i in [0,3] and j in [0,3]
1920  summer = (*alpha)[0+0*4]*x_0*y_0+(*alpha)[0+1*4]*x_0*y_1+(*alpha)[0+2*4]*x_0*y_2+(*alpha)[0+3*4]*x_0*y_3
1921  +(*alpha)[1+0*4]*x_1*y_0+(*alpha)[1+1*4]*x_1*y_1+(*alpha)[1+2*4]*x_1*y_2+(*alpha)[1+3*4]*x_1*y_3
1922  +(*alpha)[2+0*4]*x_2*y_0+(*alpha)[2+1*4]*x_2*y_1+(*alpha)[2+2*4]*x_2*y_2+(*alpha)[2+3*4]*x_2*y_3
1923  +(*alpha)[3+0*4]*x_3*y_0+(*alpha)[3+1*4]*x_3*y_1+(*alpha)[3+2*4]*x_3*y_2+(*alpha)[3+3*4]*x_3*y_3;
1924 
1925  return summer;
1926 }
1927 
1928 bool isbetween(double x1, double x2, double x)
1929 {
1930  return ((x>=x1 && x <= x2) || (x>=x2 && x <= x1));
1931 }
1932 
1933 double TTSESinglePhaseTableClass::bicubic_evaluate_one_other_input(long iInput1, double Input1, long iOther, double Other)
1934 {
1935  // My goal here is to find deltah
1936  int L,R,M,i,j,phase;
1937  double p,TL, TV, DL, DV, SL, SV, HL, HV, Q, ValV, ValL;
1938  std::vector<std::vector<double> > *mat;
1939 
1940  // Connect a pointer to the array of interest
1941  switch (iOther)
1942  {
1943  case iT:
1944  mat = &T; break;
1945  case iS:
1946  mat = &s; break;
1947  case iD:
1948  mat = &rho; break;
1949  }
1950 
1951  // One is pressure, we are getting enthalpy
1952  if (iInput1 == iP)
1953  {
1954  p = Input1;
1955 
1956  j = (int)round((log(p)-logpmin)/logpratio);
1957 
1958  if (p > pFluid->reduce.p.Pa)
1959  {
1960  // It's supercritical, we just need to iterate on the other parameter
1961 
1962  // Do interval halving over the whole range since
1963  // there can't be any saturation curve
1964  L = 0; R = Nh-1; M = (L+R)/2;
1965  while (R-L>1){
1966  if (isbetween((*mat)[M][j],(*mat)[R][j],Other)){
1967  L=M; M=(L+R)/2; continue;
1968  }
1969  else{
1970  R=M; M=(L+R)/2; continue;
1971  }
1972  }
1973  i = L;
1974  }
1975  else
1976  {
1977  // Get the saturation states
1978  switch (iOther)
1979  {
1980  case iT:
1981  TL = SatL->evaluate(iT,p); TV = SatV->evaluate(iT,p); ValL = TL; ValV = TV;
1982  if (Other > TV){ phase = iPHASE_GAS;} else if (Other < TL){ phase = iPHASE_LIQUID;} else {phase = iPHASE_TWOPHASE; throw ValueError("Two-phase T,P inputs invalid for TTSE");}
1983  break;
1984  case iD:
1985  DL = SatL->evaluate(iD,p); DV = SatV->evaluate(iD,p); ValL = DL; ValV = DV;
1986  if (Other < DV){ phase = iPHASE_GAS;} else if (Other > DL){ phase = iPHASE_LIQUID;} else {phase = iPHASE_TWOPHASE; Q = (1/Other-1/DL)/(1/DV-1/DL);}
1987  break;
1988  case iS:
1989  SL = SatL->evaluate(iS,p); SV = SatV->evaluate(iS,p); ValL = SL; ValV = SV;
1990  if (Other > SV){ phase = iPHASE_GAS;} else if (Other < SL){ phase = iPHASE_LIQUID;} else {phase = iPHASE_TWOPHASE; Q = (Other-SL)/(SV-SL);}
1991  break;
1992  }
1993 
1994  if (phase == iPHASE_GAS)
1995  {
1996  L = IV[j]+1; R = Nh - 1; M = (L+R)/2;
1997  // Check if caught between saturation curve and first point
1998  if (isbetween((*mat)[L][j],ValV,Other))
1999  {
2000  // If caught, use gas value
2001  L = R;
2002  }
2003  }
2004  else if (phase == iPHASE_LIQUID)
2005  {
2006  L = 0; R = IL[j]-1; M = (L+R)/2;
2007  // Check if caught between saturation curve and first point
2008  if (isbetween((*mat)[R][j],ValL,Other))
2009  {
2010  // If caught, use liquid value
2011  R = L;
2012  }
2013  }
2014  else if (phase == iPHASE_TWOPHASE)
2015  {
2016  // Two-phase, finished
2017  HL = SatL->evaluate(iH,p); HV = SatV->evaluate(iH,p);
2018  return HV*Q + (1-Q)*HL;
2019  }
2020  else
2021  {
2022  throw ValueError("TTSE BC error on phase");
2023  }
2024 
2025  // Interval bisection
2026  while (R-L>1)
2027  {
2028  if (isbetween((*mat)[M][j],(*mat)[R][j],Other))
2029  {
2030  L=M; M=(L+R)/2; continue;
2031  }
2032  else
2033  {
2034  R=M; M=(L+R)/2; continue;
2035  }
2036  }
2037  i = L;
2038  }
2039  // Get the nearest full cell to this point
2040  nearest_good_neighbor_ph_interpolate(&i, &j);
2041 
2042  // Invert the bicubic to find h
2043  // Old version was "summer += (*alpha)[ii+jj*4]*x^i*y^j" for i in [0,3] and j in [0,3]
2044  /*
2045  Summation function is
2046  summer = (*alpha)[0+0*4]*x_0*y_0+(*alpha)[0+1*4]*x_0*y_1+(*alpha)[0+2*4]*x_0*y_2+(*alpha)[0+3*4]*x_0*y_3
2047  +(*alpha)[1+0*4]*x_1*y_0+(*alpha)[1+1*4]*x_1*y_1+(*alpha)[1+2*4]*x_1*y_2+(*alpha)[1+3*4]*x_1*y_3
2048  +(*alpha)[2+0*4]*x_2*y_0+(*alpha)[2+1*4]*x_2*y_1+(*alpha)[2+2*4]*x_2*y_2+(*alpha)[2+3*4]*x_2*y_3
2049  +(*alpha)[3+0*4]*x_3*y_0+(*alpha)[3+1*4]*x_3*y_1+(*alpha)[3+2*4]*x_3*y_2+(*alpha)[3+3*4]*x_3*y_3;
2050 
2051  but in this case we know y, which is simply
2052 
2053  y = (p-this->p[j])/(this->p[j+1]-this->p[j]);
2054 
2055  and we solve for x, which is cubic. Convert to a form like
2056 
2057  0 = ax^3 + b*x^2 + c*x + d
2058 
2059  for use in solve_cubic where
2060 
2061  x = (h-this->h[i])/(this->h[i+1]-this->h[i]);
2062 
2063  */
2064  double y = (p-this->p[j])/(this->p[j+1]-this->p[j]);
2065  double y_0 = 1, y_1 = y, y_2 = y*y, y_3 = y*y*y;
2066 
2067  // Find the coefficients for the cubic
2068  std::vector<double> *alpha = NULL;
2069 
2070  // Get bicubic coefficients
2071  alpha = this->bicubic_cell_coeffs_ph(iOther, i, j);
2072 
2073  double a = (*alpha)[3+0*4]*y_0+(*alpha)[3+1*4]*y_1+(*alpha)[3+2*4]*y_2+(*alpha)[3+3*4]*y_3; // factors of x^3
2074  double b = (*alpha)[2+0*4]*y_0+(*alpha)[2+1*4]*y_1+(*alpha)[2+2*4]*y_2+(*alpha)[2+3*4]*y_3; // factors of x^2
2075  double c = (*alpha)[1+0*4]*y_0+(*alpha)[1+1*4]*y_1+(*alpha)[1+2*4]*y_2+(*alpha)[1+3*4]*y_3; // factors of x
2076  double d = (*alpha)[0+0*4]*y_0+(*alpha)[0+1*4]*y_1+(*alpha)[0+2*4]*y_2+(*alpha)[0+3*4]*y_3 - Other; // constant factors
2077  double x0,x1,x2;
2078  solve_cubic(a,b,c,d,&x0,&x1,&x2);
2079 
2080  double x;
2081  // Only one solution
2082  if (fabs(x0-x1) < 10*DBL_EPSILON && fabs(x0-x2) < 10*DBL_EPSILON && fabs(x1-x2) < 10*DBL_EPSILON){
2083  x = x0;
2084  }
2085  // See if first is good solution 0 < x < 1
2086  else if ( 0 <= x0 && x0 <= 1){
2087  x = x0;
2088  }
2089  else if ( 0 <= x1 && x1 <= 1){
2090  x = x1;
2091  }
2092  else if ( 0 <= x2 && x2 <= 1){
2093  x = x2;
2094  }
2095  else
2096  {
2097  throw ValueError("Multiple solutions found for cubic in TTSE-BICUBIC-P+Other");
2098  }
2099  return x*(this->h[i+1]-this->h[i])+this->h[i];
2100  }
2101  // One is enthalpy, we are getting pressure
2102  else if (iInput1 == iH)
2103  {
2104  throw ValueError("Sorry enthalpy and something else other than p is not valid for TTSE");
2105  }
2106  // Oops, neither enthalpy or pressure provided
2107  else
2108  {
2109  throw ValueError("Neither enthalpy nor pressure provided as the first parameter");
2110  }
2111 }
2112 
2113 bool TTSESinglePhaseTableClass::within_range(long iInput1, double Value1, long iInput2, double Value2)
2114 {
2115  int buf = 0;
2116  // get more for bicubic
2117  if (mode == TTSE_MODE_BICUBIC){
2118  buf += 2;
2119  }
2120 
2121  // For p,h as inputs
2122  if (match_pair(iInput1,iInput2,iP,iH)){
2123  // Sort in the order p,h
2124  sort_pair(&iInput1,&Value1,&iInput2,&Value2,iP,iH);
2125  //return (Value1 > pmin && Value1 < pmax && Value2 > hmin && Value2 < hmax);
2126  int i = (int)round((log(Value1)-logpmin)/logpratio);
2127  int j = (int)round((Value2-hmin)/(hmax-hmin)*(Nh-1));
2128  return (buf <= i && i <= (int)Np-(buf+1) && buf <= (int)j && j <= (int)Nh-(buf+1));
2129  }
2130  else if (match_pair(iInput1,iInput2,iP,iT)){
2131  // Sort in the order p, T
2132  sort_pair(&iInput1,&Value1,&iInput2,&Value2,iP,iT);
2133  return this->within_range_one_other_input(iP,Value1,iT,Value2,buf);
2134  }
2135  else if (match_pair(iInput1,iInput2,iP,iD)){
2136  // Sort in the order p, D
2137  sort_pair(&iInput1,&Value1,&iInput2,&Value2, iP, iD);
2138  return this->within_range_one_other_input(iP,Value1,iD,Value2,buf);
2139  }
2140  else if (match_pair(iInput1,iInput2,iP,iS)){
2141  // Sort in the order p, S
2142  sort_pair(&iInput1,&Value1,&iInput2,&Value2, iP, iS);
2143  return this->within_range_one_other_input(iP,Value1,iS,Value2,buf);
2144  }
2145  else if (match_pair(iInput1,iInput2,iT,iD)){
2146  // Sort in the order T, rho
2147  sort_pair(&iInput1,&Value1,&iInput2,&Value2, iT, iD);
2148  //return (Value1 > Tmin && Value1 < Tmax && Value2 > rhomin && Value2 < rhomax);
2149  int i = (int)round((Value1-Tmin)/(Tmax-Tmin)*(NT-1));
2150  int j = (int)round((log(Value2)-logrhomin)/logrhoratio);
2151  return (buf <= i && i <= (int)NT-(buf+1) && buf <= (int)j && j <= (int)Nrho-(buf+1));
2152  }
2153  //throw ValueError("Your input pair was not valid. Please supply either (P,H), (P,T), (P,D), (P,S) or (T,D) as inputs.");
2154  return true;
2155 }
2156 
2157 bool TTSESinglePhaseTableClass::within_range_one_other_input(long iInput1, double Input1, long iOther, double Other, int buf)
2158 {
2159  if (iInput1 != iP)
2160  {
2161  throw ValueError("iInput1 must be iP");
2162  }
2163 
2164  int j = (int)round((log(Input1)-logpmin)/logpratio);
2165  if (!(buf <= j && j <= (int)Np-(buf+1) )) {
2166  return false;
2167  }
2168 
2169  double low, high;
2170  if (iOther == iT)
2171  { high = this->T[Nh-(buf+1)][j]; low = this->T[buf][j]; }
2172  else if(iOther == iS)
2173  { high = this->s[Nh-(buf+1)][j]; low = this->s[buf][j]; }
2174  else if(iOther == iD)
2175  { low = this->rho[Nh-(buf+1)][j]; high = this->rho[buf][j]; }
2176  else
2177  {
2178  throw ValueError("Your input pair was not valid. Please supply either (P,T), (P,S) or (P,D) as inputs.");
2179  return false;
2180  }
2181 
2182  if (ValidNumber(high) && ValidNumber(low))
2183  { return !(low>Other || Other>high); }
2184 
2185  return false;
2186 }
2187 
2188 double TTSESinglePhaseTableClass::evaluate_one_other_input(long iInput1, double Input1, long iOther, double Other)
2189 {
2190  // Use Bicubic interpolation if requested
2191  if (mode == TTSE_MODE_BICUBIC){
2192  return bicubic_evaluate_one_other_input(iInput1,Input1,iOther,Other);
2193  }
2194  // My goal here is to find deltah
2195  int L,R,M,i;
2196  double p,dh1,dh2,a,b,c;
2197  std::vector<std::vector<double> > *mat;
2198 
2199  // Connect a pointer to the array of interest
2200  switch (iOther)
2201  {
2202  case iT:
2203  mat = &T; break;
2204  case iS:
2205  mat = &s; break;
2206  case iD:
2207  mat = &rho; break;
2208  }
2209 
2210  // One is pressure, we are getting enthalpy
2211  // within_TTSE_range() guarantees that the pressure is within range
2212  if (iInput1 == iP)
2213  {
2214  p = Input1;
2215 
2216  int j = (int)round((log(p)-logpmin)/logpratio);
2217  double deltap = p-this->p[j];
2218 
2219  // Check the bounding values for the other input to see if it is within range
2220 
2221  if (j >= jpcrit_ceil) // Is is either supercritical pressure or just a little bit below critical pressure
2222  {
2223  // Very close to the boundary of the LUT, not 1-1 relationship between p-h and other
2224  // sets of inputs, need to allow for a bit of raggedness here
2225  if ( (iOther == iT && Other < 1.1*this->T[Np-1][j] && Other > this->T[Np-1][j])
2226  || (iOther == iS && Other < this->s[Np-1][j]+0.1 && Other > this->s[Np-1][j])
2227  || (iOther == iD && Other > 0.85*this->rho[Np-1][j] && Other < this->rho[Np-1][j])
2228  )
2229  {
2230  i = Np-1;
2231  }
2232  else
2233  {
2234  // Do interval halving over the whole range since
2235  // there can't be any saturation curve
2236  L = 0; R = Nh-1; M = (L+R)/2;
2237  while (R-L>1){
2238  if (isbetween((*mat)[M][j],(*mat)[R][j],Other)){
2239  L=M; M=(L+R)/2; continue;
2240  }
2241  else{
2242  R=M; M=(L+R)/2; continue;
2243  }
2244  }
2245  // Find which one of the bounds is closer
2246  if (fabs((*mat)[L][j]-Other)<fabs((*mat)[R][j]-Other)){
2247  i = L;
2248  }
2249  else{
2250  i = R;
2251  }
2252  }
2253  }
2254  else
2255  {
2256  if ( (iOther == iT && Other > SatV->evaluate(iT,p))
2257  || (iOther == iS && Other > SatV->evaluate(iS,p))
2258  || (iOther == iD && Other < SatV->evaluate(iD,p))
2259  )
2260  {
2261  // SUPERHEATED!!
2262  //
2263  // If it is within between the saturation curve and the first point in the SH region,
2264  // just use the first point in the superheated region
2265  if ( (iOther == iT && Other < this->T[IV[j]+1][j])
2266  || (iOther == iS && Other < this->s[IV[j]+1][j])
2267  || (iOther == iD && Other > this->rho[IV[j]+1][j])
2268  )
2269  {
2270  i = IV[j]+1;
2271  }
2272  // Very close to the boundary of the LUT, not 1-1 relationship between p-h and other
2273  // sets of inputs, need to allow for a bit of raggedness here
2274  else if ( (iOther == iT && Other < 1.1*this->T[Np-1][j] && Other > this->T[Np-1][j])
2275  || (iOther == iS && Other < this->s[Np-1][j]+0.1 && Other > this->s[Np-1][j])
2276  || (iOther == iD && Other > 0.85*this->rho[Np-1][j] && Other < this->rho[Np-1][j])
2277  )
2278  {
2279  i = Np-1;
2280  }
2281  else
2282  {
2283  // SUPERHEATED!! (and away from saturation)
2284  // Make sure it is in the bounds
2285  switch (iOther)
2286  {
2287  case iT:
2288  if (Other > this->T[Nh-1][j])
2289  throw ValueError(format("Input T [%g] is greater than max T [%g] at LUT pressure [%g]",Other,this->T[Np-1][j],this->p[j]));
2290  break;
2291  case iD:
2292  if (Other < this->rho[Nh-1][j])
2293  throw ValueError(format("Input rho [%g] is less than minimum rho [%g] at LUT pressure [%g]",Other,this->rho[Np-1][j],this->p[j]));
2294  break;
2295  case iS:
2296  if (Other > this->s[Nh-1][j])
2297  throw ValueError(format("Input s [%g] is greater than max s [%g] at LUT pressure [%g]",Other,this->s[Np-1][j],this->p[j]));
2298  break;
2299  }
2300 
2301  L = IV[j]+1; R = Np-1; M = (L+R)/2;
2302  while (R-L>1)
2303  {
2304  if (isbetween((*mat)[M][j],(*mat)[R][j],Other))
2305  {
2306  L=M; M=(L+R)/2; continue;
2307  }
2308  else
2309  {
2310  R=M; M=(L+R)/2; continue;
2311  }
2312  }
2313  // Find which one of the bounds is closer
2314  if (fabs((*mat)[L][j]-Other)<fabs((*mat)[R][j]-Other))
2315  i = L;
2316  else
2317  i = R;
2318  }
2319  }
2320  else if (IL[j] < 2)
2321  {
2322  // We are at low pressure, so we don't know how to calculate, going to just use the i==1 element
2323  // if it is valid, or the i = 0 if not, otherwise, there are no values left and we have to fail
2324  if (ValidNumber(this->T[1][j])){
2325  i = 1;
2326  }
2327  else if (ValidNumber(this->T[0][j])){
2328  i = 0;
2329  }
2330  else{
2331  throw ValueError(format("Your inputs [%g,%g] do not yield a valid TTSE node",Input1,Other));
2332  }
2333  }
2334 
2335  else if ( (iOther == iT && Other < SatL->evaluate(iT,p))
2336  || (iOther == iS && Other < SatL->evaluate(iS,p))
2337  || (iOther == iD && Other > SatL->evaluate(iD,p))
2338  )
2339  {
2340  // SUBCOOLED!!
2341  //
2342  // If it is within one spacing of the outlet variable of the saturation curve,
2343  // just use the first point in the subcooled region
2344  if ( (iOther == iT && Other > this->T[IL[j]-1][j])
2345  || (iOther == iS && Other > this->s[IL[j]-1][j])
2346  || (iOther == iD && Other < this->rho[IL[j]-1][j])
2347  )
2348  {
2349  i = IL[j]-1;
2350  }
2351  else{
2352  // Make sure it is in the bounds of the LUT
2353  switch (iOther)
2354  {
2355  case iT:
2356  if (Other < this->T[0][j])
2357  throw ValueError(format("Input T [%g] is less than min T [%g] at LUT pressure [%g]",Other,this->T[0][j],this->p[j]));
2358  break;
2359  case iD:
2360  if (Other > this->rho[0][j])
2361  throw ValueError(format("Input rho [%g] is greater than max rho [%g] at LUT pressure [%g]",Other,this->rho[0][j],this->p[j]));
2362  break;
2363  case iS:
2364  if (Other < this->s[0][j])
2365  throw ValueError(format("Input s [%g] is less than min s [%g] at LUT pressure [%g]",Other,this->s[0][j],this->p[j]));
2366  break;
2367  }
2368 
2369  L = 0; R = IL[j]-1; M = (L+R)/2;
2370  // Its subcooled
2371  while (R-L>1)
2372  {
2373  if (isbetween((*mat)[M][j],(*mat)[R][j],Other))
2374  {
2375  L=M; M=(L+R)/2; continue;
2376  }
2377  else
2378  {
2379  R=M; M=(L+R)/2; continue;
2380  }
2381  }
2382  // Find which one of the bounds is closer
2383  if (fabs((*mat)[L][j]-Other)<fabs((*mat)[R][j]-Other))
2384  i = L;
2385  else
2386  i = R;
2387  }
2388  }
2389  else
2390  {
2391  // It's two-phase
2392  throw ValueError(format("It's two phase input"));
2393  }
2394  }
2395 
2396  // Now we calculate deltah
2397  switch (iOther)
2398  {
2399  // Quadratic in deltah
2400  // 0 = 0.5*deltah*deltah*d2Tdh2[i][j]+deltah*dTdh[i][j]+deltap*deltah*d2Tdhdp[i][j]+T[i][j]-T+deltap*dTdp[i][j]+0.5*deltap*deltap*d2Tdp2[i][j];
2401  // 0 = a*deltah^2+b*deltah+c
2402  case iT:
2403  a = 0.5*d2Tdh2[i][j];
2404  b = dTdh[i][j]+deltap*d2Tdhdp[i][j];
2405  c = T[i][j]-Other+deltap*dTdp[i][j]+0.5*deltap*deltap*d2Tdp2[i][j];
2406  break;
2407  case iS:
2408  a = 0.5*d2sdh2[i][j];
2409  b = dsdh[i][j]+deltap*d2sdhdp[i][j];
2410  c = s[i][j]-Other+deltap*dsdp[i][j]+0.5*deltap*deltap*d2sdp2[i][j];
2411  break;
2412  case iD:
2413  a = 0.5*d2rhodh2[i][j];
2414  b = drhodh[i][j]+deltap*d2rhodhdp[i][j];
2415  c = rho[i][j]-Other+deltap*drhodp[i][j]+0.5*deltap*deltap*d2rhodp2[i][j];
2416  break;
2417  }
2418  // Solutions from quadratic equation
2419  dh1 = (-b+sqrt(b*b-4*a*c))/(2*a);
2420  dh2 = (-b-sqrt(b*b-4*a*c))/(2*a);
2421 
2422  double hspacing = (hmax-hmin)/((double)Nh-1);
2423  // If only one is less than a multiple of enthalpy spacing, thats your solution
2424  if (fabs(dh1) < 10*hspacing && !(fabs(dh2) < 10*hspacing) )
2425  return dh1+this->h[i];
2426  else if (fabs(dh2) < 10*hspacing && !(fabs(dh1) < 10*hspacing) )
2427  return dh2+this->h[i];
2428  else{
2429  // Need to figure out which is the correct solution, try just the smaller one
2430  if (fabs(dh1)<fabs(dh2))
2431  return dh1+this->h[i];
2432  else
2433  return dh2+this->h[i];
2434  }
2435  }
2436  // One is enthalpy, we are getting pressure
2437  else if (iInput1 == iH)
2438  {
2439  throw ValueError("Sorry enthalpy and something else other than p is not valid for TTSE");
2440  }
2441  // Oops, neither enthalpy or pressure provided
2442  else
2443  {
2444  throw ValueError("Neither enthalpy nor pressure provided as the first parameter");
2445  }
2446 }
2447 
2448 double TTSESinglePhaseTableClass::evaluate_Trho(long iOutput, double T, double rho, double logrho)
2449 {
2450  if (!this->within_range(iT,T,iD,rho))
2451  {
2452  throw ValueError(format("Input to TTSE [T = %0.16g, rho = %0.16g] is out of range",T,rho));
2453  }
2454 
2455  // Use Bicubic interpolation if requested
2456  if (this->mode == TTSE_MODE_BICUBIC){
2457  return interpolate_bicubic_Trho(iOutput,T,rho,logrho);
2458  }
2459 
2460  int i = (int)round((T-Tmin)/(Tmax-Tmin)*(NT-1));
2461  int j = (int)round((logrho-logrhomin)/logrhoratio);
2462 
2463  // If the value at i,j is too close to the saturation boundary, the nearest point i,j
2464  // might be in the two-phase region which is not defined for single-phase table.
2465  // Therefore, search around its neighbors for a better choice
2466  if (!ValidNumber(mu_Trho[i][j])){
2467  nearest_good_neighbor_Trho(&i,&j);
2468  }
2469 
2470  // Distances from the node
2471  double deltaT = T - this->T_Trho[i];
2472  double deltarho = rho - this->rho_Trho[j];
2473 
2474  switch (iOutput)
2475  {
2476  case iS:
2477  return s_Trho[i][j]+deltaT*dsdT_Trho[i][j]+deltarho*dsdrho_Trho[i][j]+0.5*deltaT*deltaT*d2sdT2_Trho[i][j]+0.5*deltarho*deltarho*d2sdrho2_Trho[i][j]+deltaT*deltarho*d2sdTdrho_Trho[i][j]; break;
2478  case iP:
2479  return p_Trho[i][j]+deltaT*dpdT_Trho[i][j]+deltarho*dpdrho_Trho[i][j]+0.5*deltaT*deltaT*d2pdT2_Trho[i][j]+0.5*deltarho*deltarho*d2pdrho2_Trho[i][j]+deltaT*deltarho*d2pdTdrho_Trho[i][j]; break;
2480  case iH:
2481  return h_Trho[i][j]+deltaT*dhdT_Trho[i][j]+deltarho*dhdrho_Trho[i][j]+0.5*deltaT*deltaT*d2hdT2_Trho[i][j]+0.5*deltarho*deltarho*d2hdrho2_Trho[i][j]+deltaT*deltarho*d2hdTdrho_Trho[i][j]; break;
2482  case iV:
2483  return mu_Trho[i][j]+deltaT*dmudT_Trho[i][j]+deltarho*dmudrho_Trho[i][j]+0.5*deltaT*deltaT*d2mudT2_Trho[i][j]+0.5*deltarho*deltarho*d2mudrho2_Trho[i][j]+deltaT*deltarho*d2mudTdrho_Trho[i][j]; break;
2484  case iL:
2485  return k_Trho[i][j]+deltaT*dkdT_Trho[i][j]+deltarho*dkdrho_Trho[i][j]+0.5*deltaT*deltaT*d2kdT2_Trho[i][j]+0.5*deltarho*deltarho*d2kdrho2_Trho[i][j]+deltaT*deltarho*d2kdTdrho_Trho[i][j]; break;
2486  default:
2487  throw ValueError(format("Output key value [%d] to evaluate is invalid",iOutput));
2488  }
2489 }
2490 
2491 double TTSESinglePhaseTableClass::evaluate_first_derivative(long iOF, long iWRT, long iCONSTANT, double p, double logp, double h)
2492 {
2493 
2494  if (!this->within_range(iP,p,iH,h))
2495  {
2496  throw ValueError(format("Input to TTSE deriv [p = %0.16g, h = %0.16g] is out of range",p,h));
2497  }
2498 
2499  // Use Bicubic interpolation if requested
2500  if (mode == TTSE_MODE_BICUBIC){
2501  return bicubic_evaluate_first_derivative_ph(iOF, iWRT, iCONSTANT, p, logp, h);
2502  }
2503 
2504  int i = (int)round(((h-hmin)/(hmax-hmin)*(Nh-1)));
2505  int j = (int)round((logp-logpmin)/logpratio);
2506 
2507  // If the value at i,j is too close to the saturation boundary, the nearest point i,j
2508  // might be in the two-phase region which is not defined for single-phase table.
2509  // Therefore, search around its neighbors for a better choice
2510  if (!ValidNumber(T[i][j]))
2511  {
2512  nearest_good_neighbor(&i,&j);
2513  }
2514 
2515  // Distances from the node
2516  double deltah = h-this->h[i];
2517  double deltap = p-this->p[j];
2518 
2519  // This is a first-order expansion of the derivative around the node point.
2520  //
2521  // Derivatives for constant p
2522  if (iOF == iT && iWRT == iH && iCONSTANT == iP)
2523  {
2524  // Derivative of T w.r.t. h for p constant (for cp, the constant-pressure specific heat)
2525  return dTdh[i][j]+deltah*d2Tdh2[i][j]+deltap*d2Tdhdp[i][j];
2526  }
2527  else if (iOF == iS && iWRT == iH && iCONSTANT == iP)
2528  {
2529  // Derivative of s w.r.t. h for p constant
2530  return dsdh[i][j]+deltah*d2sdh2[i][j]+deltap*d2sdhdp[i][j];
2531  }
2532  else if (iOF == iD && iWRT == iH && iCONSTANT == iP)
2533  {
2534  // Derivative of density w.r.t. h for p constant
2535  return drhodh[i][j]+deltah*d2rhodh2[i][j]+deltap*d2rhodhdp[i][j];
2536  }
2537 
2538  // Derivatives for constant h
2539  else if (iOF == iT && iWRT == iP && iCONSTANT == iH)
2540  {
2541  // Derivative of T w.r.t. p for h constant
2542  return dTdp[i][j]+deltap*d2Tdp2[i][j]+deltah*d2Tdhdp[i][j];
2543  }
2544  else if (iOF == iS && iWRT == iP && iCONSTANT == iH)
2545  {
2546  // Derivative of s w.r.t. p for h constant
2547  return dsdp[i][j]+deltap*d2sdp2[i][j]+deltah*d2sdhdp[i][j];
2548  }
2549  else if (iOF == iD && iWRT == iP && iCONSTANT == iH)
2550  {
2551  // Derivative of density w.r.t. p for h constant
2552  return drhodp[i][j]+deltap*d2rhodp2[i][j]+deltah*d2rhodhdp[i][j];
2553  }
2554  else{
2555  throw ValueError(format("Your inputs [%d,%d,%d,%g,%g] are invalid to evaluate_first_derivative",iOF,iWRT,iCONSTANT,p,h));
2556  }
2557 }
2558 
2559 TTSETwoPhaseTableClass::TTSETwoPhaseTableClass(Fluid *pFluid, double Q)
2560 {
2561  this->pFluid = pFluid;
2562  this->Q = Q;
2563 }
2564 void TTSETwoPhaseTableClass::set_size(unsigned int N)
2565 {
2566  this->N = N;
2567 
2568  // Seed the generator
2569  srand((unsigned int)time(NULL));
2570 
2571  // Resize all the arrays
2572  h.resize(N);
2573  p.resize(N);
2574  logp.resize(N);
2575  T.resize(N);
2576  dTdp.resize(N);
2577  d2Tdp2.resize(N);
2578  rho.resize(N);
2579  logrho.resize(N);
2580  drhodp.resize(N);
2581  d2rhodp2.resize(N);
2582  s.resize(N);
2583  dsdp.resize(N);
2584  d2sdp2.resize(N);
2585  h.resize(N);
2586  dhdp.resize(N);
2587  d2hdp2.resize(N);
2588 }
2589 
2590 double TTSETwoPhaseTableClass::build(double pmin, double pmax, TTSETwoPhaseTableClass *other)
2591 {
2592  CoolPropStateClassSI CPS(pFluid);
2593 
2594  this->pmin = pmin;
2595  this->pmax = pmax;
2596  this->logpmin = log(pmin);
2597  this->logpmax = log(pmax);
2598  this->pratio = pow(pmax/pmin,1/((double)N-1));
2599  this->logpratio = log(pratio);
2600 
2601  double dlogp = (logpmax-logpmin)/(N-1);
2602  clock_t t1,t2;
2603  t1 = clock();
2604  // Logarithmic distribution of pressures
2605  for (unsigned int i = 0; i < N-1; i++)
2606  {
2607  // Calculate the pressure
2608  p[i] = exp(logpmin + i*dlogp);
2609  logp[i] = log(p[i]);
2610  // Update the class
2611  CPS.update(iP,p[i],iQ,Q);
2612 
2613  // Set the variables
2614  T[i] = CPS.T();
2615  dTdp[i] = CPS.dTdp_along_sat();
2616  d2Tdp2[i] = CPS.d2Tdp2_along_sat();
2617  h[i] = CPS.h();
2618  dhdp[i] = (this->Q>0.5) ? CPS.dhdp_along_sat_vapor() : CPS.dhdp_along_sat_liquid();
2619  d2hdp2[i] = (this->Q>0.5) ? CPS.d2hdp2_along_sat_vapor() : CPS.d2hdp2_along_sat_liquid();
2620  s[i] = CPS.s();
2621  dsdp[i] = (this->Q>0.5) ? CPS.dsdp_along_sat_vapor() : CPS.dsdp_along_sat_liquid();
2622  d2sdp2[i] = (this->Q>0.5) ? CPS.d2sdp2_along_sat_vapor() : CPS.d2sdp2_along_sat_liquid();
2623  rho[i] = CPS.rho();
2624  logrho[i] = log(CPS.rho());
2625  drhodp[i] = (this->Q>0.5) ? CPS.drhodp_along_sat_vapor() : CPS.drhodp_along_sat_liquid();
2626  d2rhodp2[i] = (this->Q>0.5) ? CPS.d2rhodp2_along_sat_vapor() : CPS.d2rhodp2_along_sat_liquid();
2627 
2628  // If other is provided
2629  if (other != NULL)
2630  {
2631  other->pmin = pmin;
2632  other->pmax = pmax;
2633  other->logpmin = log(pmin);
2634  other->logpmax = log(pmax);
2635 
2636  other->p[i] = this->p[i];
2637  other->logp[i] = log(this->p[i]);
2638 
2639  // Set the variables
2640  other->T[i] = CPS.T();
2641  other->dTdp[i] = CPS.dTdp_along_sat();
2642  other->d2Tdp2[i] = CPS.d2Tdp2_along_sat();
2643  other->h[i] = (other->Q>0.5) ? CPS.hV() : CPS.hL();
2644  other->dhdp[i] = (other->Q>0.5) ? CPS.dhdp_along_sat_vapor() : CPS.dhdp_along_sat_liquid();
2645  other->d2hdp2[i] = (other->Q>0.5) ? CPS.d2hdp2_along_sat_vapor() : CPS.d2hdp2_along_sat_liquid();
2646  other->s[i] = (other->Q>0.5) ? CPS.sV() : CPS.sL();
2647  other->dsdp[i] = (other->Q>0.5) ? CPS.dsdp_along_sat_vapor() : CPS.dsdp_along_sat_liquid();
2648  other->d2sdp2[i] = (other->Q>0.5) ? CPS.d2sdp2_along_sat_vapor() : CPS.d2sdp2_along_sat_liquid();
2649  other->rho[i] = (other->Q>0.5) ? CPS.rhoV() : CPS.rhoL();
2650  other->logrho[i] = log(other->rho[i]);
2651  other->drhodp[i] = (other->Q>0.5) ? CPS.drhodp_along_sat_vapor() : CPS.drhodp_along_sat_liquid();
2652  other->d2rhodp2[i] = (other->Q>0.5) ? CPS.d2rhodp2_along_sat_vapor() : CPS.d2rhodp2_along_sat_liquid();
2653  }
2654  }
2655  // At the last point (the critical point)
2656  CPS.flag_SinglePhase = true; // Don't have it check the state or do a saturation call
2657  CPS.update(iT,CPS.pFluid->reduce.T+1e-12,iD,CPS.pFluid->reduce.rho+1e-12);
2658  p[N-1] = CPS.p();
2659  logp[N-1] = log(p[N-1]);
2660  T[N-1] = CPS.T();
2661  dTdp[N-1] = 1/CPS.dpdT_constrho();
2662  //d2Tdp2[N-1] = CPS.d2Tdp2_along_sat();
2663  h[N-1] = CPS.h();
2664  //dhdp[N-1] = (this->Q>0.5) ? CPS.dhdp_along_sat_vapor() : CPS.dhdp_along_sat_liquid();
2665  //d2hdp2[N-1] = (this->Q>0.5) ? CPS.d2hdp2_along_sat_vapor() : CPS.d2hdp2_along_sat_liquid();
2666  s[N-1] = CPS.s();
2667  //dsdp[N-1] = (this->Q>0.5) ? CPS.dsdp_along_sat_vapor() : CPS.dsdp_along_sat_liquid();
2668  //d2sdp2[N-1] = (this->Q>0.5) ? CPS.d2sdp2_along_sat_vapor() : CPS.d2sdp2_along_sat_liquid();
2669  rho[N-1] = CPS.rho();
2670  logrho[N-1] = log(CPS.rho());
2671  //drhodp[N-1] = (this->Q>0.5) ? CPS.drhodp_along_sat_vapor() : CPS.drhodp_along_sat_liquid();
2672  //d2rhodp2[N-1] = (this->Q>0.5) ? CPS.d2rhodp2_along_sat_vapor() : CPS.d2rhodp2_along_sat_liquid();
2673 
2674  // If other is provided
2675  if (other != NULL)
2676  {
2677  other->p[N-1] = CPS.p();
2678  other->logp[N-1] = log(p[N-1]);
2679  other->T[N-1] = CPS.T();
2680  other->dTdp[N-1] = 1/CPS.dpdT_constrho();
2681  //other->d2Tdp2[N-1] = CPS.d2Tdp2_along_sat();
2682  other->h[N-1] = CPS.h();
2683  //other->dhdp[N-1] = (this->Q>0.5) ? CPS.dhdp_along_sat_vapor() : CPS.dhdp_along_sat_liquid();
2684  //other->d2hdp2[N-1] = (this->Q>0.5) ? CPS.d2hdp2_along_sat_vapor() : CPS.d2hdp2_along_sat_liquid();
2685  other->s[N-1] = CPS.s();
2686  //other->dsdp[N-1] = (this->Q>0.5) ? CPS.dsdp_along_sat_vapor() : CPS.dsdp_along_sat_liquid();
2687  //other->d2sdp2[N-1] = (this->Q>0.5) ? CPS.d2sdp2_along_sat_vapor() : CPS.d2sdp2_along_sat_liquid();
2688  other->rho[N-1] = CPS.rho();
2689  other->logrho[N-1] = log(CPS.rho());
2690  //other->drhodp[N-1] = (this->Q>0.5) ? CPS.drhodp_along_sat_vapor() : CPS.drhodp_along_sat_liquid();
2691  //other->d2rhodp2[N-1] = (this->Q>0.5) ? CPS.d2rhodp2_along_sat_vapor() : CPS.d2rhodp2_along_sat_liquid();
2692  }
2693 
2694  t2 = clock();
2695  std::cout << double(t2-t1)/CLOCKS_PER_SEC << " to build both two phase tables" << std::endl;
2696  return double(t2-t1)/CLOCKS_PER_SEC;
2697 }
2698 
2699 double TTSETwoPhaseTableClass::evaluate(long iParam, double p)
2700 {
2701  double logp = log(p);
2702  int i = (int)round(((logp-logpmin)/(logpmax-logpmin)*(N-1)));
2703  // If the value is just a little bit below the range, clip
2704  // it back to the range of the LUT
2705  if (i == -1) i = 0;
2706  // If the pressure is just barely above the critical pressure
2707  // or between the critical pressure and the next lowest point,
2708  // just just the next lowest point to avoid some of the
2709  // derivatives that are infinite at the critical point
2710  if (i == (int)N || i == (int)N-1 ) i = N-2;
2711  // If it is really out of the range, throw an error
2712  if (i<0 || i>(int)N-1)
2713  {
2714  throw ValueError(format("p [%g] is out of range[%g,%g], yielded index of: %d",p,pmin,pmax,i));
2715  }
2716 
2717  if (i == (int)N-2)
2718  {
2719  double y1, yc, dydp1;
2720  switch (iParam)
2721  {
2722  case iS:
2723  y1 = s[i]; yc = s[i+1]; dydp1 = dsdp[i]; break;
2724  case iT:
2725  y1 = T[i]; yc = T[i+1]; dydp1 = dTdp[i]; break;
2726  case iH:
2727  y1 = h[i]; yc = h[i+1]; dydp1 = dhdp[i]; break;
2728  case iD:
2729  y1 = rho[i]; yc = rho[i+1]; dydp1 = drhodp[i]; break;
2730  default:
2731  throw ValueError();
2732  }
2733 
2734  // Here we use interpolation based on a form y = a*x^(1/n)
2735  // where we develop shifted coordinates
2736  // Y = y - yc (where y is s, h, rho, etc.)
2737  // X = pc - p
2738 
2739  double X = this->p[i+1] - p;
2740  double Y1 = y1 - yc;
2741  double X1 = this->p[i+1] - this->p[i];
2742  double dYdX1 = -dydp1;
2743  double n = Y1/(dYdX1*X1);
2744  double a = Y1/pow(X1,1/n);
2745  double Y = a*pow(X,1.0/n);
2746  return Y + yc;
2747  }
2748  else
2749  {
2750  // Spline interpolation http://en.wikipedia.org/wiki/Spline_interpolation since we
2751  // know the derivatives and the values at the bounding elements
2752  // Independent variable is logp
2753  // Dependent variable is varied (entropy, enthalpy, etc...)
2754  double x1 = this->logp[i];
2755  double x2 = this->logp[i+1];
2756  double t = (logp-x1)/(x2-x1);
2757 
2758  double y1, y2, k1, k2;
2759  switch (iParam)
2760  {
2761  case iS:
2762  y1 = this->s[i]; y2 = this->s[i+1]; k1 = this->p[i]*this->dsdp[i]; k2 = this->p[i+1]*this->dsdp[i+1]; break;
2763  case iT:
2764  y1 = this->T[i]; y2 = this->T[i+1]; k1 = this->p[i]*this->dTdp[i]; k2 = this->p[i+1]*this->dTdp[i+1]; break;
2765  case iH:
2766  y1 = this->h[i]; y2 = this->h[i+1]; k1 = this->p[i]*this->dhdp[i]; k2 = this->p[i+1]*this->dhdp[i+1]; break;
2767  case iD:
2770  y1 = this->logrho[i]; y2 = this->logrho[i+1]; k1 = this->p[i]/this->rho[i]*this->drhodp[i]; k2 = this->p[i+1]/this->rho[i+1]*this->drhodp[i+1]; break;
2771  default:
2772  throw ValueError();
2773  }
2774 
2775  double a = k1*(x2-x1)-(y2-y1);
2776  double b = -k2*(x2-x1)+(y2-y1);
2777  double y = (1-t)*y1+t*y2+t*(1-t)*(a*(1-t)+b*t);
2778  if (iParam == iD){
2779  return exp(y);
2780  }
2781  else{
2782  return y;
2783  }
2784  }
2785 }
2786 double TTSETwoPhaseTableClass::evaluate_T(double T)
2787 {
2788  int L,R,M;
2789  double logp_spacing;
2790 
2791  logp_spacing = this->logp[2]-this->logp[1];
2792 
2793  // Do interval halving over the whole range to find the nearest temperature
2794  L = 0; R = N - 2; M = (L+R)/2;
2795  if (isbetween(this->T[N-2],pFluid->reduce.T,T))
2796  {
2797  // According to Matthis Thorade, dTdP|sat at the critical point is equal to dT/dP|rho evaluated at the
2798  // critical temperature and density
2799  L = N-2;
2800  // Spline interpolation http://en.wikipedia.org/wiki/Spline_interpolation since we
2801  // know the derivatives and the values at the bounding elements
2802  // Independent variable is T
2803  // Dependent variable is logp
2804  double t = (T-this->T[L])/(this->T[L+1]-this->T[L]);
2805  double x1 = this->T[L];
2806  double x2 = this->T[L+1];
2807  double y1 = this->logp[L];
2808  double y2 = this->logp[L+1];
2809  // y is log(p); d(log(p))/dT = 1/p*(dp/dT) = 1/p/(dT/dp)
2810  double k1 = 1/this->p[L]/this->dTdp[L];
2811  double k2 = 1/this->p[L+1]/this->dTdp[L+1];
2812  double a = k1*(x2-x1)-(y2-y1);
2813  double b = -k2*(x2-x1)+(y2-y1);
2814  double logp = (1-t)*y1+t*y2+t*(1-t)*(a*(1-t)+b*t);
2815  return exp(logp);
2816  }
2817  else
2818  {
2819  while (R - L > 1)
2820  {
2821  if (T > this->T[M]){
2822  L=M; M=(L+R)/2; continue;
2823  }
2824  else{
2825  R=M; M=(L+R)/2; continue;
2826  }
2827  }
2828  // Spline interpolation http://en.wikipedia.org/wiki/Spline_interpolation since we
2829  // know the derivatives and the values at the bounding elements
2830  // Independent variable is T
2831  // Dependent variable is logp
2832  double t = (T-this->T[L])/(this->T[R]-this->T[L]);
2833  double x1 = this->T[L];
2834  double x2 = this->T[R];
2835  double y1 = this->logp[L];
2836  double y2 = this->logp[R];
2837  // y is log(p); d(log(p))/dT = 1/p*(dp/dT) = 1/p/(dT/dp)
2838  double k1 = 1/this->p[L]/this->dTdp[L];
2839  double k2 = 1/this->p[R]/this->dTdp[R];
2840  double a = k1*(x2-x1)-(y2-y1);
2841  double b = -k2*(x2-x1)+(y2-y1);
2842  double logp = (1-t)*y1+t*y2+t*(1-t)*(a*(1-t)+b*t);
2843  return exp(logp);
2844  }
2845 
2846 }
2847 double TTSETwoPhaseTableClass::evaluate_sat_derivative(long iParam, double p)
2848 {
2849  double logp = log(p);
2850  int i = (int)round(((logp-logpmin)/(logpmax-logpmin)*(N-1)));
2851  // If the value is just a little bit out of the range, clip
2852  // it back to the range of the LUT
2853  if (i == -1) i = 0;
2854  if (i == (int)N) i = N-1;
2855  // If it is really out of the range, throw an error
2856  if (i<0 || i>(int)N-1)
2857  {
2858  throw ValueError(format("p [%g] is out of range",p));
2859  }
2860 
2861  double pi = this->p[i];
2862 
2863  switch (iParam)
2864  {
2865  case iT:
2866  {
2868  return dTdp[i]+(p-pi)*d2Tdp2[i];
2869  }
2870  case iH:
2871  {
2873  return dhdp[i]+(p-pi)*d2hdp2[i];
2874  }
2875  case iS:
2876  {
2878  return dsdp[i]+(p-pi)*d2sdp2[i];
2879  }
2880  case iD:
2881  {
2883  return drhodp[i]+(p-pi)*d2rhodp2[i];
2884  }
2885  default:
2886  throw ValueError(format("Cannot use the key [%d] provided in evaluate_sat_derivative",iParam));
2887  }
2888 }
2889 
2890 double TTSETwoPhaseTableClass::evaluate_randomly(long iParam, unsigned int N)
2891 {
2892  clock_t t1,t2;
2893  t1 = clock();
2894  for (unsigned int i = 0; i < N; i++)
2895  {
2896  double p1 = ((double)rand()/(double)RAND_MAX)*(pmax-pmin)+pmin;
2897 
2898  // Get the value from TTSE
2899  evaluate(iParam,p1);
2900  }
2901  t2 = clock();
2902  return (double)(t2-t1)/CLOCKS_PER_SEC/(double)N*1e6;
2903 }
2904 
2905 
2906 double TTSETwoPhaseTableClass::check_randomly(long iParam, unsigned int N, std::vector<double> *p, std::vector<double> *EOSv, std::vector<double> *TTSE)
2907 {
2908  double val=0;
2909  p->resize(N);
2910  EOSv->resize(N);
2911  TTSE->resize(N);
2912 
2913  CoolPropStateClassSI CPS(pFluid);
2914 
2915  for (unsigned int i = 0; i < N; i++)
2916  {
2917  double p1 = ((double)rand()/(double)RAND_MAX)*(pmax-pmin)+pmin;
2918 
2919  CPS.update(iP,p1,iQ,this->Q);
2920  double hEOS = CPS.h();
2921  double sEOS = CPS.s();
2922  double TEOS = CPS.T();
2923  double rhoEOS = CPS.rho();
2924 
2925  // Store the inputs
2926  (*p)[i] = p1;
2927 
2928  // Get the value from TTSE
2929  (*TTSE)[i] = evaluate(iParam,p1);
2930 
2931  // Get the value from EOS
2932  switch (iParam)
2933  {
2934  case iS:
2935  (*EOSv)[i] = sEOS; break;
2936  case iT:
2937  (*EOSv)[i] = TEOS; break;
2938  case iH:
2939  (*EOSv)[i] = hEOS; break;
2940  case iD:
2941  (*EOSv)[i] = rhoEOS; break;
2942  default:
2943  throw ValueError();
2944  }
2945 
2946  std::cout << format("%g %g %g %g TTSE (p,EOS,TTSE, diff)\n",p1,(*EOSv)[i],(*TTSE)[i],((*EOSv)[i]-(*TTSE)[i])).c_str();
2947  }
2948  return val;
2949 }
2950 
TTSESinglePhaseTableClass::d2Tdhdp
std::vector< std::vector< double > > d2Tdhdp
Definition: TTSE.h:118
TTSETwoPhaseTableClass::dTdp
std::vector< double > dTdp
Definition: TTSE.h:48
CoolPropStateClassSI::drhodp_along_sat_vapor
double drhodp_along_sat_vapor(void)
Definition: CPState.cpp:2468
TTSESinglePhaseTableClass::s
std::vector< std::vector< double > > s
Definition: TTSE.h:120
TTSESinglePhaseTableClass::root_path
std::string root_path
Definition: TTSE.h:105
make_dirs
void make_dirs(std::string file_path)
Definition: TTSE.cpp:555
TTSETwoPhaseTableClass::d2hdp2
std::vector< double > d2hdp2
Definition: TTSE.h:51
TTSESinglePhaseTableClass::nearest_neighbor_ph
void nearest_neighbor_ph(int i, int j, double *T0, double *rho0)
Definition: TTSE.cpp:425
TTSESinglePhaseTableClass::p
std::vector< double > p
Definition: TTSE.h:121
TTSESinglePhaseTableClass::d2sdh2
std::vector< std::vector< double > > d2sdh2
Definition: TTSE.h:120
TTSESinglePhaseTableClass::bicubic_cell_coeffs_Trho
std::vector< double > * bicubic_cell_coeffs_Trho(long iParam, int i, int j)
Definition: TTSE.cpp:1660
CoolProp.h
TTSESinglePhaseTableClass::nearest_good_neighbor_Trho_interpolate
void nearest_good_neighbor_Trho_interpolate(int *i, int *j)
Definition: TTSE.cpp:306
TTSETwoPhaseTableClass::rho
std::vector< double > rho
Definition: TTSE.h:49
iQ
Definition: GlobalConstants.h:13
TTSESinglePhaseTableClass::interpolate_bicubic_ph
double interpolate_bicubic_ph(long iParam, double p, double logp, double h)
Definition: TTSE.cpp:1825
TTSETwoPhaseTableClass::logpratio
double logpratio
Definition: TTSE.h:30
match_pair
bool match_pair(long iI1, long iI2, long I1, long I2)
Definition: CPState.cpp:132
TTSESinglePhaseTableClass::nearest_good_neighbor_Trho
void nearest_good_neighbor_Trho(int *i, int *j)
Definition: TTSE.cpp:396
TTSESinglePhaseTableClass::evaluate
double evaluate(long iParam, double p, double logp, double h)
Definition: TTSE.cpp:1487
isbetween
bool isbetween(double x1, double x2, double x)
Definition: TTSE.cpp:1928
CPState.h
CoolPropStateClassSI::hV
double hV(void)
Definition: CPState.cpp:1341
Fluid::temperature_ph
virtual void temperature_ph(double p, double h, double &Tout, double &rhoout, double &rhoL, double &rhoV, double &TsatLout, double &TsatVout, double T0=-1, double rho0=-1)
Definition: FluidClass.cpp:1713
get_Fluid_index
long get_Fluid_index(std::string FluidName)
Definition: CoolProp.cpp:239
TTSESinglePhaseTableClass::d2sdTdrho_Trho
std::vector< std::vector< double > > d2sdTdrho_Trho
Definition: TTSE.h:124
TTSETwoPhaseTableClass::build
double build(double pmin, double pmax, TTSETwoPhaseTableClass *other=NULL)
Definition: TTSE.cpp:2590
CoolPropStateClassSI::rhoL
double rhoL(void)
Definition: CPState.h:265
TTSESinglePhaseTableClass::h_Trho
std::vector< std::vector< double > > h_Trho
Definition: TTSE.h:126
Fluid::limits
struct FluidLimits limits
Definition: FluidClass.h:219
TTSETwoPhaseTableClass::dhdp
std::vector< double > dhdp
Definition: TTSE.h:51
TTSESinglePhaseTableClass::IL
std::vector< int > IL
Definition: TTSE.h:133
CoolPropStateClassSI::dhdp_along_sat_vapor
double dhdp_along_sat_vapor(void)
Definition: CPState.cpp:2410
CoolPropStateClassSI::dsdrho_constT
double dsdrho_constT(void)
Definition: CPState.cpp:2226
iC
Definition: GlobalConstants.h:13
TTSETwoPhaseTableClass::check_randomly
double check_randomly(long iParam, unsigned int N, std::vector< double > *p, std::vector< double > *EOSv, std::vector< double > *TTSE)
Definition: TTSE.cpp:2906
TTSESinglePhaseTableClass::dhdrho_Trho
std::vector< std::vector< double > > dhdrho_Trho
Definition: TTSE.h:126
CoolPropStateClassSI::dhdT_constp
double dhdT_constp(void)
Definition: CPState.cpp:2271
TTSESinglePhaseTableClass::nearest_good_neighbor_ph_interpolate
void nearest_good_neighbor_ph_interpolate(int *i, int *j)
Definition: TTSE.cpp:351
CoolPropStateClassSI::d2sdrho2_constT
double d2sdrho2_constT(void)
Definition: CPState.cpp:2238
CoolPropStateClassSI::d2rhodp2_along_sat_liquid
double d2rhodp2_along_sat_liquid(void)
Definition: CPState.cpp:2504
CriticalStruct::p
PressureUnit p
Definition: FluidClass.h:50
TTSETwoPhaseTableClass::dsdp
std::vector< double > dsdp
Definition: TTSE.h:50
TTSESinglePhaseTableClass::rho
std::vector< std::vector< double > > rho
Definition: TTSE.h:119
TTSESinglePhaseTableClass::T
std::vector< std::vector< double > > T
Definition: TTSE.h:118
CoolPropStateClassSI::d2hdp2_along_sat_vapor
double d2hdp2_along_sat_vapor(void)
Definition: CPState.cpp:2422
PressureUnit::Pa
double Pa
Definition: Units.h:22
TTSETwoPhaseTableClass::d2Tdp2
std::vector< double > d2Tdp2
Definition: TTSE.h:48
TTSESinglePhaseTableClass::DL
std::vector< double > DL
Definition: TTSE.h:106
TTSESinglePhaseTableClass::IV
std::vector< int > IV
Definition: TTSE.h:133
CoolPropStateClassSI::d2sdrhodT
double d2sdrhodT(void)
Definition: CPState.cpp:2242
TTSESinglePhaseTableClass::p_Trho
std::vector< std::vector< double > > p_Trho
Definition: TTSE.h:125
CoolPropStateClassSI::dpdT_consth
double dpdT_consth(void)
Definition: CPState.cpp:2258
TTSE_MODE_TTSE
Definition: TTSE.h:8
TTSESinglePhaseTableClass::bicubic_cells
BiCubicCellsContainerClass bicubic_cells
Definition: TTSE.h:115
CoolPropStateClassSI::update
void update(long iInput1, double Value1, long iInput2, double Value2, double T0=-1, double rho0=-1)
Definition: CPState.cpp:213
CoolPropStateClassSI::d2hdp2_along_sat_liquid
double d2hdp2_along_sat_liquid(void)
Definition: CPState.cpp:2430
get_home_dir
std::string get_home_dir(void)
Definition: TTSE.cpp:44
TTSESinglePhaseTableClass::bicubic_evaluate_one_other_input
double bicubic_evaluate_one_other_input(long iInput1, double Input1, long iOther, double Other)
Definition: TTSE.cpp:1933
CoolPropStateClassSI::dhdrho_constp
double dhdrho_constp(void)
Definition: CPState.cpp:2275
TTSESinglePhaseTableClass::k_Trho
std::vector< std::vector< double > > k_Trho
Definition: TTSE.h:128
TTSESinglePhaseTableClass::set_size_Trho
void set_size_Trho(unsigned int NT=100, unsigned int Nrho=100)
Set the sizes of the matrices with Trho as inputs.
Definition: TTSE.cpp:214
strsplit
std::vector< std::string > strsplit(std::string s, char del)
Definition: CoolPropTools.cpp:41
TTSESinglePhaseTableClass::TV
std::vector< double > TV
Definition: TTSE.h:106
TTSESinglePhaseTableClass::dkdrho_Trho
std::vector< std::vector< double > > dkdrho_Trho
Definition: TTSE.h:128
CoolPropStateClassSI::d2rhodp2_along_sat_vapor
double d2rhodp2_along_sat_vapor(void)
Definition: CPState.cpp:2496
TTSESinglePhaseTableClass::dpdT_Trho
std::vector< std::vector< double > > dpdT_Trho
Definition: TTSE.h:125
TTSE.h
TTSESinglePhaseTableClass::SatL
TTSETwoPhaseTableClass * SatL
Definition: TTSE.h:103
DBL_EPSILON
#define DBL_EPSILON
Definition: Helmholtz.cpp:10
Fluid::reduce
struct CriticalStruct reduce
A pointer to the point that is used to reduce the T and rho for EOS.
Definition: FluidClass.h:222
TTSESinglePhaseTableClass::interpolate_bicubic_Trho
double interpolate_bicubic_Trho(long iParam, double T, double rho, double logrho)
Definition: TTSE.cpp:1898
CoolPropStateClassSI::sL
double sL(void)
Definition: CPState.cpp:1352
TTSESinglePhaseTableClass::dTdp
std::vector< std::vector< double > > dTdp
Definition: TTSE.h:118
Fluid::get_name
std::string get_name()
Returns a std::string with the name of the fluid.
Definition: FluidClass.h:235
TTSESinglePhaseTableClass::Np
unsigned int Np
Definition: TTSE.h:94
TTSETwoPhaseTableClass::N
unsigned int N
Definition: TTSE.h:27
OtherParameters::Ttriple
double Ttriple
Definition: FluidClass.h:43
round
double round(double r)
Definition: TTSE.cpp:104
TTSESinglePhaseTableClass::dmudrho_Trho
std::vector< std::vector< double > > dmudrho_Trho
Definition: TTSE.h:127
format
std::string format(const char *fmt,...)
Definition: CoolPropTools.cpp:22
TTSESinglePhaseTableClass::bicubic_cell_coordinates_Trho
void bicubic_cell_coordinates_Trho(double Tval, double rho, double logrhoval, int *i, int *j)
Definition: TTSE.cpp:1526
CoolPropStateClassSI::cp
double cp(void)
Definition: CPState.cpp:1460
TTSESinglePhaseTableClass::SV
std::vector< double > SV
Definition: TTSE.h:106
CoolPropStateClassSI::pFluid
Fluid * pFluid
A pointer to a CoolProp fluid.
Definition: CPState.h:195
TTSESinglePhaseTableClass::mu_Trho
std::vector< std::vector< double > > mu_Trho
Definition: TTSE.h:127
TTSESinglePhaseTableClass::evaluate_first_derivative
double evaluate_first_derivative(long iOF, long iWRT, long iCONSTANT, double p, double logp, double h)
Definition: TTSE.cpp:2491
Fluid::TTSESatV
TTSETwoPhaseTableClass TTSESatV
Definition: FluidClass.h:227
Fluid
Fluid is the abstract base class that is employed by all the other fluids.
Definition: FluidClass.h:147
sort_pair
void sort_pair(long *iInput1, double *Value1, long *iInput2, double *Value2, long I1, long I2)
Definition: CPState.cpp:136
TTSESinglePhaseTableClass::bicubic_cell_coordinates_ph
void bicubic_cell_coordinates_ph(double hval, double p, double logpval, int *i, int *j)
Definition: TTSE.cpp:1548
TTSESinglePhaseTableClass::dkdT_Trho
std::vector< std::vector< double > > dkdT_Trho
Definition: TTSE.h:128
TTSESinglePhaseTableClass::d2rhodp2
std::vector< std::vector< double > > d2rhodp2
Definition: TTSE.h:119
TTSESinglePhaseTableClass::d2sdT2_Trho
std::vector< std::vector< double > > d2sdT2_Trho
Definition: TTSE.h:124
iH
Definition: GlobalConstants.h:13
TTSESinglePhaseTableClass::dmudT_Trho
std::vector< std::vector< double > > dmudT_Trho
Definition: TTSE.h:127
CoolPropStateClassSI::d2hdrhodT
double d2hdrhodT(void)
Definition: CPState.cpp:2215
TTSETwoPhaseTableClass::s
std::vector< double > s
Definition: TTSE.h:50
TTSESinglePhaseTableClass::within_range
bool within_range(long iInput1, double Input1, long iInput2, double Input2)
See if the inputs are within range.
Definition: TTSE.cpp:2113
BiCubicCellsContainerClass::cells
std::vector< std::vector< BiCubicCellClass > > cells
Definition: TTSE.h:21
CoolPropStateClassSI::dhdrho_constT
double dhdrho_constT(void)
Definition: CPState.cpp:2203
TTSESinglePhaseTableClass::dsdrho_Trho
std::vector< std::vector< double > > dsdrho_Trho
Definition: TTSE.h:124
iS
Definition: GlobalConstants.h:13
TTSESinglePhaseTableClass::dhdT_Trho
std::vector< std::vector< double > > dhdT_Trho
Definition: TTSE.h:126
CoolPropStateClassSI::dpdrho_constT
double dpdrho_constT(void)
Definition: CPState.cpp:2114
TTSESinglePhaseTableClass::dsdp
std::vector< std::vector< double > > dsdp
Definition: TTSE.h:120
TTSESinglePhaseTableClass::check_randomly
double check_randomly(long iParam, unsigned int N, std::vector< double > *h, std::vector< double > *p, std::vector< double > *EOSv, std::vector< double > *TTSE)
Definition: TTSE.cpp:1420
TTSESinglePhaseTableClass::build_ph
double build_ph(double hmin, double hmax, double pmin, double pmax, TTSETwoPhaseTableClass *SatL=NULL, TTSETwoPhaseTableClass *SatV=NULL)
Definition: TTSE.cpp:834
x
std::vector< double > x(ncmax, 0)
TTSESinglePhaseTableClass::drhodh
std::vector< std::vector< double > > drhodh
Definition: TTSE.h:119
CoolPropStateClassSI::dhdp_along_sat_liquid
double dhdp_along_sat_liquid(void)
Definition: CPState.cpp:2398
TTSESinglePhaseTableClass::h
std::vector< double > h
Definition: TTSE.h:121
TTSETwoPhaseTableClass::TTSETwoPhaseTableClass
TTSETwoPhaseTableClass()
Default Instantiator.
Definition: TTSE.h:33
TTSESinglePhaseTableClass::set_mode
int set_mode(int mode)
Definition: TTSE.cpp:160
TTSESinglePhaseTableClass::write_dotdrawing_tofile
void write_dotdrawing_tofile(char fName[])
Write a representation of the ph surface to file with O in each "good" spot and "X" in each "bad" one...
Definition: TTSE.cpp:1398
join
std::string join(std::vector< std::string > strings, char delim)
Definition: TTSE.cpp:458
CoolPropStateClassSI::drhodT_constp
double drhodT_constp(void)
Definition: CPState.cpp:2108
iD
Definition: GlobalConstants.h:13
CoolPropStateClassSI::rhoV
double rhoV(void)
Definition: CPState.h:266
TTSESinglePhaseTableClass::s_Trho
std::vector< std::vector< double > > s_Trho
Definition: TTSE.h:124
Fluid::params
struct OtherParameters params
Definition: FluidClass.h:220
TTSESinglePhaseTableClass::d2sdrho2_Trho
std::vector< std::vector< double > > d2sdrho2_Trho
Definition: TTSE.h:124
TTSESinglePhaseTableClass::alpha_bicubic
std::vector< double > alpha_bicubic
Definition: TTSE.h:113
TTSESinglePhaseTableClass::d2Tdh2
std::vector< std::vector< double > > d2Tdh2
Definition: TTSE.h:118
TTSETwoPhaseTableClass::logrho
std::vector< double > logrho
Definition: TTSE.h:49
CoolPropStateClassSI::dhdp_constT
double dhdp_constT(void)
Definition: CPState.cpp:2267
CoolPropStateClassSI::hL
double hL(void)
Definition: CPState.cpp:1330
CoolPropStateClassSI::d2hdT2_constrho
double d2hdT2_constrho(void)
Definition: CPState.cpp:2220
TTSESinglePhaseTableClass::d2kdrho2_Trho
std::vector< std::vector< double > > d2kdrho2_Trho
Definition: TTSE.h:128
pathExists
bool pathExists(std::string path)
Definition: TTSE.cpp:519
TTSESinglePhaseTableClass::DV
std::vector< double > DV
Definition: TTSE.h:106
TTSESinglePhaseTableClass::within_range_one_other_input
bool within_range_one_other_input(long iInput1, double Input1, long iOther, double Other, int buf)
Definition: TTSE.cpp:2157
TTSESinglePhaseTableClass::update_saturation_boundary_indices
void update_saturation_boundary_indices()
Definition: TTSE.cpp:1342
iPHASE_TWOPHASE
Definition: GlobalConstants.h:19
CoolPropStateClassSI::dsdp_along_sat_vapor
double dsdp_along_sat_vapor(void)
Definition: CPState.cpp:2445
TTSETwoPhaseTableClass::set_size
void set_size(unsigned int N)
Definition: TTSE.cpp:2564
TTSESinglePhaseTableClass::d2pdrho2_Trho
std::vector< std::vector< double > > d2pdrho2_Trho
Definition: TTSE.h:125
TTSESinglePhaseTableClass::rho_Trho
std::vector< double > rho_Trho
Definition: TTSE.h:129
iL
Definition: GlobalConstants.h:13
CoolPropStateClassSI::d2Tdp2_along_sat
double d2Tdp2_along_sat(void)
Second derivative of temperature w.r.t. pressure along saturation curve.
Definition: CPState.cpp:2391
iPHASE_LIQUID
Definition: GlobalConstants.h:19
TTSETwoPhaseTableClass::evaluate
double evaluate(long iParam, double p)
Definition: TTSE.cpp:2699
TTSESinglePhaseTableClass::Nrho
unsigned int Nrho
Definition: TTSE.h:94
TTSESinglePhaseTableClass::drhodp
std::vector< std::vector< double > > drhodp
Definition: TTSE.h:119
TTSESinglePhaseTableClass::d2mudTdrho_Trho
std::vector< std::vector< double > > d2mudTdrho_Trho
Definition: TTSE.h:127
TTSETwoPhaseTableClass::evaluate_T
double evaluate_T(double T)
Definition: TTSE.cpp:2786
TTSESinglePhaseTableClass::vector_from_file
void vector_from_file(std::string fName, int N, std::vector< double > *A)
Definition: TTSE.cpp:509
TTSETwoPhaseTableClass::d2rhodp2
std::vector< double > d2rhodp2
Definition: TTSE.h:49
CoolPropStateClassSI::p
double p(void)
Definition: CPState.h:291
CoolPropTools.h
CoolPropStateClassSI::d2pdrho2_constT
double d2pdrho2_constT(void)
Definition: CPState.cpp:2122
CoolPropStateClassSI::sV
double sV(void)
Definition: CPState.cpp:1363
CoolPropStateClassSI::d2sdT2_constrho
double d2sdT2_constrho(void)
Definition: CPState.cpp:2234
iV
Definition: GlobalConstants.h:13
CoolPropStateClassSI::h
double h(void)
Definition: CPState.cpp:1382
TTSETwoPhaseTableClass::h
std::vector< double > h
Definition: TTSE.h:51
TTSESinglePhaseTableClass::SatV
TTSETwoPhaseTableClass * SatV
Definition: TTSE.h:103
TTSESinglePhaseTableClass::write_all_to_file
void write_all_to_file(std::string root_path=std::string())
Definition: TTSE.cpp:745
TTSETwoPhaseTableClass::p
std::vector< double > p
Definition: TTSE.h:52
CoolPropStateClassSI::d2sdp2_along_sat_vapor
double d2sdp2_along_sat_vapor(void)
Definition: CPState.cpp:2451
OtherParameters::ptriple
double ptriple
Definition: FluidClass.h:43
TTSESinglePhaseTableClass::d2rhodhdp
std::vector< std::vector< double > > d2rhodhdp
Definition: TTSE.h:119
TTSESinglePhaseTableClass::SL
std::vector< double > SL
Definition: TTSE.h:106
TTSESinglePhaseTableClass::evaluate_randomly
double evaluate_randomly(long iParam, unsigned int N)
Definition: TTSE.cpp:1468
CoolPropStateClassSI::T
double T(void)
Definition: CPState.h:289
TTSESinglePhaseTableClass::d2mudT2_Trho
std::vector< std::vector< double > > d2mudT2_Trho
Definition: TTSE.h:127
CoolPropStateClassSI::s
double s(void)
Definition: CPState.cpp:1420
CoolPropStateClassSI::d2sdp2_along_sat_liquid
double d2sdp2_along_sat_liquid(void)
Definition: CPState.cpp:2459
TTSESinglePhaseTableClass::d2hdrho2_Trho
std::vector< std::vector< double > > d2hdrho2_Trho
Definition: TTSE.h:126
TTSETwoPhaseTableClass::evaluate_randomly
double evaluate_randomly(long iParam, unsigned int N)
Definition: TTSE.cpp:2890
CoolPropStateClassSI::rho
double rho(void)
Definition: CPState.h:290
TTSESinglePhaseTableClass::d2mudrho2_Trho
std::vector< std::vector< double > > d2mudrho2_Trho
Definition: TTSE.h:127
TTSESinglePhaseTableClass::build_Trho
double build_Trho(double Tmin, double Tmax, double rhomin, double rhomax, TTSETwoPhaseTableClass *SatL, TTSETwoPhaseTableClass *SatV)
Definition: TTSE.cpp:1051
TTSESinglePhaseTableClass::d2Tdp2
std::vector< std::vector< double > > d2Tdp2
Definition: TTSE.h:118
TTSESinglePhaseTableClass::read_all_from_file
bool read_all_from_file(std::string root_path)
Definition: TTSE.cpp:588
CoolPropStateClassSI::dsdT_constrho
double dsdT_constrho(void)
Definition: CPState.cpp:2230
CoolPropStateClassSI::dhdT_constrho
double dhdT_constrho(void)
Definition: CPState.cpp:2207
iPHASE_GAS
Definition: GlobalConstants.h:19
solve_cubic
void solve_cubic(double a, double b, double c, double d, double *x0, double *x1, double *x2)
Definition: CoolPropTools.cpp:188
TTSESinglePhaseTableClass::bicubic_evaluate_first_derivative_ph
double bicubic_evaluate_first_derivative_ph(long iOF, long iWRT, long iCONSTANT, double p, double logp, double h)
Definition: TTSE.cpp:1853
iT
Definition: GlobalConstants.h:13
TTSESinglePhaseTableClass::set_size_ph
void set_size_ph(unsigned int Np=100, unsigned int Nh=100)
Set the sizes of the matrices with h,p as inputs.
Definition: TTSE.cpp:169
TTSETwoPhaseTableClass::d2sdp2
std::vector< double > d2sdp2
Definition: TTSE.h:50
TTSESinglePhaseTableClass::dsdT_Trho
std::vector< std::vector< double > > dsdT_Trho
Definition: TTSE.h:124
CoolPropStateClassSI::dsdp_along_sat_liquid
double dsdp_along_sat_liquid(void)
Definition: CPState.cpp:2439
TTSETwoPhaseTableClass::T
std::vector< double > T
Definition: TTSE.h:48
TTSESinglePhaseTableClass::z_bicubic
std::vector< double > z_bicubic
Definition: TTSE.h:113
TTSETwoPhaseTableClass::drhodp
std::vector< double > drhodp
Definition: TTSE.h:49
Fluid::TTSESatL
TTSETwoPhaseTableClass TTSESatL
Definition: FluidClass.h:226
TTSESinglePhaseTableClass::d2sdhdp
std::vector< std::vector< double > > d2sdhdp
Definition: TTSE.h:120
TTSETwoPhaseTableClass::pFluid
Fluid * pFluid
Definition: TTSE.h:28
TTSESinglePhaseTableClass::bicubic_cell_coeffs_ph
std::vector< double > * bicubic_cell_coeffs_ph(long iParam, int i, int j)
Definition: TTSE.cpp:1570
TTSEREV
#define TTSEREV
Definition: TTSE.cpp:42
CoolPropStateClassSI::dpdrho_consth
double dpdrho_consth(void)
Definition: CPState.cpp:2262
TTSETwoPhaseTableClass::logp
std::vector< double > logp
Definition: TTSE.h:52
TTSE_MODE_BICUBIC
Definition: TTSE.h:8
iP
Definition: GlobalConstants.h:13
get_file_contents
std::string get_file_contents(const char *filename)
Definition: CoolPropTools.cpp:172
TTSESinglePhaseTableClass::TL
std::vector< double > TL
Definition: TTSE.h:106
CoolPropStateClassSI::d2pdrhodT
double d2pdrhodT(void)
Definition: CPState.cpp:2126
TTSETwoPhaseTableClass::evaluate_sat_derivative
double evaluate_sat_derivative(long iParam, double p)
Definition: TTSE.cpp:2847
TTSESinglePhaseTableClass::enable_writing_tables_to_files
bool enable_writing_tables_to_files
Definition: TTSE.h:109
TTSESinglePhaseTableClass::matrix_to_file
void matrix_to_file(std::string fName, std::vector< std::vector< double > > *A)
Definition: TTSE.cpp:468
TTSESinglePhaseTableClass::d2kdTdrho_Trho
std::vector< std::vector< double > > d2kdTdrho_Trho
Definition: TTSE.h:128
TTSESinglePhaseTableClass::d2pdTdrho_Trho
std::vector< std::vector< double > > d2pdTdrho_Trho
Definition: TTSE.h:125
TTSESinglePhaseTableClass::d2pdT2_Trho
std::vector< std::vector< double > > d2pdT2_Trho
Definition: TTSE.h:125
FluidLimits::Tmin
double Tmin
Definition: FluidClass.h:54
TTSESinglePhaseTableClass::d2hdT2_Trho
std::vector< std::vector< double > > d2hdT2_Trho
Definition: TTSE.h:126
CoolPropStateClassSI::d2pdT2_constrho
double d2pdT2_constrho(void)
Definition: CPState.cpp:2130
TTSESinglePhaseTableClass::vector_to_file
void vector_to_file(std::string fName, std::vector< double > *A)
Definition: TTSE.cpp:499
TTSESinglePhaseTableClass::Nh
unsigned int Nh
Definition: TTSE.h:94
TTSESinglePhaseTableClass::dTdh
std::vector< std::vector< double > > dTdh
Definition: TTSE.h:118
CoolPropStateClassSI::d2hdrho2_constT
double d2hdrho2_constT(void)
Definition: CPState.cpp:2211
TTSESinglePhaseTableClass::nearest_good_neighbor
void nearest_good_neighbor(int *i, int *j)
Definition: TTSE.cpp:257
TTSESinglePhaseTableClass::dsdh
std::vector< std::vector< double > > dsdh
Definition: TTSE.h:120
CoolPropStateClassSI::dTdp_along_sat
double dTdp_along_sat(void)
Derivative of temperature w.r.t. pressure along saturation curve.
Definition: CPState.cpp:2366
IPropsSI
EXPORT_CODE double CONVENTION IPropsSI(long iOutput, long iName1, double Prop1, long iName2, double Prop2, long iFluid)
Definition: CoolProp.cpp:861
TTSESinglePhaseTableClass::d2hdTdrho_Trho
std::vector< std::vector< double > > d2hdTdrho_Trho
Definition: TTSE.h:126
TTSESinglePhaseTableClass::dpdrho_Trho
std::vector< std::vector< double > > dpdrho_Trho
Definition: TTSE.h:125
ValidNumber
bool ValidNumber(double x)
Definition: CoolPropTools.h:117
TTSESinglePhaseTableClass::d2kdT2_Trho
std::vector< std::vector< double > > d2kdT2_Trho
Definition: TTSE.h:128
TTSESinglePhaseTableClass::evaluate_one_other_input
double evaluate_one_other_input(long iInput1, double Param1, long iOther, double Other)
Evaluate the TTSE using P,S or P,T.
Definition: TTSE.cpp:2188
TTSESinglePhaseTableClass::evaluate_Trho
double evaluate_Trho(long iOutput, double T, double rho, double logrho)
Definition: TTSE.cpp:2448
CoolPropStateClassSI::dpdT_constrho
double dpdT_constrho(void)
Definition: CPState.cpp:2118
matrix_vector_product
void matrix_vector_product(std::vector< double > *x, std::vector< double > *y)
Definition: TTSE.cpp:108
TTSESinglePhaseTableClass::T_Trho
std::vector< double > T_Trho
Definition: TTSE.h:129
TTSESinglePhaseTableClass::d2sdp2
std::vector< std::vector< double > > d2sdp2
Definition: TTSE.h:120
TTSESinglePhaseTableClass::d2rhodh2
std::vector< std::vector< double > > d2rhodh2
Definition: TTSE.h:119
TTSESinglePhaseTableClass::matrix_from_file
void matrix_from_file(std::string fName, std::vector< std::vector< double > > *A)
Definition: TTSE.cpp:482
TTSESinglePhaseTableClass::NT
unsigned int NT
Definition: TTSE.h:94
CoolPropStateClassSI::drhodp_along_sat_liquid
double drhodp_along_sat_liquid(void)
Definition: CPState.cpp:2482
Generated on Sun Aug 10 2014 16:02:42 for CoolProp by   doxygen 1.8.7