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