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