LORENE
eos_multi_poly.C
1 /*
2  * Methods of class Eos_multi_poly
3  *
4  * (see file eos_multi_poly.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2009 Keisuke Taniguchi
10  * Copyright (c) 2004 Keisuke Taniguchi
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 version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 char eos_multi_poly_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_multi_poly.C,v 1.10 2014/12/09 14:07:14 j_novak Exp $" ;
30 
31 /*
32  * $Id: eos_multi_poly.C,v 1.10 2014/12/09 14:07:14 j_novak Exp $
33  * $Log: eos_multi_poly.C,v $
34  * Revision 1.10 2014/12/09 14:07:14 j_novak
35  * Changed (corrected?) the formula for computing the kappa's.
36  *
37  * Revision 1.9 2014/10/13 08:52:53 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.8 2014/10/06 15:13:06 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.7 2012/10/09 16:15:26 j_novak
44  * Corrected a bug in the constructor and save into a file
45  *
46  * Revision 1.6 2009/06/23 14:34:04 k_taniguchi
47  * Completely revised.
48  *
49  * Revision 1.5 2004/06/23 15:42:08 e_gourgoulhon
50  * Replaced all "abs" by "fabs".
51  *
52  * Revision 1.4 2004/05/13 07:38:57 k_taniguchi
53  * Change the procedure for searching the baryon density from enthalpy.
54  *
55  * Revision 1.3 2004/05/09 10:43:52 k_taniguchi
56  * Change the searching method of the baryon density again.
57  *
58  * Revision 1.2 2004/05/07 11:55:59 k_taniguchi
59  * Change the searching procedure of the baryon density.
60  *
61  * Revision 1.1 2004/05/07 08:10:58 k_taniguchi
62  * Initial revision
63  *
64  *
65  *
66  * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_multi_poly.C,v 1.10 2014/12/09 14:07:14 j_novak Exp $
67  *
68  */
69 
70 // C headers
71 #include <cstdlib>
72 #include <cstring>
73 #include <cmath>
74 
75 // Lorene headers
76 #include "headcpp.h"
77 #include "eos_multi_poly.h"
78 #include "eos.h"
79 #include "utilitaires.h"
80 #include "unites.h"
81 
82 namespace Lorene {
83 double logp(double, double, double, double, double, double) ;
84 double dlpsdlh(double, double, double, double, double, double) ;
85 double dlpsdlnb(double, double, double, double, double) ;
86 
87 //*********************************************************************
88 
89  //--------------------------------------//
90  // Constructors //
91  //--------------------------------------//
92 
93 // Standard constructor
94 Eos_multi_poly::Eos_multi_poly(int npoly, double* gamma_i, double kappa0_i,
95  double logP1_i, double* logRho_i,
96  double* decInc_i)
97  : Eos("Multi-polytropic EOS"),
98  npeos(npoly), kappa0(kappa0_i), logP1(logP1_i), m0(double(1)) {
99 
100  assert(npeos > 1) ;
101 
102  gamma = new double [npeos] ;
103 
104  for (int l=0; l<npeos; l++) {
105  gamma[l] = gamma_i[l] ;
106  }
107 
108  logRho = new double [npeos-1] ;
109 
110  for (int l=0; l<npeos-1; l++) {
111  logRho[l] = logRho_i[l] ;
112  }
113 
114  decInc = new double [npeos-1] ;
115 
116  for (int l=0; l<npeos-1; l++) {
117  decInc[l] = decInc_i[l] ;
118  }
119 
120  set_auxiliary() ;
121 
122 }
123 
124 
125 // Copy constructor
127  : Eos(eosmp),
128  npeos(eosmp.npeos), kappa0(eosmp.kappa0), logP1(eosmp.logP1),
129  m0(eosmp.m0) {
130 
131  gamma = new double [npeos] ;
132 
133  for (int l=0; l<npeos; l++) {
134  gamma[l] = eosmp.gamma[l] ;
135  }
136 
137  logRho = new double [npeos-1] ;
138 
139  for (int l=0; l<npeos-1; l++) {
140  logRho[l] = eosmp.logRho[l] ;
141  }
142 
143  kappa = new double [npeos] ;
144 
145  for (int l=0; l<npeos; l++) {
146  kappa[l] = eosmp.kappa[l] ;
147  }
148 
149  nbCrit = new double [npeos-1] ;
150 
151  for (int l=0; l<npeos-1; l++) {
152  nbCrit[l] = eosmp.nbCrit[l] ;
153  }
154 
155  entCrit = new double [npeos-1] ;
156 
157  for (int l=0; l<npeos-1; l++) {
158  entCrit[l] = eosmp.entCrit[l] ;
159  }
160 
161  decInc = new double [npeos-1] ;
162 
163  for (int l=0; l<npeos-1; l++) {
164  decInc[l] = eosmp.decInc[l] ;
165  }
166 
167  mu0 = new double [npeos] ;
168 
169  for (int l=0; l<npeos; l++) {
170  mu0[l] = eosmp.mu0[l] ;
171  }
172 
173 }
174 
175 // Constructor from a binary file
176 Eos_multi_poly::Eos_multi_poly(FILE* fich) : Eos(fich) {
177 
178  fread_be(&npeos, sizeof(int), 1, fich) ;
179 
180  gamma = new double [npeos] ;
181 
182  for (int l=0; l<npeos; l++) {
183  fread_be(&gamma[l], sizeof(double), 1, fich) ;
184  }
185 
186  fread_be(&kappa0, sizeof(double), 1, fich) ;
187  fread_be(&logP1, sizeof(double), 1, fich) ;
188 
189  logRho = new double [npeos-1] ;
190 
191  for (int l=0; l<npeos-1; l++) {
192  fread_be(&logRho[l], sizeof(double), 1, fich) ;
193  }
194 
195  decInc = new double [npeos-1] ;
196 
197  for (int l=0; l<npeos-1; l++) {
198  fread_be(&decInc[l], sizeof(double), 1, fich) ;
199  }
200 
201  m0 = double(1) ;
202 
203  set_auxiliary() ;
204 
205 }
206 
207 // Constructor from a formatted file
208 Eos_multi_poly::Eos_multi_poly(ifstream& fich) : Eos(fich) {
209 
210  char blabla[80] ;
211 
212  fich >> npeos ; fich.getline(blabla, 80) ;
213 
214  gamma = new double [npeos] ;
215 
216  for (int l=0; l<npeos; l++) {
217  fich >> gamma[l] ; fich.getline(blabla, 80) ;
218  }
219 
220  fich >> kappa0 ; fich.getline(blabla, 80) ;
221  fich >> logP1 ; fich.getline(blabla, 80) ;
222 
223  logRho = new double [npeos-1] ;
224 
225  for (int l=0; l<npeos-1; l++) {
226  fich >> logRho[l] ; fich.getline(blabla, 80) ;
227  }
228 
229  decInc = new double [npeos-1] ;
230 
231  for (int l=0; l<npeos-1; l++) {
232  fich >> decInc[l] ; fich.getline(blabla, 80) ;
233  }
234 
235  m0 = double(1) ;
236 
237  set_auxiliary() ;
238 
239 }
240 
241 // Destructor
243 
244  delete [] gamma ;
245  delete [] logRho ;
246  delete [] kappa ;
247  delete [] nbCrit ;
248  delete [] entCrit ;
249  delete [] decInc ;
250  delete [] mu0 ;
251 
252 }
253 
254  //--------------//
255  // Assignment //
256  //--------------//
257 
259 
260  cout << "Eos_multi_poly::operator= : not implemented yet !" << endl ;
261  abort() ;
262 
263 }
264 
265 
266  //-----------------------//
267  // Miscellaneous //
268  //-----------------------//
269 
271 
272  using namespace Unites ;
273 
274  double* kappa_cgs = new double [npeos] ;
275 
276  kappa_cgs[0] = kappa0 ;
277 
278  kappa_cgs[1] = pow(10., logP1-logRho[0]*gamma[1]) ;
279 
280  if (npeos > 2) {
281 
282  kappa_cgs[2] = kappa_cgs[1]
283  * pow(10., logRho[1]*(gamma[1]-gamma[2])) ;
284 
285  if (npeos > 3) {
286 
287  for (int l=3; l<npeos; l++) {
288 
289  kappa_cgs[l] = kappa_cgs[l-1]
290  * pow(10., logRho[l-1]*(gamma[l-1]-gamma[l])) ;
291 
292  }
293 
294  }
295 
296  }
297 
298  kappa = new double [npeos] ;
299 
300  double rhonuc_cgs = rhonuc_si * 1.e-3 ;
301 
302  for (int l=0; l<npeos; l++) {
303  kappa[l] = kappa_cgs[l] * pow( rhonuc_cgs, gamma[l] - double(1) ) ;
304  // Conversion from cgs units to Lorene units
305  }
306 
307  delete [] kappa_cgs ;
308 
309  mu0 = new double [npeos] ;
310  mu0[0] = double(1) ; // We define
311 
312  entCrit = new double [npeos-1] ;
313 
314  nbCrit = new double [npeos-1] ;
315 
316  for (int l=0; l<npeos-1; l++) {
317 
318  nbCrit[l] =
319  pow(kappa[l]/kappa[l+1], double(1)/(gamma[l+1]-gamma[l])) ;
320 
321  mu0[l+1] = mu0[l]
322  + ( kappa[l] * pow(nbCrit[l], gamma[l]-double(1))
323  / (gamma[l]-double(1))
324  - kappa[l+1] * pow(nbCrit[l], gamma[l+1]-double(1))
325  / (gamma[l+1]-double(1)) ) ;
326 
327  entCrit[l] = log ( mu0[l] / m0
328  + kappa[l] * gamma[l]
329  * pow(nbCrit[l], gamma[l]-double(1))
330  / (gamma[l]-double(1)) / m0 ) ;
331 
332  }
333 
334 }
335 
336 
337  //------------------------------//
338  // Comparison operators //
339  //------------------------------//
340 
341 bool Eos_multi_poly::operator==(const Eos& eos_i) const {
342 
343  bool resu = true ;
344 
345  if ( eos_i.identify() != identify() ) {
346  cout << "The second EOS is not of type Eos_multi_poly !" << endl ;
347  resu = false ;
348  }
349  else{
350 
351  const Eos_multi_poly& eos
352  = dynamic_cast<const Eos_multi_poly&>(eos_i) ;
353 
354  if (eos.get_npeos() != npeos) {
355  cout << "The two Eos_multi_poly have "
356  << "different number of polytropes : "
357  << npeos << " <-> " << eos.get_npeos() << endl ;
358  resu = false ;
359  }
360 
361  for (int l=0; l<npeos; l++) {
362  if (eos.get_gamma(l) != gamma[l]) {
363  cout << "The two Eos_multi_poly have different gamma "
364  << gamma[l] << " <-> "
365  << eos.get_gamma(l) << endl ;
366  resu = false ;
367  }
368  }
369 
370  for (int l=0; l<npeos; l++) {
371  if (eos.get_kappa(l) != kappa[l]) {
372  cout << "The two Eos_multi_poly have different kappa "
373  << kappa[l] << " <-> "
374  << eos.get_kappa(l) << endl ;
375  resu = false ;
376  }
377  }
378 
379  }
380 
381  return resu ;
382 
383 }
384 
385 bool Eos_multi_poly::operator!=(const Eos& eos_i) const {
386 
387  return !(operator==(eos_i)) ;
388 
389 }
390 
391  //--------------------------//
392  // Outputs //
393  //--------------------------//
394 
395 void Eos_multi_poly::sauve(FILE* fich) const {
396 
397  Eos::sauve(fich) ;
398 
399  fwrite_be(&npeos, sizeof(int), 1, fich) ;
400 
401  for (int l=0; l<npeos; l++) {
402  fwrite_be(&gamma[l], sizeof(double), 1, fich) ;
403  }
404 
405  fwrite_be(&kappa0, sizeof(double), 1, fich) ;
406  fwrite_be(&logP1, sizeof(double), 1, fich) ;
407 
408  for (int l=0; l<npeos-1; l++) {
409  fwrite_be(&logRho[l], sizeof(double), 1, fich) ;
410  }
411 
412  for (int l=0; l<npeos-1; l++) {
413  fwrite_be(&decInc[l], sizeof(double), 1, fich) ;
414  }
415 
416 }
417 
418 
419 ostream& Eos_multi_poly::operator>>(ostream & ost) const {
420 
421  using namespace Unites ;
422 
423  ost << "EOS of class Eos_multi_poly "
424  << "(multiple polytropic equation of state) : " << endl ;
425 
426  ost << " Number of polytropes : "
427  << npeos << endl << endl ;
428 
429  double rhonuc_cgs = rhonuc_si * 1.e-3 ;
430 
431  ost.precision(16) ;
432  for (int l=0; l<npeos; l++) {
433  ost << " EOS in region " << l << " : " << endl ;
434  ost << " ---------------" << endl ;
435  ost << " gamma : " << gamma[l] << endl ;
436  ost << " kappa : " << kappa[l]
437  << " [Lorene units: rho_nuc c^2 / n_nuc^gamma]" << endl ;
438 
439  double kappa_cgs = kappa[l]
440  * pow( rhonuc_cgs, double(1) - gamma[l] ) ;
441 
442  ost << " : " << kappa_cgs
443  << " [(g/cm^3)^{1-gamma}]" << endl ;
444  }
445 
446  ost << endl ;
447  ost << " Exponent of the pressure at the fiducial density rho_1"
448  << endl ;
449  ost << " ------------------------------------------------------"
450  << endl ;
451  ost << " log P1 : " << logP1 << endl ;
452 
453  ost << endl ;
454  ost << " Exponent of fiducial densities" << endl ;
455  ost << " ------------------------------" << endl ;
456  for (int l=0; l<npeos-1; l++) {
457  ost << " log rho[" << l << "] : " << logRho[l] << endl ;
458  }
459 
460  ost << endl ;
461  for (int l=0; l<npeos-1; l++) {
462  ost << " Critical density and enthalpy between domains "
463  << l << " and " << l+1 << " : " << endl ;
464  ost << " -----------------------------------------------------"
465  << endl ;
466  ost << " num. dens. : " << nbCrit[l] << " [Lorene units: n_nuc]"
467  << endl ;
468  ost << " density : " << nbCrit[l] * rhonuc_cgs << " [g/cm^3]"
469  << endl ;
470 
471  ost << " ln(ent) : " << entCrit[l] << endl ;
472  }
473 
474  ost << endl ;
475  for (int l=0; l<npeos; l++) {
476  ost << " Relat. chem. pot. in region " << l << " : " << endl ;
477  ost << " -----------------------------" << endl ;
478  ost << " mu : " << mu0[l] << " [m_B c^2]" << endl ;
479  }
480 
481  return ost ;
482 
483 }
484 
485 
486  //------------------------------------------//
487  // Computational routines //
488  //------------------------------------------//
489 
490 // Baryon rest-mass density from enthalpy
491 //----------------------------------------
492 
493 double Eos_multi_poly::nbar_ent_p(double ent, const Param* ) const {
494 
495  int i = 0 ; // "i" corresponds to the number of domain
496  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
497  // The buffer zone is included in the next zone.
498  for (int l=0; l<npeos-1; l++) {
499  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
500  i++ ;
501  }
502  }
503 
504  double mgam1skapgam = 1. ; // Initialization
505  if (i == 0) {
506 
507  if ( ent > double(0) ) {
508 
509  mgam1skapgam = m0*(gamma[0]-double(1))/kappa[0]/gamma[0] ;
510 
511  return pow( mgam1skapgam*(exp(ent)-double(1)),
512  double(1)/(gamma[0]-double(1)) ) ; // mu0[0]/m0=1
513 
514  }
515  else {
516  return double(0) ;
517  }
518 
519  }
520  else {
521 
522  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
523  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
524 
525  if ( ent < entLarge ) {
526 
527  double log10H = log10( ent ) ;
528  double log10HSmall = log10( entSmall ) ;
529  double log10HLarge = log10( entLarge ) ;
530  double dH = log10HLarge - log10HSmall ;
531  double uu = (log10H - log10HSmall) / dH ;
532 
533  double mgam1skapgamSmall = m0*(gamma[i-1]-double(1))
534  /kappa[i-1]/gamma[i-1] ;
535  double mgam1skapgamLarge = m0*(gamma[i]-double(1))
536  /kappa[i]/gamma[i] ;
537 
538  double nnSmall = pow( mgam1skapgamSmall
539  *(exp(entSmall)-mu0[i-1]/m0),
540  double(1)/(gamma[i-1]-double(1)) ) ;
541  double nnLarge = pow( mgam1skapgamLarge
542  *(exp(entLarge)-mu0[i]/m0),
543  double(1)/(gamma[i]-double(1)) ) ;
544 
545  double ppSmall = kappa[i-1] * pow( nnSmall, gamma[i-1] ) ;
546  double ppLarge = kappa[i] * pow( nnLarge, gamma[i] ) ;
547 
548  double eeSmall = mu0[i-1] * nnSmall
549  + ppSmall / (gamma[i-1] - double(1)) ;
550  double eeLarge = mu0[i] * nnLarge
551  + ppLarge / (gamma[i] - double(1)) ;
552 
553  double log10PSmall = log10( ppSmall ) ;
554  double log10PLarge = log10( ppLarge ) ;
555 
556  double dlpsdlhSmall = entSmall
557  * (double(1) + eeSmall / ppSmall) ;
558  double dlpsdlhLarge = entLarge
559  * (double(1) + eeLarge / ppLarge) ;
560 
561  double log10PInterpolate = logp(log10PSmall, log10PLarge,
562  dlpsdlhSmall, dlpsdlhLarge,
563  dH, uu) ;
564 
565  double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
566  dlpsdlhSmall, dlpsdlhLarge,
567  dH, uu) ;
568 
569  double pp = pow(double(10), log10PInterpolate) ;
570 
571  return pp / ent * dlpsdlhInterpolate * exp(-ent) / m0 ;
572  // Is m0 necessary?
573 
574  }
575  else {
576 
577  mgam1skapgam = m0*(gamma[i]-double(1))/kappa[i]/gamma[i] ;
578 
579  return pow( mgam1skapgam*(exp(ent)-mu0[i]/m0),
580  double(1)/(gamma[i]-double(1)) ) ;
581 
582  }
583 
584  }
585 
586 }
587 
588 // Energy density from enthalpy
589 //------------------------------
590 
591 double Eos_multi_poly::ener_ent_p(double ent, const Param* ) const {
592 
593  int i = 0 ; // "i" corresponds to the number of domain
594  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
595  // The buffer zone is included in the next zone.
596  for (int l=0; l<npeos-1; l++) {
597  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
598  i++ ;
599  }
600  }
601 
602  double mgam1skapgam = 1. ; // Initialization
603  double nn = 0. ; // Initialization
604  double pp = 0. ; // Initialization
605 
606  if (i == 0) {
607 
608  if ( ent > double(0) ) {
609 
610  mgam1skapgam = m0*(gamma[0]-double(1))/kappa[0]/gamma[0] ;
611 
612  nn = pow( mgam1skapgam*(exp(ent)-double(1)),
613  double(1)/(gamma[0]-double(1)) ) ; // mu0[0]/m0=1
614 
615  pp = kappa[0] * pow( nn, gamma[0] ) ;
616 
617  return mu0[0] * nn + pp / (gamma[0] - double(1)) ;
618 
619  }
620  else {
621  return double(0) ;
622  }
623 
624  }
625  else {
626 
627  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
628  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
629 
630  if ( ent < entLarge ) {
631 
632  double log10H = log10( ent ) ;
633  double log10HSmall = log10( entSmall ) ;
634  double log10HLarge = log10( entLarge ) ;
635  double dH = log10HLarge - log10HSmall ;
636  double uu = (log10H - log10HSmall) / dH ;
637 
638  double mgam1skapgamSmall = m0*(gamma[i-1]-double(1))
639  /kappa[i-1]/gamma[i-1] ;
640  double mgam1skapgamLarge = m0*(gamma[i]-double(1))
641  /kappa[i]/gamma[i] ;
642 
643  double nnSmall = pow( mgam1skapgamSmall
644  *(exp(entSmall)-mu0[i-1]/m0),
645  double(1)/(gamma[i-1]-double(1)) ) ;
646  double nnLarge = pow( mgam1skapgamLarge
647  *(exp(entLarge)-mu0[i]/m0),
648  double(1)/(gamma[i]-double(1)) ) ;
649 
650  double ppSmall = kappa[i-1] * pow( nnSmall, gamma[i-1] ) ;
651  double ppLarge = kappa[i] * pow( nnLarge, gamma[i] ) ;
652 
653  double eeSmall = mu0[i-1] * nnSmall
654  + ppSmall / (gamma[i-1] - double(1)) ;
655  double eeLarge = mu0[i] * nnLarge
656  + ppLarge / (gamma[i] - double(1)) ;
657 
658  double log10PSmall = log10( ppSmall ) ;
659  double log10PLarge = log10( ppLarge ) ;
660 
661  double dlpsdlhSmall = entSmall
662  * (double(1) + eeSmall / ppSmall) ;
663  double dlpsdlhLarge = entLarge
664  * (double(1) + eeLarge / ppLarge) ;
665 
666  double log10PInterpolate = logp(log10PSmall, log10PLarge,
667  dlpsdlhSmall, dlpsdlhLarge,
668  dH, uu) ;
669 
670  double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
671  dlpsdlhSmall, dlpsdlhLarge,
672  dH, uu) ;
673 
674  pp = pow(double(10), log10PInterpolate) ;
675 
676  return pp / ent * dlpsdlhInterpolate - pp ;
677 
678  }
679  else {
680 
681  mgam1skapgam = m0*(gamma[i]-double(1))/kappa[i]/gamma[i] ;
682 
683  nn = pow( mgam1skapgam*(exp(ent)-mu0[i]/m0),
684  double(1)/(gamma[i]-double(1)) ) ;
685 
686  pp = kappa[i] * pow( nn, gamma[i] ) ;
687 
688  return mu0[i] * nn + pp / (gamma[i] - double(1)) ;
689 
690  }
691 
692  }
693 
694 }
695 
696 
697 // Pressure from enthalpy
698 //------------------------
699 
700 double Eos_multi_poly::press_ent_p(double ent, const Param* ) const {
701 
702  int i = 0 ; // "i" corresponds to the number of domain
703  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
704  // The buffer zone is included in the next zone.
705  for (int l=0; l<npeos-1; l++) {
706  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
707  i++ ;
708  }
709  }
710 
711  double mgam1skapgam = 1. ; // Initialization
712  double nn = 0. ; // Initialization
713 
714  if (i == 0) {
715 
716  if ( ent > double(0) ) {
717 
718  mgam1skapgam = m0*(gamma[0]-double(1))/kappa[0]/gamma[0] ;
719 
720  nn = pow( mgam1skapgam*(exp(ent)-double(1)),
721  double(1)/(gamma[0]-double(1)) ) ; // mu0[0]/m0=1
722 
723  return kappa[0] * pow( nn, gamma[0] ) ;
724 
725  }
726  else {
727  return double(0) ;
728  }
729 
730  }
731  else {
732 
733  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
734  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
735 
736  if ( ent < entLarge ) {
737 
738  double log10H = log10( ent ) ;
739  double log10HSmall = log10( entSmall ) ;
740  double log10HLarge = log10( entLarge ) ;
741  double dH = log10HLarge - log10HSmall ;
742  double uu = (log10H - log10HSmall) / dH ;
743 
744  double mgam1skapgamSmall = m0*(gamma[i-1]-double(1))
745  /kappa[i-1]/gamma[i-1] ;
746  double mgam1skapgamLarge = m0*(gamma[i]-double(1))
747  /kappa[i]/gamma[i] ;
748 
749  double nnSmall = pow( mgam1skapgamSmall
750  *(exp(entSmall)-mu0[i-1]/m0),
751  double(1)/(gamma[i-1]-double(1)) ) ;
752  double nnLarge = pow( mgam1skapgamLarge
753  *(exp(entLarge)-mu0[i]/m0),
754  double(1)/(gamma[i]-double(1)) ) ;
755 
756  double ppSmall = kappa[i-1] * pow( nnSmall, gamma[i-1] ) ;
757  double ppLarge = kappa[i] * pow( nnLarge, gamma[i] ) ;
758 
759  double eeSmall = mu0[i-1] * nnSmall
760  + ppSmall / (gamma[i-1] - double(1)) ;
761  double eeLarge = mu0[i] * nnLarge
762  + ppLarge / (gamma[i] - double(1)) ;
763 
764  double log10PSmall = log10( ppSmall ) ;
765  double log10PLarge = log10( ppLarge ) ;
766 
767  double dlpsdlhSmall = entSmall
768  * (double(1) + eeSmall / ppSmall) ;
769  double dlpsdlhLarge = entLarge
770  * (double(1) + eeLarge / ppLarge) ;
771 
772  double log10PInterpolate = logp(log10PSmall, log10PLarge,
773  dlpsdlhSmall, dlpsdlhLarge,
774  dH, uu) ;
775 
776  return pow(double(10), log10PInterpolate) ;
777 
778  }
779  else {
780 
781  mgam1skapgam = m0*(gamma[i]-double(1))/kappa[i]/gamma[i] ;
782 
783  nn = pow( mgam1skapgam*(exp(ent)-mu0[i]/m0),
784  double(1)/(gamma[i]-double(1)) ) ;
785 
786  return kappa[i] * pow( nn, gamma[i] ) ;
787 
788  }
789 
790  }
791 
792 }
793 
794 
795 // dln(n)/dln(H) from enthalpy
796 //----------------------------
797 
798 double Eos_multi_poly::der_nbar_ent_p(double ent, const Param* ) const {
799 
800  int i = 0 ; // "i" corresponds to the number of domain
801  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
802  // The buffer zone is included in the next zone.
803  for (int l=0; l<npeos-1; l++) {
804  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
805  i++ ;
806  }
807  }
808 
809  if (i == 0) {
810 
811  if ( ent > double(0) ) {
812 
813  if ( ent < 1.e-13 ) {
814 
815  return ( double(1) + ent/double(2) + ent*ent/double(12) )
816  / (gamma[0] - double(1)) ;
817 
818  }
819  else {
820 
821  return ent / (double(1) - exp(-ent))
822  / (gamma[0] - double(1)) ; // mu0[0]/m0=1
823 
824  }
825 
826  }
827  else {
828 
829  return double(1) / (gamma[0] - double(1)) ;
830 
831  }
832 
833  }
834  else {
835 
836  if ( ent < entCrit[i-1]*(double(1)+decInc[i-1]) ) {
837 
838  double zeta = der_press_ent_p(ent) / der_press_nbar_p(ent) ;
839 
840  return zeta ;
841 
842  }
843  else {
844 
845  return ent / (double(1) - (mu0[i]/m0) * exp(-ent))
846  / (gamma[i] - double(1)) ;
847 
848  }
849 
850  }
851 }
852 
853 
854 // dln(e)/dln(H) from enthalpy
855 //----------------------------
856 
857 double Eos_multi_poly::der_ener_ent_p(double ent, const Param* ) const {
858 
859  int i = 0 ; // "i" corresponds to the number of domain
860  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
861  // The buffer zone is included in the next zone.
862  for (int l=0; l<npeos-1; l++) {
863  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
864  i++ ;
865  }
866  }
867 
868  double mgam1skapgam = 1. ; // Initialization
869  double nn = 0. ; // Initialization
870  double pp = 0. ; // Initialization
871  double ee = 0. ; // Initialization
872 
873  if (i == 0) {
874 
875  if ( ent > double(0) ) {
876 
877  mgam1skapgam = m0*(gamma[0]-double(1))/kappa[0]/gamma[0] ;
878 
879  nn = pow( mgam1skapgam*(exp(ent)-double(1)),
880  double(1)/(gamma[0]-double(1)) ) ; // mu0[0]/m0=1
881 
882  pp = kappa[0] * pow( nn, gamma[0] ) ;
883 
884  ee = mu0[0] * nn + pp / (gamma[0] - double(1)) ;
885 
886  if ( ent < 1.e-13 ) {
887 
888  return ( double(1) + ent/double(2) + ent*ent/double(12) )
889  / (gamma[0] - double(1)) * (double(1) + pp / ee) ;
890 
891  }
892  else {
893 
894  return ent / (double(1) - exp(-ent))
895  / (gamma[0] - double(1)) * (double(1) + pp / ee) ;
896  // mu0[0]/m0=1
897 
898  }
899 
900  }
901  else {
902 
903  return double(1) / (gamma[0] - double(1)) ;
904 
905  }
906 
907  }
908  else {
909 
910  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
911  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
912 
913  if ( ent < entLarge ) {
914 
915  double log10H = log10( ent ) ;
916  double log10HSmall = log10( entSmall ) ;
917  double log10HLarge = log10( entLarge ) ;
918  double dH = log10HLarge - log10HSmall ;
919  double uu = (log10H - log10HSmall) / dH ;
920 
921  double mgam1skapgamSmall = m0*(gamma[i-1]-double(1))
922  /kappa[i-1]/gamma[i-1] ;
923  double mgam1skapgamLarge = m0*(gamma[i]-double(1))
924  /kappa[i]/gamma[i] ;
925 
926  double nnSmall = pow( mgam1skapgamSmall
927  *(exp(entSmall)-mu0[i-1]/m0),
928  double(1)/(gamma[i-1]-double(1)) ) ;
929  double nnLarge = pow( mgam1skapgamLarge
930  *(exp(entLarge)-mu0[i]/m0),
931  double(1)/(gamma[i]-double(1)) ) ;
932 
933  double ppSmall = kappa[i-1] * pow( nnSmall, gamma[i-1] ) ;
934  double ppLarge = kappa[i] * pow( nnLarge, gamma[i] ) ;
935 
936  double eeSmall = mu0[i-1] * nnSmall
937  + ppSmall / (gamma[i-1] - double(1)) ;
938  double eeLarge = mu0[i] * nnLarge
939  + ppLarge / (gamma[i] - double(1)) ;
940 
941  double log10PSmall = log10( ppSmall ) ;
942  double log10PLarge = log10( ppLarge ) ;
943 
944  double dlpsdlhSmall = entSmall
945  * (double(1) + eeSmall / ppSmall) ;
946  double dlpsdlhLarge = entLarge
947  * (double(1) + eeLarge / ppLarge) ;
948 
949  double log10PInterpolate = logp(log10PSmall, log10PLarge,
950  dlpsdlhSmall, dlpsdlhLarge,
951  dH, uu) ;
952 
953  double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
954  dlpsdlhSmall, dlpsdlhLarge,
955  dH, uu) ;
956 
957  pp = pow(double(10), log10PInterpolate) ;
958 
959  ee = pp / ent * dlpsdlhInterpolate - pp ;
960 
961  double zeta = (double(1) + pp / ee) * der_press_ent_p(ent)
962  / der_press_nbar_p(ent) ;
963 
964  return zeta ;
965 
966  }
967  else {
968 
969  mgam1skapgam = m0*(gamma[i]-double(1))/kappa[i]/gamma[i] ;
970 
971  nn = pow( mgam1skapgam*(exp(ent)-mu0[i]/m0),
972  double(1)/(gamma[i]-double(1)) ) ;
973 
974  pp = kappa[i] * pow( nn, gamma[i] ) ;
975 
976  ee = mu0[i] * nn + pp / (gamma[i] - double(1)) ;
977 
978  return ent / (double(1) - (mu0[i]/m0) * exp(-ent))
979  / (gamma[i] - double(1)) * (double(1) + pp / ee) ;
980 
981  }
982 
983  }
984 }
985 
986 
987 // dln(p)/dln(H) from enthalpy
988 //----------------------------
989 
990 double Eos_multi_poly::der_press_ent_p(double ent, const Param* ) const {
991 
992  int i = 0 ; // "i" corresponds to the number of domain
993  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
994  // The buffer zone is included in the next zone.
995  for (int l=0; l<npeos-1; l++) {
996  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
997  i++ ;
998  }
999  }
1000 
1001  if (i == 0) {
1002 
1003  if ( ent > double(0) ) {
1004 
1005  if ( ent < 1.e-13 ) {
1006 
1007  return gamma[0]
1008  * ( double(1) + ent/double(2) + ent*ent/double(12) )
1009  / (gamma[0] - double(1)) ;
1010 
1011  }
1012  else {
1013 
1014  return gamma[0] * ent / (double(1) - exp(-ent))
1015  / (gamma[0] - double(1)) ; // mu0[0]/m0=1
1016 
1017  }
1018 
1019  }
1020  else {
1021 
1022  return gamma[0] / (gamma[0] - double(1)) ;
1023 
1024  }
1025 
1026  }
1027  else {
1028 
1029  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
1030  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
1031 
1032  if ( ent < entLarge ) {
1033 
1034  double log10H = log10( ent ) ;
1035  double log10HSmall = log10( entSmall ) ;
1036  double log10HLarge = log10( entLarge ) ;
1037  double dH = log10HLarge - log10HSmall ;
1038  double uu = (log10H - log10HSmall) / dH ;
1039 
1040  double mgam1skapgamSmall = m0*(gamma[i-1]-double(1))
1041  /kappa[i-1]/gamma[i-1] ;
1042  double mgam1skapgamLarge = m0*(gamma[i]-double(1))
1043  /kappa[i]/gamma[i] ;
1044 
1045  double nnSmall = pow( mgam1skapgamSmall
1046  *(exp(entSmall)-mu0[i-1]/m0),
1047  double(1)/(gamma[i-1]-double(1)) ) ;
1048  double nnLarge = pow( mgam1skapgamLarge
1049  *(exp(entLarge)-mu0[i]/m0),
1050  double(1)/(gamma[i]-double(1)) ) ;
1051 
1052  double ppSmall = kappa[i-1] * pow( nnSmall, gamma[i-1] ) ;
1053  double ppLarge = kappa[i] * pow( nnLarge, gamma[i] ) ;
1054 
1055  double eeSmall = mu0[i-1] * nnSmall
1056  + ppSmall / (gamma[i-1] - double(1)) ;
1057  double eeLarge = mu0[i] * nnLarge
1058  + ppLarge / (gamma[i] - double(1)) ;
1059 
1060  double log10PSmall = log10( ppSmall ) ;
1061  double log10PLarge = log10( ppLarge ) ;
1062 
1063  double dlpsdlhSmall = entSmall
1064  * (double(1) + eeSmall / ppSmall) ;
1065  double dlpsdlhLarge = entLarge
1066  * (double(1) + eeLarge / ppLarge) ;
1067 
1068  double dlpsdlhInterpolate = dlpsdlh(log10PSmall, log10PLarge,
1069  dlpsdlhSmall, dlpsdlhLarge,
1070  dH, uu) ;
1071  return dlpsdlhInterpolate ;
1072 
1073  }
1074  else {
1075 
1076  return gamma[i] * ent / (double(1) - (mu0[i]/m0) * exp(-ent))
1077  / (gamma[i] - double(1)) ;
1078 
1079  }
1080 
1081  }
1082 }
1083 
1084 
1085 // dln(p)/dln(n) from enthalpy
1086 //----------------------------
1087 
1088 double Eos_multi_poly::der_press_nbar_p(double ent, const Param* ) const {
1089 
1090  int i = 0 ; // "i" corresponds to the number of domain
1091  // i=0: gamma[0], kappa[0], i=1: gamma[1], kappa[1], .....
1092  // The buffer zone is included in the next zone.
1093  for (int l=0; l<npeos-1; l++) {
1094  if ( ent > entCrit[l]*(double(1)-decInc[l]) ) {
1095  i++ ;
1096  }
1097  }
1098 
1099  if (i == 0) {
1100 
1101  return gamma[0] ;
1102 
1103  }
1104  else {
1105 
1106  double entSmall = entCrit[i-1]*(double(1)-decInc[i-1]) ;
1107  double entLarge = entCrit[i-1]*(double(1)+decInc[i-1]) ;
1108 
1109  if ( ent < entLarge ) {
1110 
1111  double log10H = log10( ent ) ;
1112  double log10HSmall = log10( entSmall ) ;
1113  double log10HLarge = log10( entLarge ) ;
1114 
1115  double dlpsdlnbInterpolate = dlpsdlnb(log10HSmall, log10HLarge,
1116  gamma[i-1], gamma[i],
1117  log10H) ;
1118 
1119  return dlpsdlnbInterpolate ;
1120 
1121  }
1122  else {
1123 
1124  return gamma[i] ;
1125 
1126  }
1127 
1128  }
1129 }
1130 
1131 
1132 //***************************************************//
1133 // Functions which appear in the computation //
1134 //***************************************************//
1135 
1136 double logp(double log10PressSmall, double log10PressLarge,
1137  double dlpsdlhSmall, double dlpsdlhLarge,
1138  double dx, double u) {
1139 
1140  double resu = log10PressSmall * (double(2) * u * u * u
1141  - double(3) * u * u + double(1))
1142  + log10PressLarge * (double(3) * u * u - double(2) * u * u * u)
1143  + dlpsdlhSmall * dx * (u * u * u - double(2) * u * u + u)
1144  - dlpsdlhLarge * dx * (u * u - u * u * u) ;
1145 
1146  return resu ;
1147 
1148 }
1149 
1150 double dlpsdlh(double log10PressSmall, double log10PressLarge,
1151  double dlpsdlhSmall, double dlpsdlhLarge,
1152  double dx, double u) {
1153 
1154  double resu = double(6) * (log10PressSmall - log10PressLarge)
1155  * (u * u - u) / dx
1156  + dlpsdlhSmall * (double(3) * u * u - double(4) * u + double(1))
1157  + dlpsdlhLarge * (double(3) * u * u - double(2) * u) ;
1158 
1159  return resu ;
1160 
1161 }
1162 
1163 double dlpsdlnb(double log10HSmall, double log10HLarge,
1164  double dlpsdlnbSmall, double dlpsdlnbLarge,
1165  double log10H) {
1166 
1167  double resu = log10H * (dlpsdlnbSmall - dlpsdlnbLarge)
1168  / (log10HSmall - log10HLarge)
1169  + (log10HSmall * dlpsdlnbLarge - log10HLarge * dlpsdlnbSmall)
1170  / (log10HSmall - log10HLarge) ;
1171 
1172  return resu ;
1173 
1174 }
1175 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
virtual ~Eos_multi_poly()
Destructor.
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
Lorene prototypes.
Definition: app_hor.h:64
virtual bool operator==(const Eos &) const
Read/write kappa.
double * kappa
Array (size: npeos) of pressure coefficient [unit: ], where and .
Standard units of space, time and mass.
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 double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
const int & get_npeos() const
Returns the number of polytropes npeos.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
void operator=(const Eos_multi_poly &)
Assignment to another Eos_multi_poly.
Eos_multi_poly(int npoly, double *gamma_i, double kappa0_i, double logP1_i, double *logRho_i, double *decInc_i)
Standard constructor (sets m0 to 1).
virtual void sauve(FILE *) const
Save in a file.
Definition: eos.C:179
double kappa0
Pressure coefficient for the crust [unit: ].
double * gamma
Array (size: npeos) of adiabatic index .
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
double * entCrit
Array (size npeos - 1) of the critical enthalpy at which the polytropic EOS changes its index and con...
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 logP1
Exponent of the pressure at the fiducial density .
double * nbCrit
Array (size npeos - 1) of the number density at which the polytropic EOS changes its index and consta...
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
const double & get_kappa(int n) const
Returns the pressure coefficient [unit: ], where and .
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
double * mu0
Array (size: npeos) of the relativistic chemical potential at zero pressure [unit: ...
Base class for a multiple polytropic equation of state.
virtual ostream & operator>>(ostream &) const
Operator >>
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
double * decInc
Array (size npeos - 1) of the percentage which detemines the terminating enthalpy for lower density a...
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
int npeos
Number of polytropic equations of state.
void set_auxiliary()
Computes the auxiliary quantities.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
double m0
Individual particule mass [unit: ].
virtual void sauve(FILE *) const
Save in a file.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
double * logRho
Array (size: npeos - 1) of the exponent of fiducial densities.
const double & get_gamma(int n) const
Returns the adiabatic index .
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.