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
29char 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
82namespace Lorene {
83double logp(double, double, double, double, double, double) ;
84double dlpsdlh(double, double, double, double, double, double) ;
85double dlpsdlnb(double, double, double, double, double) ;
86
87//*********************************************************************
88
89 //--------------------------------------//
90 // Constructors //
91 //--------------------------------------//
92
93// Standard constructor
94Eos_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
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
208Eos_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
341bool 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
385bool Eos_multi_poly::operator!=(const Eos& eos_i) const {
386
387 return !(operator==(eos_i)) ;
388
389}
390
391 //--------------------------//
392 // Outputs //
393 //--------------------------//
394
395void 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
419ostream& 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
493double 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
591double 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
700double 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
798double 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
857double 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
990double 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
1088double 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
1136double 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
1150double 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
1163double 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}
Base class for a multiple polytropic equation of state.
virtual void sauve(FILE *) const
Save in a file.
const double & get_gamma(int n) const
Returns the adiabatic index .
int npeos
Number of polytropic equations of state.
double * decInc
Array (size npeos - 1) of the percentage which detemines the terminating enthalpy for lower density a...
double m0
Individual particule mass [unit: ].
void operator=(const Eos_multi_poly &)
Assignment to another Eos_multi_poly.
double * nbCrit
Array (size npeos - 1) of the number density at which the polytropic EOS changes its index and consta...
const double & get_kappa(int n) const
Returns the pressure coefficient [unit: ], where and .
virtual double der_press_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
double * mu0
Array (size: npeos) of the relativistic chemical potential at zero pressure [unit: ,...
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
double * gamma
Array (size: npeos) of adiabatic index .
const int & get_npeos() const
Returns the number of polytropes npeos.
double * logRho
Array (size: npeos - 1) of the exponent of fiducial densities.
virtual double der_ener_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
double logP1
Exponent of the pressure at the fiducial density .
virtual bool operator!=(const Eos &) const
Comparison operator (difference)
void set_auxiliary()
Computes the auxiliary quantities.
virtual double der_press_nbar_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
virtual ostream & operator>>(ostream &) const
Operator >>
double * kappa
Array (size: npeos) of pressure coefficient [unit: ], where and .
double kappa0
Pressure coefficient for the crust [unit: ].
virtual ~Eos_multi_poly()
Destructor.
double * entCrit
Array (size npeos - 1) of the critical enthalpy at which the polytropic EOS changes its index and con...
virtual double der_nbar_ent_p(double ent, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy.
virtual bool operator==(const Eos &) const
Read/write kappa.
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).
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
Parameter storage.
Definition param.h:125
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
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.