LORENE
eos_fitting.C
1 /*
2  * Method of class Eos_fitting
3  *
4  * (see file eos_fitting.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Keisuke Taniguchi
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char eos_fitting_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_fitting.C,v 1.5 2014/10/13 08:52:53 j_novak Exp $" ;
29 
30 /*
31  * $Id: eos_fitting.C,v 1.5 2014/10/13 08:52:53 j_novak Exp $
32  * $Log: eos_fitting.C,v $
33  * Revision 1.5 2014/10/13 08:52:53 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.4 2014/10/06 15:13:06 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.3 2005/05/23 14:14:00 k_taniguchi
40  * Minor modification.
41  *
42  * Revision 1.2 2005/05/22 20:53:06 k_taniguchi
43  * Modify the method to calculate baryon number density from enthalpy.
44  *
45  * Revision 1.1 2004/09/26 18:53:53 k_taniguchi
46  * Initial revision
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_fitting.C,v 1.5 2014/10/13 08:52:53 j_novak Exp $
50  *
51  */
52 
53 // C headers
54 #include <cstdlib>
55 #include <cstring>
56 #include <cmath>
57 
58 // Lorene headers
59 #include "headcpp.h"
60 #include "eos_fitting.h"
61 #include "eos.h"
62 #include "utilitaires.h"
63 #include "unites.h"
64 
65 namespace Lorene {
66 double fc(double) ;
67 double gc(double) ;
68 double hc(double) ;
69 
70 //************************************************************************
71 
72  //--------------------------------//
73  // Constructors //
74  //--------------------------------//
75 
76 // Standard constructor
77 // --------------------
78 Eos_fitting::Eos_fitting(const char* name_i, const char* data,
79  const char* path) : Eos(name_i) {
80 
81  strcpy(dataname, path) ;
82  strcat(dataname, "/") ;
83  strcat(dataname, data) ;
84 
85  read_coef() ;
86 
87 }
88 
89 // Constructor from a binary file
90 // ------------------------------
91 Eos_fitting::Eos_fitting(FILE* fich) : Eos(fich) {
92 
93  fread(dataname, sizeof(char), 160, fich) ;
94 
95  read_coef() ;
96 
97 }
98 
99 // Constructor from a formatted file
100 // ---------------------------------
101 Eos_fitting::Eos_fitting(ifstream& fich, const char* data) : Eos(fich) {
102 
103  char path[160] ;
104 
105  fich.getline(path, 160) ;
106 
107  strcpy(dataname, path) ;
108  strcat(dataname, "/") ;
109  strcat(dataname, data) ;
110 
111  read_coef() ;
112 
113 }
114 
115 // Destructor
117 
118  delete [] pp ;
119 
120 }
121 
122 
123  //---------------------------------------//
124  // Outputs //
125  //---------------------------------------//
126 
127 void Eos_fitting::sauve(FILE* fich) const {
128 
129  Eos::sauve(fich) ;
130 
131  fwrite(dataname, sizeof(char), 160, fich) ;
132 
133 }
134  //-----------------------//
135  // Miscellaneous //
136  //-----------------------//
137 
139 
140  char blabla[120] ;
141 
142  ifstream fich(dataname) ;
143 
144  for (int i=0; i<3; i++) { // Jump over the file header
145  fich.getline(blabla, 120) ;
146  }
147 
148  int nb_coef ;
149  fich >> nb_coef ; fich.getline(blabla, 120) ; // Number of coefficients
150 
151  for (int i=0; i<3; i++) { // Jump over the table header
152  fich.getline(blabla, 120) ;
153  }
154 
155  pp = new double[nb_coef] ;
156 
157  for (int i=0; i<nb_coef; i++) {
158  fich >> pp[i] ; fich.getline(blabla, 120) ;
159  }
160 
161  fich.close() ;
162 
163 }
164 
165 
166  //------------------------------//
167  // Computational routines //
168  //------------------------------//
169 
170 // Baryon density from enthalpy
171 //------------------------------
172 
173 double Eos_fitting::nbar_ent_p(double ent, const Param* ) const {
174 
175  using namespace Unites ;
176 
177  if ( ent > double(0) ) {
178 
179  double aa = 0. ;
180  double xx = 0.01 ;
181  int m ;
182  double yy ;
183  double ent_value ;
184  double rhob ;
185  double nb ;
186  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
187  / rhonuc_si ;
188 
189  while (xx > 1.e-15) {
190 
191  ent_value = 1. ; // Initialization
192  xx = 0.1 * xx ;
193  m = 0 ;
194 
195  while (ent_value > 1.e-15) {
196 
197  m++ ;
198  yy = aa + m * xx ;
199 
200  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
201  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
202  double ccc = pp[0]*pp[1]*pow(yy,pp[1])
203  +pp[2]*pp[3]*pow(yy,pp[3]) ;
204  double ddd = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-1.) ;
205  double eee = -pp[7]*yy+pp[9] ;
206  double fff = -pp[8]*yy+pp[10] ;
207  double ggg = pp[4]*pp[5]*pp[6]*pow(yy,pp[5]) ;
208 
209  ent_value = exp(ent) - 1.0
210  -( aaa*bbb - 1. ) * fc(eee)
211  -pp[11]*pow(yy,pp[12])*fc(-eee)*fc(fff)
212  -pp[13]*pow(yy,pp[14])*fc(-fff)
213  -( ccc*bbb + aaa*ddd*ggg )*fc(eee)
214  -yy*( aaa*bbb - 1. )*pp[7]*gc(eee)
215  -pp[11]*pp[12]*pow(yy,pp[12])*fc(-eee)*fc(fff)
216  +pp[11]*pow(yy,pp[12]+1.)*(pp[7]*gc(-eee)*fc(fff)
217  -pp[8]*fc(-eee)*gc(fff))
218  -pp[13]*pow(yy,pp[14])*(pp[14]*fc(-fff)
219  -pp[8]*yy*gc(-fff)) ;
220 
221  }
222  aa += (m - 1) * xx ;
223  }
224  rhob = aa ;
225 
226  // The transformation from rhob to nb
227  nb = rhob * trans_dens ;
228 
229  return nb ;
230  }
231  else {
232  return 0 ;
233  }
234 
235 }
236 
237 // Energy density from enthalpy
238 //------------------------------
239 
240 double Eos_fitting::ener_ent_p(double ent, const Param* ) const {
241 
242  using namespace Unites ;
243 
244  if ( ent > double(0) ) {
245 
246  // Number density in the unit of [n_nuc]
247  double nb = nbar_ent_p(ent) ;
248 
249  // The transformation from nb to yy
250  // --------------------------------
251 
252  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
253  / rhonuc_si ;
254 
255  // Baryon density in the unit of c=G=Msol=1
256  double yy = nb / trans_dens ;
257 
258  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
259  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
260  double eee = -pp[7]*yy+pp[9] ;
261  double fff = -pp[8]*yy+pp[10] ;
262 
263  double epsil = ( aaa*bbb - 1. ) * fc(eee)
264  +pp[11]*pow(yy,pp[12])*fc(-eee)*fc(fff)
265  +pp[13]*pow(yy,pp[14])*fc(-fff) ;
266 
267  // The transformation from epsil to ee
268  // -----------------------------------
269 
270  // Energy density in the unit of [rho_nuc * c^2]
271  double ee = nb * (1. + epsil) ;
272 
273  return ee ;
274  }
275  else {
276  return 0 ;
277  }
278 
279 }
280 
281 // Pressure from enthalpy
282 //------------------------
283 
284 double Eos_fitting::press_ent_p(double ent, const Param* ) const {
285 
286  using namespace Unites ;
287 
288  if ( ent > double(0) ) {
289 
290  // Number density in the unit of [n_nuc]
291  double nb = nbar_ent_p(ent) ;
292 
293  // The transformation from nb to yy
294  // --------------------------------
295 
296  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
297  / rhonuc_si ;
298 
299  // Baryon density in the unit of c=G=Msol=1
300  double yy = nb / trans_dens ;
301 
302  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
303  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
304  double ccc = pp[0]*pp[1]*pow(yy,pp[1])+pp[2]*pp[3]*pow(yy,pp[3]) ;
305  double ddd = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-1.) ;
306  double eee = -pp[7]*yy+pp[9] ;
307  double fff = -pp[8]*yy+pp[10] ;
308  double ggg = pp[4]*pp[5]*pp[6]*pow(yy,pp[5]) ;
309 
310  double ppp = yy*( ccc*bbb + aaa*ddd*ggg )*fc(eee)
311  +yy*yy*( aaa*bbb - 1. )*pp[7]*gc(eee)
312  +pp[11]*pp[12]*pow(yy,pp[12]+1.)*fc(-eee)*fc(fff)
313  -pp[11]*pow(yy,pp[12]+2.)*(pp[7]*gc(-eee)*fc(fff)
314  -pp[8]*fc(-eee)*gc(fff))
315  +pp[13]*pow(yy,pp[14]+1.)*(pp[14]*fc(-fff)-pp[8]*yy*gc(-fff)) ;
316 
317  // The transformation from ppp to pres
318  // -----------------------------------
319 
320  // Pressure in the unit of [rho_nuc * c^2]
321  double pres = ppp * trans_dens ;
322 
323  return pres ;
324  }
325  else {
326  return 0 ;
327  }
328 
329 }
330 
331 // dln(n)/dln(H) from enthalpy
332 //----------------------------
333 
334 double Eos_fitting::der_nbar_ent_p(double ent, const Param* ) const {
335 
336  using namespace Unites ;
337 
338  if ( ent > double(0) ) {
339 
340  // Number density in the unit of [n_nuc]
341  double nb = nbar_ent_p(ent) ;
342 
343  // The transformation from nb to yy
344  // --------------------------------
345 
346  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
347  / rhonuc_si ;
348 
349  // Baryon density in the unit of c=G=Msol=1
350  double yy = nb / trans_dens ;
351 
352  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
353  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
354  double ccc = pp[0]*pp[1]*pow(yy,pp[1]) + pp[2]*pp[3]*pow(yy,pp[3]) ;
355  double ddd = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-1.) ;
356  double eee = -pp[7]*yy+pp[9] ;
357  double fff = -pp[8]*yy+pp[10] ;
358  double ggg = pp[4]*pp[5]*pp[6]*pow(yy,pp[5]) ;
359  double jjj = pp[0]*pp[1]*pp[1]*pow(yy,pp[1])
360  +pp[2]*pp[3]*pp[3]*pow(yy,pp[3]) ;
361 
362  double dlnsdlh = exp(ent) * ent /
363  ( ( ccc*bbb + aaa*ddd*ggg )*( fc(eee) + 2.*yy*pp[7]*gc(eee) )
364  +yy*pp[7]*( aaa*bbb - 1. )*( 2.*gc(eee) - yy*pp[7]*hc(eee) )
365  +pp[11]*pp[12]*(pp[12]+1.)*pow(yy,pp[12])*fc(-eee)*fc(fff)
366  +2.*pp[11]*(pp[12]+1.)*pow(yy,pp[12]+1.)
367  *( -pp[7]*gc(-eee)*fc(fff) + pp[8]*fc(-eee)*gc(fff) )
368  -pp[11]*pow(yy,pp[12]+2.)*( pp[7]*pp[7]*hc(-eee)*fc(fff)
369  +2.*pp[7]*pp[8]*gc(-eee)*gc(fff)
370  +pp[8]*pp[8]*fc(-eee)*hc(fff) )
371  +pp[13]*pp[14]*(pp[14]+1.)*pow(yy,pp[14])*fc(-fff)
372  -2.*pp[8]*pp[13]*(pp[14]+1.)*pow(yy,pp[14]+1.)*gc(-fff)
373  -pp[8]*pp[8]*pp[13]*pow(yy,pp[14]+2.)*hc(-fff)
374  +( jjj*bbb + 2.*ccc*ddd*ggg
375  + aaa*pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-2.)
376  *ggg*(ggg+pp[5]) )*fc(eee)
377  ) ;
378 
379  return dlnsdlh ;
380 
381  }
382  else {
383 
384  return double(1) / pp[12] ; // To ensure continuity at ent=0
385 
386  }
387 
388 }
389 
390 // dln(e)/dln(H) from enthalpy
391 //----------------------------
392 
393 double Eos_fitting::der_ener_ent_p(double ent, const Param* ) const {
394 
395  using namespace Unites ;
396 
397  if ( ent > double(0) ) {
398 
399  // Number density in the unit of [n_nuc]
400  double nb = nbar_ent_p(ent) ;
401 
402  // The transformation from nb to yy
403  // --------------------------------
404 
405  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
406  / rhonuc_si ;
407 
408  // Baryon density in the unit of c=G=Msol=1
409  double yy = nb / trans_dens ;
410 
411  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
412  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
413  double ccc = pp[0]*pp[1]*pow(yy,pp[1]) + pp[2]*pp[3]*pow(yy,pp[3]) ;
414  double ddd = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-1.) ;
415  double eee = -pp[7]*yy+pp[9] ;
416  double fff = -pp[8]*yy+pp[10] ;
417  double ggg = pp[4]*pp[5]*pp[6]*pow(yy,pp[5]) ;
418 
419  double dlnsdlh = der_nbar_ent_p(ent) ;
420 
421  double dlesdlh = dlnsdlh
422  * (1. + ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
423  +yy*pp[7]*( aaa*bbb - 1. )*gc(eee)
424  +pp[11]*pp[12]*pow(yy,pp[12])*fc(-eee)*fc(fff)
425  +pp[11]*pow(yy,pp[12]+1.)*( -pp[7]*gc(-eee)*fc(fff)
426  +pp[8]*fc(-eee)*gc(fff) )
427  +pp[13]*pp[14]*pow(yy,pp[14])*fc(-fff)
428  -pp[8]*pp[13]*pow(yy,pp[14]+1.)*gc(-fff) )
429  / ( 1. + ( aaa*bbb - 1. )*fc(eee)
430  + pp[11]*pow(yy,pp[12])*fc(-eee)*fc(fff)
431  + pp[13]*pow(yy,pp[14])*fc(-fff) ) ) ;
432 
433  return dlesdlh ;
434 
435  }
436  else {
437 
438  return double(1) / pp[12] ; // To ensure continuity at ent=0
439 
440  }
441 
442 }
443 
444 // dln(p)/dln(H) from enthalpy
445 //----------------------------
446 
447 double Eos_fitting::der_press_ent_p(double ent, const Param* ) const {
448 
449  using namespace Unites ;
450 
451  if ( ent > double(0) ) {
452 
453  // Number density in the unit of [n_nuc]
454  double nb = nbar_ent_p(ent) ;
455 
456  // The transformation from nb to yy
457  // --------------------------------
458 
459  double trans_dens = msol_si / pow(g_si*msol_si/c_si/c_si,3.)
460  / rhonuc_si ;
461 
462  // Baryon density in the unit of c=G=Msol=1
463  double yy = nb / trans_dens ;
464 
465  double aaa = 1.+pp[0]*pow(yy,pp[1])+pp[2]*pow(yy,pp[3]) ;
466  double bbb = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]) ;
467  double ccc = pp[0]*pp[1]*pow(yy,pp[1]) + pp[2]*pp[3]*pow(yy,pp[3]) ;
468  double ddd = pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-1.) ;
469  double eee = -pp[7]*yy+pp[9] ;
470  double fff = -pp[8]*yy+pp[10] ;
471  double ggg = pp[4]*pp[5]*pp[6]*pow(yy,pp[5]) ;
472  double jjj = pp[0]*pp[1]*pp[1]*pow(yy,pp[1])
473  +pp[2]*pp[3]*pp[3]*pow(yy,pp[3]) ;
474 
475  double dlnsdlh = der_nbar_ent_p(ent) ;
476 
477  double dlpsdlh = dlnsdlh
478  * ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
479  +( jjj*bbb + 2.*ccc*ddd*ggg
480  + aaa*pow(1.+pp[4]*pow(yy,pp[5]),pp[6]-2.)*ggg*(ggg+pp[5])
481  )*fc(eee)
482  +2.*yy*pp[7]*( ccc*bbb + aaa*ddd*ggg )*gc(eee)
483  +yy*pp[7]*( aaa*bbb - 1. )*(2.*gc(eee) - yy*pp[7]*hc(eee))
484  +pp[11]*pp[12]*(pp[12]+1.)*pow(yy,pp[12])*fc(-eee)*fc(fff)
485  +2.*pp[11]*(pp[12]+1.)*pow(yy,pp[12]+1.)
486  *( -pp[7]*gc(-eee)*fc(fff) + pp[8]*fc(-eee)*gc(fff) )
487  -pp[11]*pow(yy,pp[12]+2.)*( pp[7]*pp[7]*hc(-eee)*fc(fff)
488  +2.*pp[7]*pp[8]*gc(-eee)*gc(fff)
489  +pp[8]*pp[8]*fc(-eee)*hc(fff) )
490  +pp[13]*(pp[14]+1.)*pow(yy,pp[14])*( pp[14]*fc(-fff)
491  -2.*pp[8]*yy*gc(-fff) )
492  -pp[8]*pp[8]*pp[13]*pow(yy,pp[14]+2.)*hc(-fff) )
493  / ( ( ccc*bbb + aaa*ddd*ggg )*fc(eee)
494  +yy*pp[7]*( aaa*bbb - 1. )*gc(eee)
495  +pp[11]*pow(yy,pp[12])*( pp[12]*fc(-eee)*fc(fff)
496  -yy*pp[7]*gc(-eee)*fc(fff)
497  +yy*pp[8]*fc(-eee)*gc(fff) )
498  +pp[13]*pow(yy,pp[14])*( pp[14]*fc(-fff)
499  -yy*pp[8]*gc(-fff) ) ) ;
500 
501  return dlpsdlh ;
502 
503  }
504  else {
505 
506  return (pp[12] + 1.) / pp[12] ; // To ensure continuity at ent=0
507 
508  }
509 
510 }
511 
512 //********************************************************
513 // Functions which appear in the fitting formula
514 //********************************************************
515 
516 double fc(double xx) {
517 
518  double resu = double(1) / (double(1) + exp(xx)) ;
519 
520  return resu ;
521 
522 }
523 
524 double gc(double xx) {
525 
526  double resu ;
527 
528  if (xx > 0.) {
529  resu = exp(-xx) / pow(exp(-xx)+double(1),double(2)) ;
530  }
531  else {
532  resu = exp(xx) / pow(double(1)+exp(xx),double(2)) ;
533  }
534 
535  return resu ;
536 
537 }
538 
539  double hc(double xx) {
540 
541  double resu ;
542 
543  if (xx > 0.) {
544  resu = exp(-xx) * (exp(-xx)-double(1)) /
545  pow(exp(-xx)+double(1),double(3)) ;
546  }
547  else {
548  resu = exp(xx) * (double(1)-exp(xx)) /
549  pow(double(1)+exp(xx),double(3)) ;
550  }
551 
552  return resu ;
553 
554  }
555 }
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_fitting.C:393
virtual ~Eos_fitting()
Destructor.
Definition: eos_fitting.C:116
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_fitting.C:334
void read_coef()
Reading coefficients of the fitting equation for the energy density, the pressure,...
Definition: eos_fitting.C:138
char dataname[160]
Name of the file containing the fitting data.
Definition: eos_fitting.h:86
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Definition: eos_fitting.C:447
double * pp
Array of the coefficients of the fitting data.
Definition: eos_fitting.h:89
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition: eos_fitting.C:173
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition: eos_fitting.C:240
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_fitting.C:284
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_fitting.C:127
Eos_fitting(const char *name_i, const char *data, const char *path)
Standard constructor.
Definition: eos_fitting.C:78
Equation of state base class.
Definition: eos.h:190
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:179
Parameter storage.
Definition: param.h:125
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.