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
68namespace Lorene {
69double fc(double) ;
70double gc(double) ;
71double hc(double) ;
72
73//************************************************************************
74
75 //--------------------------------//
76 // Constructors //
77 //--------------------------------//
78
79// Standard constructor
80// --------------------
81Eos_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// ------------------------------
94Eos_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// ---------------------------------
104Eos_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
119Eos_fitting::~Eos_fitting() {
120
121 delete [] pp ;
122
123}
124
125
126 //---------------------------------------//
127 // Outputs //
128 //---------------------------------------//
129
130void 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
141void 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
176double 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
243double 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
287double 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
337double 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
396double 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
450double 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
519double fc(double xx) {
520
521 double resu = double(1) / (double(1) + exp(xx)) ;
522
523 return resu ;
524
525}
526
527double 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::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::~Eos_fitting
virtual ~Eos_fitting()
Destructor.
Definition eos_fitting.C:119
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::Eos_fitting::read_coef
void read_coef()
Reading coefficients of the fitting equation for the energy density, the pressure,...
Definition eos_fitting.C:141
Lorene::Eos_fitting::dataname
char dataname[160]
Name of the file containing the fitting data.
Definition eos_fitting.h:89
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
Lorene::Eos_fitting::pp
double * pp
Array of the coefficients of the fitting data.
Definition eos_fitting.h:92
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::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::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::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos_fitting.C:130
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::Eos::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:189
Lorene::Eos::Eos
Eos()
Standard constructor.
Definition eos.C:118
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene
Lorene prototypes.
Definition app_hor.h:67
Unites
Standard units of space, time and mass.