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
28char 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
65namespace Lorene {
66double fc(double) ;
67double gc(double) ;
68double hc(double) ;
69
70//************************************************************************
71
72 //--------------------------------//
73 // Constructors //
74 //--------------------------------//
75
76// Standard constructor
77// --------------------
78Eos_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// ------------------------------
91Eos_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// ---------------------------------
101Eos_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
127void 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
173double 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
240double 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
284double 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
334double 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
393double 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
447double 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
516double fc(double xx) {
517
518 double resu = double(1) / (double(1) + exp(xx)) ;
519
520 return resu ;
521
522}
523
524double 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.
virtual ~Eos_fitting()
Destructor.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void read_coef()
Reading coefficients of the fitting equation for the energy density, the pressure,...
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.
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.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
virtual void sauve(FILE *) const
Save in a file.
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.