LORENE
eos_strange_cr.C
1/*
2 * Methods for the class Eos_strange_cr_cr
3 *
4 * (see file eos.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2000 J. Leszek Zdunik
10 * Copyright (c) 2000-2001 Eric Gourgoulhon
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30
31char Eos_strange_cr_cr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_strange_cr.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $" ;
32
33/*
34 * $Id: eos_strange_cr.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $
35 * $Log: eos_strange_cr.C,v $
36 * Revision 1.7 2014/10/13 08:52:54 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.6 2014/10/06 15:13:07 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.5 2004/03/25 10:29:02 j_novak
43 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
44 *
45 * Revision 1.4 2002/10/16 14:36:35 j_novak
46 * Reorganization of #include instructions of standard C++, in order to
47 * use experimental version 3 of gcc.
48 *
49 * Revision 1.3 2002/04/09 14:32:15 e_gourgoulhon
50 * 1/ Added extra parameters in EOS computational functions (argument par)
51 * 2/ New class MEos for multi-domain EOS
52 *
53 * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
54 *
55 * All writing/reading to a binary file are now performed according to
56 * the big endian convention, whatever the system is big endian or
57 * small endian, thanks to the functions fwrite_be and fread_be
58 *
59 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
60 * LORENE
61 *
62 * Revision 2.2 2001/02/07 09:49:26 eric
63 * Suppression de la fonction derent_ent_p.
64 * Ajout des fonctions donnant les derivees de l'EOS:
65 * der_nbar_ent_p
66 * der_ener_ent_p
67 * der_press_ent_p
68 * /
69 *
70 * Revision 2.1 2000/11/24 13:26:50 eric
71 * Correction erreur dans set_auxillary(): rho_nd_nucl est un membre !
72 *
73 * Revision 2.0 2000/11/23 14:46:35 eric
74 * *** empty log message ***
75 *
76 *
77 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_strange_cr.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $
78 *
79 */
80
81
82// Headers C
83#include <cstdlib>
84#include <cstring>
85#include <cmath>
86
87// Headers Lorene
88#include "eos.h"
89#include "cmp.h"
90#include "utilitaires.h"
91#include "unites.h"
92
93 //------------------------------------//
94 // Constructors //
95 //------------------------------------//
96
97// Standard constructor
98// --------------------
99namespace Lorene {
100Eos_strange_cr::Eos_strange_cr(double n0_b60_i, double b60_i, double ent0_i,
101 double eps_fit_i, double rho0_b60_i,
102 double ent_nd_i, double rho_nd_i, double gam_i) :
103 Eos("Strange matter EOS from Zdunik (2000) with crust"),
104 n0_b60(n0_b60_i),
105 b60(b60_i),
106 ent0(ent0_i),
107 eps_fit(eps_fit_i),
108 rho0_b60(rho0_b60_i),
109 ent_nd(ent_nd_i),
110 rho_nd(rho_nd_i),
111 gam(gam_i) {
112
113 set_auxiliary() ;
114
115}
116
117// Copy constructor
118// -----------------
119
121 Eos(eos_i),
122 n0_b60(eos_i.n0_b60),
123 b60(eos_i.b60),
124 ent0(eos_i.ent0),
125 eps_fit(eos_i.eps_fit),
126 rho0_b60(eos_i.rho0_b60),
127 ent_nd(eos_i.ent_nd),
128 rho_nd(eos_i.rho_nd),
129 gam(eos_i.gam) {
130
131 set_auxiliary() ;
132
133}
134
135
136// Constructor from binary file
137// ----------------------------
139 Eos(fich) {
140
141 fread_be(&n0_b60, sizeof(double), 1, fich) ;
142 fread_be(&b60, sizeof(double), 1, fich) ;
143 fread_be(&ent0, sizeof(double), 1, fich) ;
144 fread_be(&eps_fit, sizeof(double), 1, fich) ;
145 fread_be(&rho0_b60, sizeof(double), 1, fich) ;
146 fread_be(&ent_nd, sizeof(double), 1, fich) ;
147 fread_be(&rho_nd, sizeof(double), 1, fich) ;
148 fread_be(&gam, sizeof(double), 1, fich) ;
149
150 set_auxiliary() ;
151
152}
153
154// Constructor from a formatted file
155// ---------------------------------
157 Eos(fich) {
158
159 char blabla[80] ;
160
161 fich >> n0_b60 ; fich.getline(blabla, 80) ;
162 fich >> b60 ; fich.getline(blabla, 80) ;
163 fich >> ent0 ; fich.getline(blabla, 80) ;
164 fich >> eps_fit ; fich.getline(blabla, 80) ;
165 fich >> rho0_b60 ; fich.getline(blabla, 80) ;
166 fich >> ent_nd ; fich.getline(blabla, 80) ;
167 fich >> rho_nd ; fich.getline(blabla, 80) ;
168 fich >> gam ; fich.getline(blabla, 80) ;
169
170 set_auxiliary() ;
171
172}
173 //--------------//
174 // Destructor //
175 //--------------//
176
178
179 // does nothing
180
181}
182
183 //--------------//
184 // Assignment //
185 //--------------//
186
188
189 set_name(eosi.name) ;
190
191 n0_b60 = eosi.n0_b60 ;
192 b60 = eosi.b60 ;
193 ent0 = eosi.ent0 ;
194 eps_fit = eosi.eps_fit ;
195 rho0_b60 = eosi.rho0_b60 ;
196 ent_nd = eosi.ent_nd ;
197 rho_nd = eosi.rho_nd ;
198 gam = eosi.gam ;
199
200 set_auxiliary() ;
201
202}
203
204
205 //-----------------------//
206 // Miscellaneous //
207 //-----------------------//
208
210
211 using namespace Unites ;
212
213 rho_nd_nucl = rho_nd * mevpfm3 ;
214
215 rho0 = b60 * rho0_b60 * mevpfm3 ;
216
217 b34 = pow(b60, double(0.75)) ;
218
219 n0 = b34 * n0_b60 * double(10) ; // 10 : fm^{-3} --> 0.1 fm^{-3}
220
221 fach = (double(4) + eps_fit) / (double(1) + eps_fit) ;
222
223 x_nd = (exp (ent_nd) - 1. )/( 1. +exp(ent_nd)/( gam - 1. ) );
224
226
227 ncr_nd = (1. + x_nd) * n0 * rho_nd_nucl/rho0
228 / exp(ent_nd) / exp(delent);
229
230 unsgam1 = double(1) / ( gam - double(1) ) ;
231
232 gam1sx = ( gam - double(1) - x_nd) / gam / x_nd ;
233
234
235 cout << "ncr_nd =" << ncr_nd << endl ;
236
237 cout << "eos ent_nd, x_nd, delent, ent_nd +delent " << ent_nd << " "<< x_nd << " "
238 << delent <<" "<< (ent_nd+delent) << endl ;
239
240 double pcr = x_nd * rho_nd_nucl *
241 pow(( gam - 1. - x_nd)/gam/x_nd *
242 fabs((exp(ent_nd)-1.)), gam/(gam-1.)) ;
243
244 double psq = ( exp(fach * (ent_nd + delent)) - 1) / fach
245 * rho0 ;
246
247 cout << " eos sqm fach =" << fach << " PNS_s " <<
248 psq << " PNS_cr = " << pcr << endl ;
249 //double relp=psq/pcr ;
250 cout << " 1+fach ... " << (1+fach*x_nd*rho_nd_nucl/rho0) << endl;
251 cout << " log(miu) " <<
252 log(pow((1+fach*x_nd*rho_nd_nucl/rho0),1./fach)) << endl;
253 double miu_rate = rho0/n0*pow((1+fach*x_nd*rho_nd_nucl/rho0),1./fach)/
254 rho_nd_nucl/(1-x_nd/(gam-1.))*ncr_nd/exp(ent_nd) ;
255 cout << " miu_rate -1 " << miu_rate -1. ;
256
257 double aa0 = log(1+fach*x_nd*rho_nd_nucl/rho0) / fach ;
258 double aa1 = log(pow(1+fach*x_nd*rho_nd_nucl/rho0, 1./fach)) ;
259 cout << "aa0, aa1, aa0-aa1 : " << aa0 << " " << aa1 << " " <<
260 aa0-aa1 << endl ;
261 cout << "fach : " << fach << endl ;
262 cout << "1+fach*x_nd*rho_nd_nucl/rho0 : " << 1+fach*x_nd*rho_nd_nucl/rho0 << endl ;
263
264}
265
266
267 //------------------------//
268 // Comparison operators //
269 //------------------------//
270
271
272bool Eos_strange_cr::operator==(const Eos& eos_i) const {
273
274 bool resu = true ;
275
276 if ( eos_i.identify() != identify() ) {
277 cout << "The second EOS is not of type Eos_strange_cr !" << endl ;
278 resu = false ;
279 }
280 else{
281
282 const Eos_strange_cr& eos = dynamic_cast<const Eos_strange_cr&>( eos_i ) ;
283
284 if (eos.n0_b60 != n0_b60) {
285 cout
286 << "The two Eos_strange_cr have different n0_b60 : " << n0_b60 << " <-> "
287 << eos.n0_b60 << endl ;
288 resu = false ;
289 }
290
291 if (eos.b60 != b60) {
292 cout
293 << "The two Eos_strange_cr have different b60 : " << b60 << " <-> "
294 << eos.b60 << endl ;
295 resu = false ;
296 }
297
298 if (eos.ent0 != ent0) {
299 cout
300 << "The two Eos_strange_cr have different ent0 : " << ent0 << " <-> "
301 << eos.ent0 << endl ;
302 resu = false ;
303 }
304
305 if (eos.eps_fit != eps_fit) {
306 cout
307 << "The two Eos_strange_cr have different eps_fit : " << eps_fit
308 << " <-> " << eos.eps_fit << endl ;
309 resu = false ;
310 }
311
312 if (eos.rho0_b60 != rho0_b60) {
313 cout
314 << "The two Eos_strange_cr have different rho0_b60 : " << rho0_b60
315 << " <-> " << eos.rho0_b60 << endl ;
316 resu = false ;
317 }
318
319
320 if (eos.ent_nd != ent_nd) {
321 cout
322 << "The two Eos_strange_cr have different ent_nd : "
323 << ent_nd
324 << " <-> " << eos.ent_nd << endl ;
325 resu = false ;
326 }
327
328
329 if (eos.rho_nd != rho_nd) {
330 cout
331 << "The two Eos_strange_cr have different rho_nd : "
332 << rho_nd
333 << " <-> " << eos.rho_nd << endl ;
334 resu = false ;
335 }
336
337
338 if (eos.gam != gam) {
339 cout
340 << "The two Eos_strange_cr have different gam : " << gam
341 << " <-> " << eos.gam << endl ;
342 resu = false ;
343 }
344
345
346
347 }
348
349 return resu ;
350
351}
352
353bool Eos_strange_cr::operator!=(const Eos& eos_i) const {
354
355 return !(operator==(eos_i)) ;
356
357}
358
359 //------------//
360 // Outputs //
361 //------------//
362
363void Eos_strange_cr::sauve(FILE* fich) const {
364
365 Eos::sauve(fich) ;
366
367 fwrite_be(&n0_b60, sizeof(double), 1, fich) ;
368 fwrite_be(&b60, sizeof(double), 1, fich) ;
369 fwrite_be(&ent0, sizeof(double), 1, fich) ;
370 fwrite_be(&eps_fit, sizeof(double), 1, fich) ;
371 fwrite_be(&rho0_b60, sizeof(double), 1, fich) ;
372 fwrite_be(&ent_nd, sizeof(double), 1, fich) ;
373 fwrite_be(&rho_nd, sizeof(double), 1, fich) ;
374 fwrite_be(&gam, sizeof(double), 1, fich) ;
375
376}
377
378ostream& Eos_strange_cr::operator>>(ostream & ost) const {
379
380 ost <<
381 "EOS of class Eos_strange_cr " << endl
382 << " (Strange matter EOS from Zdunik (2000) with crust) "
383 << endl ;
384 ost << " Baryon density at zero pressure : " << n0_b60
385 << " * B_{60}^{3/4}" << endl ;
386 ost << " Bag constant B : " << b60 << " * 60 MeV/fm^3"<< endl ;
387 ost <<
388 " Log-enthalpy threshold for setting the energy density to non-zero: "
389 << endl << " " << ent0 << endl ;
390 ost << " Fitting parameter eps_fit : " << eps_fit << endl ;
391 ost << " Energy density at zero pressure : " << rho0_b60
392 << " * B_{60} MeV/fm^3" << endl ;
393 ost << " Log-enthalpy at neutron drip point : " << ent_nd << endl ;
394 ost << " Energy density at neutron drip point : " << rho_nd
395 << " MeV/fm^3" << endl ;
396 ost << " Adiabatic index for the crust : " << gam << endl ;
397
398 return ost ;
399
400}
401
402
403 //------------------------------//
404 // Computational routines //
405 //------------------------------//
406
407// Baryon density from enthalpy
408//------------------------------
409
410double Eos_strange_cr::nbar_ent_p(double ent, const Param* ) const {
411
412 if ( ent > ent0 ) {
413
414 if (ent > ent_nd) { // Strange quark matter
415
416 double entsqm = ent + delent ;
417
418 return n0 * exp( double(3) * entsqm / (double(1) + eps_fit)) ;
419 }
420 else { // crust
421
422 double fn_au = pow( gam1sx * ( exp(ent)- double(1) ), unsgam1) ;
423
424 return ncr_nd * fn_au ;
425
426 }
427
428 }
429 else{
430 return 0 ;
431 }
432}
433
434// Energy density from enthalpy
435//------------------------------
436
437double Eos_strange_cr::ener_ent_p(double ent, const Param* ) const {
438
439
440 if ( ent > ent0 ) {
441
442 if (ent > ent_nd) { // Strange quark matter
443
444 double entsqm = ent + delent ;
445
446 double pp = ( exp(fach * entsqm) - double(1)) / fach * rho0 ;
447
448 return rho0 + double(3) * pp / (double(1) + eps_fit) ;
449
450 }
451 else { // crust
452
453 double fn_au = pow( gam1sx *(exp(ent) - double(1) ), unsgam1) ;
454
455 return rho_nd_nucl * ( (double(1) - x_nd * unsgam1 ) * fn_au
456 + x_nd * unsgam1 * pow(fn_au, gam)) ;
457 }
458
459 }
460 else{
461 return 0 ;
462 }
463}
464
465// Pressure from enthalpy
466//------------------------
467
468double Eos_strange_cr::press_ent_p(double ent, const Param* ) const {
469
470 if ( ent > ent0 ) {
471
472 if (ent > ent_nd) { // Strange quark matter
473
474 double entsqm = ent + delent ;
475
476 return ( exp(fach * entsqm) - double(1) ) / fach * rho0 ;
477 }
478 else{ // crust
479
480 double fn_au = pow( gam1sx *(exp(ent) - double(1) ), unsgam1) ;
481
482 return x_nd * rho_nd_nucl * pow(fn_au, gam) ;
483
484 }
485
486
487 }
488 else{
489 return 0 ;
490 }
491}
492
493
494// dln(n)/ln(H) from enthalpy
495//---------------------------
496
497double Eos_strange_cr::der_nbar_ent_p(double ent, const Param* ) const {
498
499 if ( ent > ent0 ) {
500
501 if (ent > ent_nd) { // Strange quark matter
502
503 double entsqm = ent + delent ;
504
505 return double(3) * entsqm / ( double(1) + eps_fit ) ;
506 }
507 else{ // crust
508
509 cout << "Eos_strange_cr::der_nbar_ent_p not ready yet !" << endl ;
510 abort() ;
511 return 0 ;
512 }
513
514 }
515 else{
516 return 0 ;
517 }
518}
519
520// dln(e)/ln(H) from enthalpy
521//---------------------------
522
523double Eos_strange_cr::der_ener_ent_p(double ent, const Param* ) const {
524
525 if ( ent > ent0 ) {
526
527 if (ent > ent_nd) { // Strange quark matter
528
529 double entsqm = ent + delent ;
530
531 double xx = fach * entsqm ;
532
533 return xx / ( double(1) +
534 ( double(1) + eps_fit ) / double(3) * exp(-xx) ) ;
535 }
536 else{ // crust
537
538 cout << "Eos_strange_cr::der_ener_ent_p not ready yet !" << endl ;
539 abort() ;
540 return 0 ;
541 }
542 }
543 else{
544 return 0 ;
545 }
546}
547
548// dln(p)/ln(H) from enthalpy
549//---------------------------
550
551double Eos_strange_cr::der_press_ent_p(double ent, const Param* ) const {
552
553 if ( ent > ent0 ) {
554
555
556 if (ent > ent_nd) { // Strange quark matter
557
558 double entsqm = ent + delent ;
559
560 double xx = fach * entsqm ;
561
562 return xx / ( double(1) - exp(-xx) ) ;
563 }
564 else{ // crust
565
566 cout << "Eos_strange_cr::der_press_ent_p not ready yet !" << endl ;
567 abort() ;
568 return 0 ;
569 }
570
571 }
572 else{
573 return 0 ;
574 }
575}
576
577
578
579
580}
Strange matter EOS (MIT Bag model) with crust.
Definition eos.h:1779
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
double delent
Enthalpy shift in quark phase.
Definition eos.h:1867
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double rho0
Energy density at zero pressure.
Definition eos.h:1839
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double rho_nd_nucl
Energy density at neutron drip point, defining the boundary between crust and core [unit: rho_unit ].
Definition eos.h:1856
virtual ~Eos_strange_cr()
Destructor.
double rho_nd
Energy density at neutron drip point, defining the boundary between crust and core [unit: ].
Definition eos.h:1820
double rho0_b60
Energy density at zero pressure divided by .
Definition eos.h:1807
double ent0
Log-enthalpy threshold for setting the energy density to a non zero value (should be negative).
Definition eos.h:1796
double ent_nd
Log-enthalpy at neutron drip point, defining the boundary between crust and core.
Definition eos.h:1813
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
double n0_b60
Baryon density at zero pressure divided by .
Definition eos.h:1788
double b60
Bag constant [unit: ].
Definition eos.h:1791
Eos_strange_cr(double n0_b60_i, double b60_i, double ent0_i, double eps_fit_i, double rho0_b60_i, double ent_nd_i, double rho_nd_i, double gam_i)
Standard constructor.
double ncr_nd
Rescaled number density at neutron drip point.
Definition eos.h:1864
double x_nd
Ratio of pressure to energy density at neutron drip point.
Definition eos.h:1861
virtual ostream & operator>>(ostream &) const
Operator >>
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void operator=(const Eos_strange_cr &)
Assignment to another Eos_strange.
void set_auxiliary()
Computes the auxiliary quantities n0 , rh0 , b34 and fach from the values of the other parameters.
virtual bool operator==(const Eos &) const
Comparison operator (egality)
virtual void sauve(FILE *) const
Save in a file.
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
double eps_fit
Fitting parameter related to the square of sound velocity by .
Definition eos.h:1802
double n0
Baryon density at zero pressure.
Definition eos.h:1833
double fach
Factor .
Definition eos.h:1849
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
double gam
Adiabatic index for the crust model.
Definition eos.h:1825
Equation of state base class.
Definition eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual void sauve(FILE *) const
Save in a file.
Definition eos.C:179
char name[100]
EOS name.
Definition eos.h:196
void set_name(const char *name_i)
Sets the EOS name.
Definition eos.C:163
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.