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