LORENE
eos_bf_tabul.C
1/*
2 * Methods of class Eos_bf_tabul.
3 *
4 * (see file eos_bifluid.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2015 Aurelien Sourie
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
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/*
33 * $Id: eos_bf_tabul.C,v 1.5 2017/10/06 12:36:33 a_sourie Exp $
34 * $Log: eos_bf_tabul.C,v $
35 * Revision 1.5 2017/10/06 12:36:33 a_sourie
36 * Cleaning of tabulated 2-fluid EoS class + superfluid rotating star model.
37 *
38 * Revision 1.4 2016/12/05 16:17:51 j_novak
39 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40 *
41 * Revision 1.3 2015/06/11 14:41:59 a_sourie
42 * Corrected minor bug
43 *
44 * Revision 1.2 2015/06/11 13:50:19 j_novak
45 * Minor corrections
46 *
47 * Revision 1.1 2015/06/10 14:39:17 a_sourie
48 * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
49 *
50 *
51 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_tabul.C,v 1.5 2017/10/06 12:36:33 a_sourie Exp $
52 *
53 */
54
55// headers C
56#include <cstdlib>
57#include <cstring>
58#include <cmath>
59
60// Headers Lorene
61#include "eos_bifluid.h"
62#include "param.h"
63#include "eos.h"
64#include "tbl.h"
65#include "utilitaires.h"
66#include "unites.h"
67#include "cmp.h"
68#include "nbr_spx.h"
69
70namespace Lorene {
71
72
73void interpol_herm(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
74 double x, int& i, double& y, double& dy) ;
75
76
77 //----------------------------//
78 // Constructors //
79 //----------------------------//
80
81
82// Standard constructor
83// --------------------
84Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* table,
85 const char* path, double mass1, double mass2) :
86Eos_bifluid(name_i, mass1, mass2)
87{
88 tablename = path ;
89 tablename += "/" ;
90 tablename += table ;
91
92 read_table() ;
93}
94
95// Standard constructor with full filename
96// ---------------------------------------
97Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* file_name,
98 double mass1, double mass2) :
99Eos_bifluid(name_i, mass1, mass2)
100{
101 tablename = file_name ;
102
103 read_table() ;
104}
105
106// Constructor from binary file
107// ----------------------------
108Eos_bf_tabul::Eos_bf_tabul(FILE* fich) :
109Eos_bifluid(fich)
110{
111 char tmp_string[160] ;
112 fread(tmp_string, sizeof(char), 160, fich) ;
113 tablename = tmp_string ;
114
115 read_table() ;
116}
117
118// Constructor from a formatted file
119// ---------------------------------
120Eos_bf_tabul::Eos_bf_tabul(ifstream& fich, const char* table) :
121Eos_bifluid(fich)
122{
123 fich >> tablename ;
124 tablename += "/" ;
125 tablename += table ;
126
127 read_table() ;
128}
129
130Eos_bf_tabul::Eos_bf_tabul(ifstream& fich) :
131Eos_bifluid(fich)
132{
133 fich.ignore(1000, '\n') ;
134 fich >> tablename ; // directory in which the EoSs are stored
135 tablename += name ; // the name of the EoS is set thanks to the call to Eos_bifluid(fich)
136 tablename += "_2f.d" ; // table for the two-fluid EoS
137 read_table() ;
138
139}
140
141
142 //--------------//
143 // Destructor //
144 //--------------//
145
146Eos_bf_tabul::~Eos_bf_tabul(){
147
148 delete mu1_tab ;
149 delete mu2_tab ;
150 delete delta_car_tab ;
151 delete press_tab ;
152 delete n1_tab;
153 delete n2_tab ;
154 delete d2psdmu1dmu2_tab ;
155 delete dpsddelta_car_tab ;
156 delete dn1sddelta_car_tab ;
157 delete dn2sddelta_car_tab ;
158 delete delta_car_n0 ;
159 delete mu1_n0 ;
160 delete mu2_n0 ;
161 delete delta_car_p0 ;
162 delete mu1_p0 ;
163 delete mu2_p0 ;
164 delete mu1_N ;
165 delete n_n_N ;
166 delete press_N ;
167 delete mu2_P ;
168 delete n_p_P ;
169 delete press_P ;
170
171}
172
173 //--------------//
174 // Assignment //
175 //--------------//
176
177void Eos_bf_tabul::operator=(const Eos_bf_tabul& eosi) {
178
179 // Assignment of quantities common to all the derived classes of Eos_bifluid
180 Eos_bifluid::operator=(eosi) ;
181
182 tablename = eosi.tablename;
183 mu1_min = eosi.mu1_min ;
184 mu1_max = eosi.mu1_max ;
185 mu2_min = eosi.mu2_min ;
186 mu2_max = eosi.mu2_max ;
187 delta_car_min = eosi.delta_car_min ;
188 delta_car_max = eosi.delta_car_max ;
189 mu1_tab = eosi.mu1_tab ;
190 mu2_tab = eosi.mu2_tab ;
191 delta_car_tab = eosi.delta_car_tab ;
192 press_tab = eosi.press_tab ;
193 n1_tab = eosi.n1_tab;
194 n2_tab = eosi.n2_tab ;
195 d2psdmu1dmu2_tab = eosi.d2psdmu1dmu2_tab ;
196 dpsddelta_car_tab = eosi.dpsddelta_car_tab ;
197 dn1sddelta_car_tab = eosi.dn1sddelta_car_tab ;
198 dn2sddelta_car_tab = eosi.dn2sddelta_car_tab;
199 delta_car_n0 = eosi.delta_car_n0 ;
200 mu1_n0 = eosi.mu1_n0 ;
201 mu2_n0 = eosi.mu2_n0 ;
202 delta_car_p0 = eosi.delta_car_p0 ;
203 mu1_p0 = eosi.mu1_p0;
204 mu2_p0 = eosi.mu2_p0 ;
205 mu1_N = eosi.mu1_N ;
206 n_n_N = eosi.n_n_N;
207 press_N = eosi.press_N;
208 mu2_P = eosi.mu2_P ;
209 n_p_P = eosi.n_p_P;
210 press_P = eosi.press_P;
211
212}
213
214 //------------//
215 // Outputs //
216 //------------//
217
218void Eos_bf_tabul::sauve(FILE* fich) const {
219
220 Eos_bifluid::sauve(fich) ;
221
222 char tmp_string[160] ;
223 strcpy(tmp_string, tablename.c_str()) ;
224 fwrite(tmp_string, sizeof(char), 160, fich) ;
225
226}
227
228
229ostream& Eos_bf_tabul::operator>>(ostream & ost) const {
230
231 ost <<
232 "EOS of class Eos_bf_tabul : tabulated EOS depending on three variables (relative velocity and chemical potentials of neutrons and protons)."
233 << '\n' ;
234 ost << "Table read from " << tablename << endl ;
235
236 return ost ;
237}
238
239 //------------------------//
240 // Comparison operators //
241 //------------------------//
242
243
244bool Eos_bf_tabul::operator==(const Eos_bifluid& eos_i) const {
245
246 bool resu = true ;
247
248 if ( eos_i.identify() != identify() ) {
249 cout << "The second EOS is not of type Eos_bf_tabul !" << endl ;
250 resu = false ;
251 }
252 else {
253 const Eos_bf_tabul& eos = dynamic_cast<const Eos_bf_tabul&>( eos_i ) ;
254
255 if (eos.tablename != tablename){
256 cout <<
257 "The two Eos_bf_tabul have different names and not refer to the same tables! " << endl ;
258 resu = false ;
259 }
260 }
261
262 return resu ;
263
264}
265
266
267bool Eos_bf_tabul::operator!=(const Eos_bifluid& eos_i) const {
268
269 return !(operator==(eos_i)) ;
270
271}
272
273
274 //------------------------//
275 // Reading of the tables //
276 //------------------------//
277
278
279void Eos_bf_tabul::read_table() {
280
281 using namespace Unites ;
282 double m_b_si = rhonuc_si / 1e44;
283
284 //------------------------------------
285 //---------- 2 fluid table -----------
286 //------------------------------------
287
288 ifstream fich(tablename.data()) ;
289
290 if (!fich) {
291 cerr << "Eos_bf_tabul::read_table(): " << endl ;
292 cerr << "Problem in opening the EOS file!" << endl ;
293 cerr << "While trying to open " << tablename << endl ;
294 cerr << "Aborting..." << endl ;
295 abort() ;
296 }
297
298 fich.ignore(1000, '\n') ; // Jump over the first header
299 fich.ignore(2) ;
300 getline(fich, authors, '\n') ;
301
302 for (int i=0; i<5; i++) { // Jump over the file header
303 fich.ignore(1000, '\n') ; //
304 } //
305
306 int nbp ;
307 fich >> nbp ; fich.ignore(1000, '\n') ; // number of data
308
309 int n_delta, n_mu1, n_mu2; // number of values in delta_car, mu_n and mu_p
310 fich >> n_delta; fich.ignore(1000, '\n') ;
311 fich >> n_mu1; fich.ignore(1000, '\n') ;
312 fich >> n_mu2; fich.ignore(1000, '\n') ;
313
314 if (nbp<=0) {
315 cerr << "Eos_bf_tabul::read_table(): " << endl ;
316 cerr << "Wrong value for the number of lines!" << endl ;
317 cerr << "nbp = " << nbp << endl ;
318 cerr << "Aborting..." << endl ;
319 abort() ;
320 }
321 if (n_delta<=0) {
322 cerr << "Eos_bf_tabul::read_table(): " << endl ;
323 cerr << "Wrong value for the number of values of delta_car!" << endl ;
324 cerr << "n_delta = " << n_delta << endl ;
325 cerr << "Aborting..." << endl ;
326 abort() ;
327 }
328 if (n_mu1<=0) {
329 cerr << "Eos_bf_tabul::read_table(): " << endl ;
330 cerr << "Wrong value for the number of values of mu_n!" << endl ;
331 cerr << "n_mu1 = " << n_mu1 << endl ;
332 cerr << "Aborting..." << endl ;
333 abort() ;
334 }
335 if (n_mu2<=0) {
336 cerr << "Eos_bf_tabul::read_table(): " << endl ;
337 cerr << "Wrong value for the number of values of mu_p!" << endl ;
338 cerr << "n_mu2 = " << n_mu2 << endl ;
339 cerr << "Aborting..." << endl ;
340 abort() ;
341 }
342 if (n_mu2*n_mu1*n_delta != nbp ) {
343 cerr << "Eos_bf_tabul::read_table(): " << endl ;
344 cerr << "Wrong value for the number of lines!" << endl ;
345 cerr << "Aborting..." << endl ;
346 abort() ;
347 }
348
349 for (int i=0; i<3; i++) { // jump over the table
350 fich.ignore(1000, '\n') ; // description
351 }
352
353 mu1_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
354 mu2_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
355 delta_car_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
356 press_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
357 n1_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
358 n2_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
359 d2psdmu1dmu2_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
360 dpsddelta_car_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
361 dn1sddelta_car_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
362 dn2sddelta_car_tab = new Tbl(n_delta, n_mu1, n_mu2) ;
363
364 mu1_tab->set_etat_qcq() ;
365 mu2_tab->set_etat_qcq() ;
366 delta_car_tab->set_etat_qcq() ;
367 press_tab->set_etat_qcq() ;
368 n1_tab->set_etat_qcq() ;
369 n2_tab->set_etat_qcq() ;
370 d2psdmu1dmu2_tab ->set_etat_qcq() ;
371 dpsddelta_car_tab->set_etat_qcq() ;
372 dn1sddelta_car_tab->set_etat_qcq() ;
373 dn2sddelta_car_tab->set_etat_qcq() ;
374
375 //------------------------------------------------------
376 // We have to choose the right unites (SI, cgs , LORENE)
377 //------------------------------------------------------
378
379 // Quantities given by the tabulated EoS (be careful to the units!)
380 double mu1_MeV, mu2_MeV, n1_fm3, n2_fm3;
381 double Knp_Mev_2, press_MeV_fm3;
382 double d2press_dmu1dmu2_MeV_fm3, dn1_ddelta_car_fm3, dn2_ddelta_car_fm3;
383
384 // Stored quantities (in Lorene units)
385 double delta_car_adim, mu1_adim, mu2_adim, n1_01fm3, n2_01fm3, Knp_adim ;
386 double press_adim, dpress_ddelta_car_adim, dn1_ddelta_car_adim, dn2_ddelta_car_adim ;
387 double d2press_dmu1dmu2_adim;
388
389 double hbarc = 197.33 ; // Mev*fm
390 double hbarc3 = hbarc * hbarc * hbarc ;
391
392 // reading of the table.
393 for (int j=0 ; j < n_delta ; j++) {
394 for (int k=0 ; k < n_mu1 ; k++) {
395 for (int l=0 ; l < n_mu2 ; l++) {
396
397 fich >> delta_car_adim ;
398 fich >> mu1_MeV ;
399 mu1_adim = mu1_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
400 fich >> mu2_MeV ;
401 mu2_adim = mu2_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
402 fich >> n1_fm3 ;
403 n1_01fm3 = 10. * n1_fm3 ;
404 fich >> n2_fm3 ;
405 n2_01fm3 = 10. * n2_fm3 ;
406 fich >> Knp_Mev_2 ;
407 Knp_adim = Knp_Mev_2 / ( m_b_si * v_unit * v_unit *10. ) * (mev_si *hbarc3 ) ;
408 dpress_ddelta_car_adim = - Knp_adim * n1_01fm3 * n2_01fm3 * pow( 1.-delta_car_adim, -1.5) / 2. ;
409 fich >> press_MeV_fm3 ;
410 press_adim = press_MeV_fm3 * mevpfm3 ;
411 fich >> d2press_dmu1dmu2_MeV_fm3 ;
412 d2press_dmu1dmu2_adim = d2press_dmu1dmu2_MeV_fm3 * (10. * m_b_si * v_unit * v_unit ) / mev_si ;
413 fich >> dn1_ddelta_car_fm3 ;
414 dn1_ddelta_car_adim = dn1_ddelta_car_fm3 * 10. ;
415 fich >> dn2_ddelta_car_fm3 ;
416 dn2_ddelta_car_adim = dn2_ddelta_car_fm3 * 10. ;
417
418 fich.ignore(1000,'\n') ;
419
420 mu1_tab->set(j,k,l) = mu1_adim ;
421 mu2_tab->set(j,k,l) = mu2_adim ;
422 delta_car_tab->set(j,k,l) = delta_car_adim ;
423 if ((n1_01fm3 <=0) && (n2_01fm3 <=0)) {press_adim = 0. ;}
424 press_tab->set(j,k,l) = press_adim ;
425 n1_tab->set(j,k,l) = n1_01fm3 ;
426 n2_tab->set(j,k,l) = n2_01fm3 ;
427 if ((n1_01fm3 <= 0 ) || (n2_01fm3 <=0)) {
428 d2press_dmu1dmu2_adim = 0. ;
429 dpress_ddelta_car_adim = 0. ;
430 dn1_ddelta_car_adim = 0. ;
431 dn2_ddelta_car_adim = 0. ;
432 }
433 d2psdmu1dmu2_tab ->set(j,k,l) = d2press_dmu1dmu2_adim ;
434 dpsddelta_car_tab->set(j,k,l) = dpress_ddelta_car_adim ;
435 dn1sddelta_car_tab->set(j,k,l) = dn1_ddelta_car_adim ;
436 dn2sddelta_car_tab->set(j,k,l) = dn2_ddelta_car_adim ;
437 }
438
439 fich.ignore(1000, '\n') ;
440
441 }
442
443 fich.ignore(1000, '\n') ;
444
445 }
446
447 delta_car_min = (*delta_car_tab)(0, 0, 0) ;
448 delta_car_max = (*delta_car_tab)(n_delta-1, 0, 0) ;
449 mu1_min = (*mu1_tab)(0, 0, 0) ;
450 mu1_max = (*mu1_tab)(0, n_mu1-1, 0) ;
451 mu2_min = (*mu2_tab)(0, 0, 0) ;
452 mu2_max = (*mu2_tab)(0, 0, n_mu2-1);
453
454 fich.close();
455
456 //---------------------------------------------------------
457 //---------- Table with a single fluid: fluid N -----------
458 //---------------------------------------------------------
459
460 string tablename_1f_n = tablename.c_str() ;
461 tablename_1f_n.replace(tablename_1f_n.end()-5,tablename_1f_n.end(),"_1f_n.d");
462
463 ifstream fichN ;
464 fichN.open(tablename_1f_n.c_str()) ;
465
466 if (!fichN) {
467 cerr << "Eos_bf_tabul::read_table(): " << endl ;
468 cerr << "Problem in opening the EOS file!" << endl ;
469 cerr << "While trying to open " << tablename << endl ;
470 cerr << "Aborting..." << endl ;
471 abort() ;
472 }
473
474 fichN.ignore(1000, '\n') ; // Jump over the first header
475 fichN.ignore(2) ;
476 fichN.ignore(1000, '\n') ;
477
478 for (int i=0; i<5; i++) {
479 fichN.ignore(1000, '\n') ;
480 }
481
482 int nbp_N ; // number of data
483 fichN >> nbp_N ;
484 fichN.ignore(1000, '\n') ;
485 int n_mu1_N; // number of different mu_n
486 fichN >> n_mu1_N ;
487 fichN.ignore(1000, '\n') ;
488
489 if (nbp_N<=0) {
490 cerr << "Eos_bf_tabul::read_table(): " << endl ;
491 cerr << "Wrong value for the number of lines!" << endl ;
492 cerr << "nbp = " << nbp << endl ;
493 cerr << "Aborting..." << endl ;
494 abort() ;
495 }
496 if (n_mu1_N<=0) {
497 cerr << "Eos_bf_tabul::read_table(): " << endl ;
498 cerr << "Wrong value for the number of values of mu_n!" << endl ;
499 cerr << "n_mu1 = " << n_mu1 << endl ;
500 cerr << "Aborting..." << endl ;
501 abort() ;
502 }
503 for (int i=0; i<3; i++) { // jump over the table
504 fichN.ignore(1000, '\n') ; // description
505 }
506
507 mu1_N = new Tbl(n_mu1_N) ;
508 n_n_N = new Tbl(n_mu1_N) ;
509 press_N = new Tbl(n_mu1_N) ;
510
511 mu1_N ->set_etat_qcq() ;
512 n_n_N->set_etat_qcq() ;
513 press_N ->set_etat_qcq() ;
514
515 // Quantities given by the tabulated EoS (be careful to the units!)
516 double mu1_MeV_N, n1_fm3_N, press_MeV_fm3_N;
517
518 // Stored quantities (in Lorene units)
519 double mu1_adimN, n1_01fm3N, press_adimN;
520
521 for (int k=0 ; k < n_mu1_N ; k++) {
522
523 fichN >> mu1_MeV_N ;
524 mu1_adimN = mu1_MeV_N * mev_si / (m_b_si * v_unit * v_unit ) ;
525 fichN >> n1_fm3_N ;
526 n1_01fm3N = 10. * n1_fm3_N ;
527 fichN >> press_MeV_fm3_N;
528 press_adimN = press_MeV_fm3_N * mevpfm3 ;
529 fichN.ignore(1000,'\n') ;
530
531 if ( (n1_01fm3N<0) || (press_adimN < 0)){
532 cout << "Eos_tabul::read_table(): " << endl ;
533 cout << "Negative value in table!" << endl ;
534 cout << "n_neutrons = " << n1_01fm3N << ", p = " << press_adimN << ", "<< endl ;
535 cout << "Aborting..." << endl ;
536 abort() ;
537 }
538
539 mu1_N ->set(k) = mu1_adimN ;
540 n_n_N->set(k) = n1_01fm3N ;
541 press_N ->set(k) = press_adimN ;
542 }
543
544 fichN.close();
545
546 //---------------------------------------------------------
547 //---------- Table with a single fluid: fluid P -----------
548 //---------------------------------------------------------
549
550 string tablename_1f_p = tablename.c_str() ;
551 tablename_1f_p.replace(tablename_1f_p.end()-5,tablename_1f_p.end(),"_1f_p.d");
552
553 ifstream fichP ;
554 fichP.open(tablename_1f_p.c_str()) ;
555
556 if (!fichP) {
557 cerr << "Eos_bf_tabul::read_table(): " << endl ;
558 cerr << "Problem in opening the EOS file!" << endl ;
559 cerr << "While trying to open " << tablename << endl ;
560 cerr << "Aborting..." << endl ;
561 abort() ;
562 }
563
564 fichP.ignore(1000, '\n') ; // Jump over the first header
565 fichP.ignore(2) ;
566 fichP.ignore(1000, '\n') ;
567
568 for (int i=0; i<5; i++) { // jump over the file
569 fichP.ignore(1000, '\n') ; // header
570 }
571
572 int nbp_P ; // number of data
573 fichP >> nbp_P ;
574 fichP.ignore(1000, '\n') ;
575 int n_mu2_P; // number of different values in mu_2
576 fichP >> n_mu2_P ;
577 fichP.ignore(1000, '\n') ;
578
579 if (nbp_P<=0) {
580 cerr << "Eos_bf_tabul::read_table(): " << endl ;
581 cerr << "Wrong value for the number of lines!" << endl ;
582 cerr << "nbp = " << nbp << endl ;
583 cerr << "Aborting..." << endl ;
584 abort() ;
585 }
586 if (n_mu2_P<=0) {
587 cerr << "Eos_bf_tabul::read_table(): " << endl ;
588 cerr << "Wrong value for the number of values of mu_p!" << endl ;
589 cerr << "n_mu2 = " << n_mu2 << endl ;
590 cerr << "Aborting..." << endl ;
591 abort() ;
592 }
593
594 for (int i=0; i<3; i++) { // jump over the table
595 fichP.ignore(1000, '\n') ; // description
596 }
597
598 mu2_P = new Tbl(n_mu2_P) ;
599 n_p_P = new Tbl(n_mu2_P) ;
600 press_P = new Tbl(n_mu2_P) ;
601
602 mu2_P ->set_etat_qcq() ;
603 n_p_P->set_etat_qcq() ;
604 press_P ->set_etat_qcq() ;
605
606 // Quantities given by the tabulated EoS (be careful to the units!)
607 double mu2_MeV_P, n2_fm3_P, press_MeV_fm3_P;
608
609 // Stored quantities (in Lorene units)
610 double mu2_adimP, n2_01fm3P, press_adimP;
611
612 for (int l=0 ; l < n_mu2_P ; l++) {
613
614 fichP >> mu2_MeV_P ;
615 mu2_adimP = mu2_MeV_P * mev_si / (m_b_si * v_unit * v_unit ) ;
616 fichP >> n2_fm3_P ;
617 n2_01fm3P = 10. * n2_fm3_P ;
618 fichP >> press_MeV_fm3_P;
619 press_adimP = press_MeV_fm3_P * mevpfm3 ;
620 fichP.ignore(1000,'\n') ;
621
622 if ( (n2_01fm3P<0) || (press_adimP < 0)){
623 cout << "Eos_tabul::read_table(): " << endl ;
624 cout << "Pegative value in table!" << endl ;
625 cout <<", n_protons " << n2_01fm3P << ", p = " << press_adimP << ", "<< endl ;
626 cout << "Aborting..." << endl ;
627 abort() ;
628 }
629
630 mu2_P ->set(l) = mu2_adimP ;
631 n_p_P->set(l) = n2_01fm3P ;
632 press_P ->set(l) = press_adimP ;
633
634 }
635
636 fichP.close();
637
638 //----------------------------------------------------
639 //---------- Curve corresponding to np = 0 -----------
640 //----------------------------------------------------
641 // Rk: the table is sorted with increasing values of mu_n
642
643 string tablename_np_0 = tablename.c_str() ;
644 tablename_np_0.replace(tablename_np_0.end()-5,tablename_np_0.end(),"_np=0.d");
645
646 ifstream fich1 ;
647 fich1.open(tablename_np_0.c_str()) ;
648
649 if (!fich1) {
650 cerr << "Eos_bf_tabul::read_table(): " << endl ;
651 cerr << "Problem in opening the EOS file!" << endl ;
652 cerr << "While trying to open " << tablename << endl ;
653 cerr << "Aborting..." << endl ;
654 abort() ;
655 }
656
657 int n_delta_n0, n_mu1_n0;
658 fich1 >> n_delta_n0;fich1.ignore(1000, '\n') ;
659 fich1 >> n_mu1_n0;fich1.ignore(1000, '\n') ;
660 fich1.ignore(1000, '\n') ; // Jump over the first header
661
662 delta_car_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
663 mu1_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
664 mu2_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
665
666 delta_car_n0 ->set_etat_qcq() ;
667 mu1_n0->set_etat_qcq() ;
668 mu2_n0->set_etat_qcq() ;
669
670 double delta_car_nn0, mu1_MeV_nn0, mu2_MeV_nn0;
671
672 for (int o = 0; o < n_delta_n0 ; o++ ) {
673
674 for (int p = 0 ; p < n_mu1_n0 ; p++) {
675
676 fich1 >> delta_car_nn0 ;
677 fich1 >> mu1_MeV_nn0 ;
678 fich1 >> mu2_MeV_nn0 ;
679 fich1.ignore(1000,'\n') ;
680
681 delta_car_n0 ->set(o,p) = delta_car_nn0;
682 mu1_n0 ->set(o,p) = mu1_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
683 mu2_n0 ->set(o,p) = mu2_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
684
685 }
686 fich1.ignore(1000,'\n') ;
687
688 }
689
690 fich1.close();
691
692 //----------------------------------------------------
693 //---------- Curve corresponding to nn = 0 -----------
694 //----------------------------------------------------
695 // Rk: the table is sorted with increasing values of mu_p
696
697 string tablename_nn_0 = tablename.c_str() ;
698 tablename_nn_0.replace(tablename_nn_0.end()-5,tablename_nn_0.end(),"_nn=0.d");
699
700 ifstream fich2 ;
701 fich2.open(tablename_nn_0.c_str()) ;
702
703 if (!fich2) {
704 cerr << "Eos_bf_tabul::read_table(): " << endl ;
705 cerr << "Problem in opening the EOS file!" << endl ;
706 cerr << "While trying to open " << tablename << endl ;
707 cerr << "Aborting..." << endl ;
708 abort() ;
709 }
710
711 int n_delta_p0, n_mu2_p0;
712 fich2 >> n_delta_p0;fich2.ignore(1000, '\n') ;
713 fich2 >> n_mu2_p0;fich2.ignore(1000, '\n') ;
714 fich2.ignore(1000, '\n') ;
715
716 delta_car_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
717 mu1_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
718 mu2_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
719
720 delta_car_p0 ->set_etat_qcq() ;
721 mu1_p0->set_etat_qcq() ;
722 mu2_p0 ->set_etat_qcq() ;
723
724 double delta_car_np0, mu1_MeV_np0, mu2_MeV_np0;
725
726 for (int o = 0; o < n_delta_p0 ; o++ ) {
727
728 for (int p = 0 ; p < n_mu2_p0 ; p++) {
729
730 fich2 >> delta_car_np0 ;
731 fich2 >> mu1_MeV_np0 ;
732 fich2 >> mu2_MeV_np0 ;
733 fich2.ignore(1000,'\n') ;
734
735 delta_car_p0 ->set(o,p) = delta_car_np0;
736 mu1_p0 ->set(o,p) = mu1_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
737 mu2_p0 ->set(o,p) = mu2_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
738
739 }
740 fich2.ignore(1000,'\n') ;
741
742 }
743
744 fich2.close();
745
746}
747
748
749 //------------------------------//
750 // Computational routines //
751 //------------------------------//
752
753
754// Complete computational routine giving all thermo variables
755//-----------------------------------------------------------
756
757void Eos_bf_tabul::calcule_interpol(const Cmp& ent1, const Cmp& ent2, const Cmp& delta2,
758 Cmp& nbar1, Cmp& nbar2, Cmp& ener, Cmp& press,
759 Cmp& K_nn, Cmp& K_np, Cmp& K_pp, Cmp& alpha_eos,
760 int nzet, int l_min) const {
761
762 assert(ent1.get_etat() != ETATNONDEF) ;
763 assert(ent2.get_etat() != ETATNONDEF) ;
764 assert(delta2.get_etat() != ETATNONDEF) ;
765
766 const Map* mp = ent1.get_mp() ; // Mapping
767 assert(mp == ent2.get_mp()) ;
768 assert(mp == delta2.get_mp()) ;
769 assert(mp == ener.get_mp()) ;
770
771 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
772
773 nbar1.set_etat_zero() ;
774 nbar2.set_etat_zero() ;
775 ener.set_etat_zero() ;
776 press.set_etat_zero() ;
777 K_nn.set_etat_zero() ;
778 K_np.set_etat_zero() ;
779 K_pp.set_etat_zero() ;
780 alpha_eos.set_etat_zero() ;
781
782 return ;
783
784 }
785
786 nbar1.allocate_all() ;
787 nbar2.allocate_all() ;
788 ener.allocate_all() ;
789 press.allocate_all() ;
790 K_nn.allocate_all() ;
791 K_np.allocate_all() ;
792 K_pp.allocate_all() ;
793 alpha_eos.allocate_all() ;
794
795 const Mg3d* mg = mp->get_mg() ; // Multi-grid
796
797 int nz = mg->get_nzone() ; // total number of domains
798
799 // resu is set to zero in the other domains :
800 if (l_min > 0) {
801
802 nbar1.annule(0, l_min-1) ;
803 nbar2.annule(0, l_min-1) ;
804 ener.annule(0, l_min-1) ;
805 press.annule(0, l_min-1) ;
806 K_nn.annule(0, l_min-1) ;
807 K_np.annule(0, l_min-1) ;
808 K_pp.annule(0, l_min-1) ;
809 alpha_eos.annule(0, l_min-1) ;
810
811 }
812
813 if (l_min + nzet < nz) {
814
815 nbar1.annule(l_min + nzet, nz - 1) ;
816 nbar2.annule(l_min + nzet, nz - 1) ;
817 ener.annule(l_min + nzet, nz - 1) ;
818 press.annule(l_min + nzet, nz - 1) ;
819 K_nn.annule(l_min + nzet, nz - 1) ;
820 K_np.annule(l_min + nzet, nz - 1) ;
821 K_pp.annule(l_min + nzet, nz - 1) ;
822 alpha_eos.annule(l_min + nzet, nz - 1) ;
823
824 }
825
826 int np0 = mp->get_mg()->get_np(0) ;
827 int nt0 = mp->get_mg()->get_nt(0) ;
828
829 for (int l=l_min; l<l_min+nzet; l++) {
830 assert(mp->get_mg()->get_np(l) == np0) ;
831 assert(mp->get_mg()->get_nt(l) == nt0) ;
832 }
833
834 for (int k=0; k<np0; k++) {
835 for (int j=0; j<nt0; j++) {
836 for (int l=l_min; l< l_min + nzet; l++) {
837 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
838
839 double xx, xx1, xx2; // xx1 = H1 = ln(mu1/m1)
840 xx1 = ent1(l,k,j,i) ;
841 xx2 = ent2(l,k,j,i) ;
842 xx = delta2(l,k,j,i) ;
843
844 if (xx < delta_car_min) {
845 cout << "Eos_bf_tabul::calcule_tout : delta2 < delta_car_min !" << endl ;
846 abort() ;
847 }
848 if (xx > delta_car_max) {
849 cout << "Eos_bf_tabul::calcule_tout : delta2 > delta_car_max !" << endl ;
850 abort() ;
851 }
852 if (m_1 * exp(xx1) > mu1_max) {
853 cout << "Eos_bf_tabul::calcule_tout : ent1 > mu1_max !" << endl ;
854 abort() ;
855 }
856 if (m_2 *exp(xx2) > mu2_max) {
857 cout << "Eos_bf_tabul::calcule_tout : ent2 > mu2_max !" << endl ;
858 abort() ;
859 }
860
861 double n1 = 0 , n2 = 0, pressure = 0 ;
862 double alpha_int = 0, K11 = 0, K12 = 0, K22 = 0 ;
863
864 double mu1 = m_1 * exp(xx1);
865 double mu2 = m_2 * exp(xx2);
866
867 if ( (exp(xx1) < 1.) && (exp(xx2) < 1.) ) { // no fluds
868 n1 = 0. ;
869 n2 = 0. ;
870 pressure = 0.;
871 alpha_int = 0 ;
872 K11 = 0 ;
873 K12 = 0 ;
874 K22 = 0 ;
875 }
876
877 else { // at least one fluid is present !
878
879 double p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha_interpol ;
880
881 Eos_bf_tabul::interpol_3d_bifluid(xx, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha_interpol) ;
882
883 n1 = dpsdent1_interpol ;
884 n2 = dpsdent2_interpol ;
885 pressure = p_interpol;
886 alpha_int = alpha_interpol ;
887 }
888
889 if (n1 < 0 ) {
890 cout << "Eos_bf_tabul::calcule_tout : nbar1<0 !" << endl ;
891 // abort() ;
892 n1 = 0 ;
893 }
894 if (n2 < 0 ) {
895 cout << "Eos_bf_tabul::calcule_tout : nbar2<0 !" << endl ;
896 // abort() ;
897 n2 = 0 ;
898 }
899 if (pressure < 0 ) {
900 cout << "Eos_bf_tabul::calcule_tout : pressure < 0 !" << endl ;
901 // abort() ;
902 pressure = 0 ;
903 }
904
905 // Knn
906 if (n1 >0.) {
907 K11 = mu1 / n1 - double(2) * alpha_int * ( 1. - xx) / ( n1 * n1 ) ;
908 }
909 // Kpp
910 if (n2 >0.) {
911 K22 = mu2 / n2 - double(2) * alpha_int * ( 1. - xx) / ( n2 * n2 ) ;
912 }
913 // Knp
914 if ((n1 <= 0.) || (n2 <= 0.) ) {
915 K12 = 0. ;
916 alpha_int = 0. ;
917 }
918 else {
919 K12 = double(2) * alpha_int * pow(1.-xx, 1.5)/ ( n1 * n2 );
920 }
921
922 nbar1.set(l,k,j,i) = n1 ;
923 nbar2.set(l,k,j,i) = n2 ;
924 press.set(l,k,j,i) = pressure ;
925 ener.set(l,k,j,i) = mu1 * n1 + mu2 * n2 - pressure ;
926 K_nn.set(l,k,j,i) = K11 ;
927 K_np.set(l,k,j,i) = K12;
928 K_pp.set(l,k,j,i) = K22 ;
929 alpha_eos.set(l,k,j,i) = alpha_int ;
930
931 }
932 }
933 }
934 }
935
936}
937
938
939// Baryon densities from enthalpies
940//---------------------------------
941
942// this function is not called anymore but should be implemented (virtual function)
943bool Eos_bf_tabul::nbar_ent_p(const double ent1, const double ent2,
944 const double delta2, double& nbar1,
945 double& nbar2) const {
946
947 bool one_fluid = false;
948
949 if (delta2 < delta_car_min) {
950 cout << "Eos_bf_tabul::nbar_ent_p : delta2 < delta_car_min !" << endl ;
951 abort() ;
952 }
953 if (delta2 > delta_car_max) {
954 cout << "Eos_bf_tabul::nbar_ent_p : delta2 > delta_car_max !" << endl ;
955 abort() ;
956 }
957 if (m_1 * exp(ent1) > mu1_max) {
958 cout << "Eos_bf_tabul::nbar_ent_p : ent1 > mu1_max !" << endl ;
959 abort() ;
960 }
961 if (m_2 *exp(ent2) > mu2_max) {
962 cout << "Eos_bf_tabul::nbar_ent_p : ent2 > mu2_max !" << endl ;
963 abort() ;
964 }
965
966 if ( (exp(ent1) < 1.) && (exp(ent2) < 1.) ) {
967 nbar1 = 0. ;
968 nbar2 = 0. ;
969 }
970 else {
971
972 double mu1 = m_1 * exp(ent1);
973 double mu2 = m_2 * exp(ent2);
974
975 double p_interpol ;
976 double dpsdent1_interpol ;
977 double dpsdent2_interpol ;
978 double alpha ;
979
980 Eos_bf_tabul::interpol_3d_bifluid(delta2, mu1, mu2,
981 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
982
983 nbar1 = dpsdent1_interpol ;
984 nbar2 = dpsdent2_interpol ;
985
986 }
987
988 if (nbar1 < 0 ) {
989 cout << "Eos_bf_tabul::nbar_ent_p : nbar1<0 !" << endl ;
990 // abort() ;
991 nbar1 = 0 ;
992 }
993 if (nbar2 < 0 ) {
994 cout << "Eos_bf_tabul::nbar_ent_p : nbar2<0 !" << endl ;
995 // abort() ;
996 nbar2 = 0 ;
997 }
998
999 one_fluid = ((nbar1 <= 0.)||(nbar2 <= 0.)) ;
1000
1001 return one_fluid;
1002
1003}
1004
1005// One fluid sub-EOSs
1006//-------------------
1007
1008// this function is not called anymore but should be implemented (virtual function)
1009double Eos_bf_tabul::nbar_ent_p1(const double ent1) const {
1010
1011 c_est_pas_fait("Eos_bf_tabul::nbar_ent_p1" ) ;
1012
1013 return ent1 ;
1014
1015 /*
1016 double pressN_interpol ;
1017 double n_n_N_interpol ;
1018 double mu1 = m_1 * exp(ent1);
1019 int i =0;
1020
1021 if (exp(ent1) < 1. ) {
1022 n_n_N_interpol = 0. ;
1023 }
1024 else {
1025 interpol_herm( *mu1_N, *press_N, *n_n_N,
1026 mu1, i, pressN_interpol , n_n_N_interpol) ;
1027 }
1028 return n_n_N_interpol ;
1029 */
1030
1031}
1032
1033// this function is not called anymore but should be implemented (virtual function)
1034 double Eos_bf_tabul::nbar_ent_p2(const double ent2) const {
1035
1036 c_est_pas_fait("Eos_bf_tabul::nbar_ent_p2" ) ;
1037
1038 return ent2 ;
1039
1040 /*
1041 double pressP_interpol ;
1042 double n_p_P_interpol ;
1043 double mu2 = m_2 * exp(ent2);
1044 int i =0;
1045 if (exp(ent2) < 1. ) {
1046 n_p_P_interpol = 0. ;
1047 }
1048 else {
1049 interpol_herm( *mu2_P, *press_P, *n_p_P,
1050 mu2, i, pressP_interpol , n_p_P_interpol) ;
1051 }
1052 return n_p_P_interpol ;
1053 */
1054}
1055
1056 // Energy density from baryonic densities
1057//-----------------------------------------
1058
1059// this function is not called anymore but should be implemented (virtual function)
1060double Eos_bf_tabul::ener_nbar_p(const double nbar1, const double nbar2,
1061 const double delta2) const{
1062
1063 c_est_pas_fait("Eos_bf_tabul::ener_nbar_p" ) ;
1064
1065 return nbar1 + nbar2 + delta2;
1066
1067 }
1068
1069// Pressure from baryonic densities
1070//----------------------------------
1071
1072// this function is not called anymore but should be implemented (virtual function)
1073double Eos_bf_tabul::press_nbar_p(const double nbar1, const double nbar2,
1074 const double delta2) const {
1075
1076 c_est_pas_fait("Eos_bf_tabul::press_nbar_p" ) ;
1077
1078 return nbar1 + nbar2 + delta2;
1079
1080}
1081
1082
1083// Pressure from baryonic log-enthalpies
1084//--------------------------------------
1085
1086// this function is not called but can be useful if necessary
1087double Eos_bf_tabul::press_ent_p(const double ent1, const double ent2, const double delta2) const {
1088
1089 if (delta2 < delta_car_min) {
1090 cout << "Eos_bf_tabul::press_ent_p : delta2 < delta_car_min !" << endl ;
1091 abort() ;
1092 }
1093 if (delta2 > delta_car_max) {
1094 cout << "Eos_bf_tabul::press_ent_p : delta2 > delta_car_max !" << endl ;
1095 abort() ;
1096 }
1097 if (m_1 * exp(ent1) > mu1_max) {
1098 cout << "Eos_bf_tabul::press_ent_p : ent1 > mu1_max !" << endl ;
1099 abort() ;
1100 }
1101 if (m_2 * exp(ent2) > mu2_max) {
1102 cout << "Eos_bf_tabul::press_ent_p : ent2 > mu2_max !" << endl ;
1103 abort() ;
1104 }
1105
1106 double pressure ;
1107
1108 if ( (exp(ent1) < 1.) && (exp(ent2) < 1.)) {
1109 //abort();
1110 pressure = 0. ;
1111 }
1112 else {
1113
1114 double mu1 = m_1 * exp(ent1);
1115 double mu2 = m_2 * exp(ent2);
1116
1117 double p_interpol ;
1118 double dpsdent1_interpol ;
1119 double dpsdent2_interpol ;
1120 double alpha ;
1121
1122 Eos_bf_tabul::interpol_3d_bifluid(delta2, mu1, mu2,
1123 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1124
1125 pressure = p_interpol;
1126
1127 }
1128
1129 if (pressure < 0 ) {
1130 cout << "Eos_bf_tabul::press_ent_p : pressure < 0 !" << endl ;
1131 // abort() ;
1132 pressure = 0 ;
1133 }
1134 return pressure ;
1135 }
1136
1137
1138// Energy from baryonic log - enthalpies
1139//--------------------------------------
1140
1141// this function is not called but can be useful if necessary
1142 double Eos_bf_tabul::ener_ent_p(const double ent1, const double ent2,
1143 const double delta2) const {
1144 double energy= 0.;
1145
1146 if ( (exp(ent1) < 1.) && ( exp(ent2) < 1.)) {
1147 energy = 0. ;
1148 }
1149 else {
1150 double mu1 = m_1 * exp(ent1);
1151 double mu2 = m_2 * exp(ent2);
1152
1153 double p_interpol ;
1154 double dpsdent1_interpol ;
1155 double dpsdent2_interpol ;
1156 double alpha ;
1157
1158 Eos_bf_tabul::interpol_3d_bifluid( delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1159
1160 energy = mu1 * dpsdent1_interpol + mu2 * dpsdent2_interpol - p_interpol ;
1161 }
1162
1163 if (energy < 0 ) {
1164 cout << "Eos_bf_tabul::ener_ent_p : energy < 0 !" << endl ;
1165 // abort() ;
1166 energy = 0 ;
1167 }
1168 return energy;
1169
1170}
1171
1172
1173// Alpha from baryonic log - enthalpies
1174//---------------------------------------
1175
1176// this function is not called but can be useful if necessary
1177double Eos_bf_tabul::alpha_ent_p(const double ent1, const double ent2,
1178 const double delta2) const {
1179
1180 if (delta2 < delta_car_min) {
1181 cout << "Eos_bf_tabul::alpha_ent_p : delta2 < delta_car_min !"
1182 << endl ;
1183 abort() ;
1184 }
1185 if (delta2 > delta_car_max) {
1186 cout << "Eos_bf_tabul::alpha_ent_p : delta2 > delta_car_max !"
1187 << endl ;
1188 abort() ;
1189 }
1190 if (m_1 * exp(ent1) > mu1_max) {
1191 cout << "Eos_bf_tabul::alpha_ent_p : ent1 > mu1_max !" << endl ;
1192 abort() ;
1193 }
1194 if (m_2 * exp(ent2) > mu2_max) {
1195 cout << "Eos_bf_tabul::alpha_ent_p : ent2 > mu2_max !" << endl ;
1196 abort() ;
1197 }
1198
1199 double alpha;
1200
1201 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ) {
1202 alpha = 0. ;
1203 }
1204 else {
1205 double mu1 = m_1 * exp(ent1);
1206 double mu2 = m_2 * exp(ent2);
1207 double p_interpol ;
1208 double dpsdent1_interpol ;
1209 double dpsdent2_interpol ;
1210
1211 Eos_bf_tabul::interpol_3d_bifluid( delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1212
1213 }
1214
1215 return alpha;
1216}
1217
1218
1219// Derivatives of energy
1220//----------------------
1221
1222// this function is not called but can be useful if necessary
1223double Eos_bf_tabul::get_K11(const double delta2, const double ent1, const double ent2) const {
1224
1225 double xx = 0.; // K_nn
1226 double mu_1 = m_1 * exp(ent1);
1227 double nbar1;
1228 double nbar2;
1229
1230 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1231 xx = 0. ;
1232 }
1233 else {
1234
1235 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1, nbar2) ;
1236
1237 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1238 if (nbar1 >0.) {
1239 xx = mu_1 / nbar1 - double(2) * alpha * ( 1. - delta2) / ( nbar1 * nbar1 ) ;
1240 }
1241 }
1242
1243 return xx;
1244}
1245
1246// this function is not called but can be useful if necessary
1247double Eos_bf_tabul::get_K22( const double delta2, const double ent1, const double ent2) const {
1248
1249 double xx=0.;
1250 double mu_2 = m_2 * exp (ent2) ;
1251 double nbar1;
1252 double nbar2;
1253
1254 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1255 xx = 0. ;
1256 }
1257 else {
1258
1259 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1,nbar2) ;
1260
1261 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1262 if (nbar2 >0.) {
1263 xx = mu_2 / nbar2 - double(2) * alpha * ( 1. - delta2) / ( nbar2 * nbar2 ) ;
1264 }
1265 }
1266
1267 return xx;
1268}
1269
1270double Eos_bf_tabul::get_K12(const double delta2, const double ent1, const double ent2) const {
1271 double xx =0.;
1272 double nbar1;
1273 double nbar2;
1274
1275 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1276 xx = 0. ;
1277 }
1278 else {
1279
1280 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1,nbar2) ;
1281
1282 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1283 if ((nbar1 <= 0.) || (nbar2 <= 0.) ) {
1284 xx = 0. ;
1285 }
1286 else {
1287 xx = double(2) * alpha * pow(1.-delta2, 1.5)/ ( nbar1 * nbar2 );
1288 }
1289 }
1290
1291 return xx;
1292}
1293
1294// Computes the interpolated values of the pressure, the baryon densities and alpha at the point under consideration from tabulated EoSs.
1295// This routine uses the following 3D interpolation scheme : Hermite interpolation in the chemical potentials
1296// and linear interpolation in the relative speed.
1297// --------------------------------------------------------------------------------------------------------------------------------------
1298void Eos_bf_tabul::interpol_3d_bifluid(const double delta2, const double mu1, const double mu2, double& press, double& nbar1, double& nbar2, double& alpha) const
1299{
1300
1301 assert((*mu1_tab).dim == (*delta_car_tab).dim) ;
1302 assert((*mu2_tab).dim == (*delta_car_tab).dim) ;
1303 assert((*press_tab).dim == (*delta_car_tab).dim) ;
1304 assert((*n1_tab).dim == (*delta_car_tab).dim) ;
1305 assert((*n2_tab).dim == (*delta_car_tab).dim) ;
1306 assert((*d2psdmu1dmu2_tab ).dim == (*delta_car_tab).dim) ;
1307
1308 int nbp1, nbp2, nbp3;
1309 nbp1 = (*delta_car_tab).get_dim(2) ; // number of values of \Delta^2 in the table
1310 nbp2 = (*delta_car_tab).get_dim(1) ; // number of values of \mu_n in the table
1311 nbp3 = (*delta_car_tab).get_dim(0) ; // number of values of \mu_p in the table
1312
1313 Tbl* null_tab = new Tbl(nbp1,nbp2,nbp3) ; // Table whose components are all equal to zero
1314 null_tab->set_etat_zero() ;
1315
1316 int i_near = 0 ;
1317 int j_near = 0 ;
1318 int k_near = 0 ;
1319
1320 // looking for the positions of (delta2,mu1,mu2) in the tables
1321 while ( ( (*delta_car_tab)(i_near,0,0) <= delta2 ) && ( ( nbp1-1 ) > i_near ) ) {
1322 i_near++ ;
1323 }
1324 if (i_near != 0) {
1325 i_near -- ;
1326 }
1327 while ( ( (*mu1_tab)(i_near,j_near, 0) <= mu1 ) && ( ( nbp2-1 ) > j_near ) ) {
1328 j_near++ ;
1329 }
1330 if (j_near != 0) {
1331 j_near -- ;
1332 }
1333 while ( ( (*mu2_tab)( i_near, j_near, k_near) <= mu2) && ( ( nbp3-1 ) > k_near ) ) {
1334 k_near++ ;
1335 }
1336 if (k_near != 0) {
1337 k_near-- ;
1338 }
1339 int i1 = i_near + 1 ;
1340 int j1 = j_near + 1 ;
1341 int k1 = k_near + 1 ;
1342
1343 // The location in the table is refined if necessary
1344 if ( ( (*delta_car_tab)( i_near, j_near, k_near) > delta2 ) && (i_near !=0 ) ) {
1345 i_near -= 1;
1346 i1 -= 1;
1347 }
1348 if ( (delta2 > (*delta_car_tab)( i1, j_near, k_near) ) && (i_near != nbp1 ) ) {
1349 i_near += 1;
1350 i1 += 1;
1351 }
1352 if ( ( (*mu1_tab)( i_near, j_near, k_near) > mu1 ) && (j_near !=0 ) ) {
1353 j_near -= 1;
1354 j1 -= 1;
1355 }
1356 if ( ( mu1 > (*mu1_tab)( i1, j1, k_near) ) && ( j_near != nbp2) ) {
1357 j_near += 1;
1358 j1 += 1;
1359 }
1360 if ( ( (*mu2_tab)( i_near, j_near, k_near) > mu2 ) && (k_near !=0 ) ) {
1361 k_near -= 1;
1362 k1 -= 1;
1363 }
1364 if ( ( mu2 > (*mu2_tab)( i1, j_near, k1) ) && ( k_near != nbp3) ) {
1365 k_near += 1;
1366 k1 += 1;
1367 }
1368
1369 // Check of the location
1370 if ( ( (*delta_car_tab)( i_near, j_near, k_near) > delta2 ) || (delta2 > (*delta_car_tab)( i1, j_near, k_near) ) ) {
1371 cout << "bad location of delta2 in *delta_car_tab " << endl ;
1372 cout << (*delta_car_tab)( i_near, j_near, k_near) << " " << delta2 << " " << (*delta_car_tab)( i1, j_near, k_near) << endl;
1373 abort();
1374 }
1375 if ( ( (*mu1_tab)( i_near, j_near, k_near) > mu1 ) || (mu1 > (*mu1_tab)( i1, j1, k_near) ) ) {
1376 cout << "bad location of mu1 in *mu1_tab " << endl ;
1377 cout << (*mu1_tab)( i_near, j_near, k_near) << " " << mu1 << " " << (*mu1_tab)( i1, j1, k_near) << endl;
1378 abort();
1379 }
1380 if ( ( (*mu2_tab)( i_near, j_near, k_near) > mu2 ) || ( mu2 > (*mu2_tab)( i1, j_near, k1) ) ){
1381 cout << "bad location of mu2 in *mu2_tab "<< endl ;
1382 cout << (*mu2_tab)( i_near, j_near, k_near) << " " << mu2 << " " << (*mu2_tab)( i1, j_near, k1) << endl;
1383 abort();
1384 }
1385
1386 // Values in the slice i
1387 double press_i_near = 0. ;
1388 double nbar1_i_near = 0. ;
1389 double nbar2_i_near = 0. ;
1390 double malpha_i_near = 0. ;
1391 // Values in the slice i+1
1392 double press_i1 = 0. ;
1393 double nbar1_i1 = 0. ;
1394 double nbar2_i1 = 0. ;
1395 double malpha_i1 = 0. ;
1396 // -alpha
1397 double malpha = 0. ;
1398
1399 int n_deltaN = (*delta_car_n0).get_dim(1) ;
1400 int n_mu1N = (*delta_car_n0).get_dim(0) ;
1401 int n_deltaP = (*delta_car_p0).get_dim(1) ;
1402 int n_mu2P = (*delta_car_p0).get_dim(0) ;
1403
1404
1405 /*******************************
1406 * 2D interpolation on Slice i *
1407 *******************************/
1408
1409 // Looking for the table to be used (concerning protons only, neutrons only or both fluids).
1410 // -----------------------------------------------------------------------------------------
1411
1412 int Placei = 0 ; // 0 = two fluids, 1 = only neutrons (fluid 1), 2 = only protons (fluid 2), 3 = no fluids
1413
1414 int i_nearN_i = 0;
1415 int j_nearN_i = 0;
1416 int i_nearP_i = 0;
1417 int j_nearP_i = 0;
1418
1419 if ( mu1 > m_1 ) // both fluids are present or only the neutron fluid (fluid 1)
1420 {
1421 // to find if one or two fluid(s) is (are) present, we compare the position under consideration
1422 // with the curve n_p = 0 for the EoS which is used.
1423 // Note that the following procedure is adapted to DDH and DDHdelta EoSs (close to beta eq. and corotation)
1424 // but the curve n_p = 0 could be possibly much more complicated depending on the EoS...
1425 while ( ( (*delta_car_n0)(i_nearN_i,0) <= (*delta_car_tab)(i_near, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i ) ) {
1426 i_nearN_i++ ;
1427 }
1428 if (i_nearN_i != 0) {
1429 i_nearN_i -- ;
1430 }
1431 while ( ( (*mu1_n0)(i_nearN_i,j_nearN_i) <= mu1 ) && ( ( n_mu1N-1 ) > j_nearN_i ) ) {
1432 j_nearN_i++ ;
1433 }
1434 if (j_nearN_i != 0) {
1435 j_nearN_i -- ;
1436 }
1437
1438 // some checks
1439 if ( ( (*delta_car_n0)(i_nearN_i,0) > (*delta_car_tab)(i_near, j_near, k_near) ) || ((*delta_car_tab)(i_near, j_near, k_near) > (*delta_car_n0)(i_nearN_i+1,0) ) )
1440 {
1441 cout << " bad location of delta_car_tab_i in *delta_car_n0 (courbe limite np = 0) " << endl ;
1442 cout << (*delta_car_n0)(i_nearN_i,0) << " " << (*delta_car_tab)(i_near, j_near, k_near) << " " << (*delta_car_n0)(i_nearN_i+1,0) << endl;
1443 abort();
1444 }
1445 if ( ( (*mu1_n0)(i_nearN_i,j_nearN_i) > mu1 ) || (mu1 > (*mu1_n0)(i_nearN_i,j_nearN_i+1) ) )
1446 {
1447 cout << " bad location of mu_n in *mu1_n0 (limit curve np = 0) " << endl ;
1448 cout << (*mu1_n0)(i_nearN_i,j_nearN_i) << " " << mu1 << " " << (*mu1_n0)(i_nearN_i,j_nearN_i+1) << endl;
1449 abort();
1450 }
1451
1452 double aN_i, bN_i;
1453 aN_i = ((*mu2_n0)(i_nearN_i,j_nearN_i+1) - (*mu2_n0)(i_nearN_i,j_nearN_i) ) / ((*mu1_n0)(i_nearN_i,j_nearN_i+1) - (*mu1_n0)(i_nearN_i,j_nearN_i) ) ;
1454 bN_i = (*mu2_n0)(i_nearN_i,j_nearN_i) - aN_i * (*mu1_n0)(i_nearN_i,j_nearN_i) ;
1455 double zN_i = aN_i * mu1 + bN_i ;
1456
1457 if (zN_i < mu2)
1458 {
1459 Placei = 0; // two fluids
1460 }
1461 else
1462 {
1463 Placei = 1 ; // fluid 1 only
1464 }
1465 }
1466
1467 else // both fluids are present or only the proton fluid (fluid 2) or no fluids !
1468 {
1469 if ( mu2 <= m_2)
1470 {
1471 Placei = 3 ; // no fluids
1472 }
1473 else
1474 {
1475
1476 // to find if one or two fluid(s) is (are) present, we compare the position under consideration
1477 // with the curve n_n = 0 for the EoS which is used.
1478 // Note that the following procedure is adapted to DDH and DDHdelta EoSs (close to beta eq. and corotation)
1479 // but the curve n_n = 0 could be possibly much more complicated depending on the EoS...
1480 while ( ( (*delta_car_p0)(i_nearP_i,0) <= (*delta_car_tab)(i_near, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i ) ) {
1481 i_nearP_i++ ;
1482 }
1483 if (i_nearP_i != 0) {
1484 i_nearP_i -- ;
1485 }
1486 while ( ( (*mu2_p0)(i_nearP_i,j_nearP_i) <= mu2 ) && ( ( n_mu2P-1 ) > j_nearP_i ) ) {
1487 j_nearP_i++ ;
1488 }
1489 if (j_nearP_i != 0) {
1490 j_nearP_i -- ;
1491 }
1492
1493 // some checks
1494 if ( ( (*delta_car_p0)(i_nearP_i,0) > (*delta_car_tab)(i_near, j_near, k_near) ) || ((*delta_car_tab)(i_near, j_near, k_near) > (*delta_car_p0)(i_nearP_i+1,0) ) )
1495 {
1496 cout << " bad location of delta_car_tab_i in *delta_car_p0 (courbe limite nn = 0) " << endl ;
1497 cout << (*delta_car_p0)(i_nearP_i,0) << " " << (*delta_car_tab)(i_near, j_near, k_near) << " " << (*delta_car_p0)(i_nearP_i+1,0) << endl;
1498 abort();
1499 }
1500 if ( ( (*mu2_p0)(i_nearP_i,j_nearP_i) > mu2 ) || (mu2 > (*mu2_p0)(i_nearP_i,j_nearP_i+1) ) ) {
1501 cout << " bad location of mu_p in *mu2_p0 (limit curve nn = 0) " << endl ;
1502 cout << (*mu2_p0)(i_nearP_i,j_nearP_i) << " " << mu2 << " " << (*mu2_p0)(i_nearP_i,j_nearP_i+1) << endl;
1503 abort();
1504 }
1505
1506 double aP_i, bP_i;
1507 aP_i = ( (*mu2_p0)(i_nearP_i,j_nearP_i+1) - (*mu2_p0)(i_nearP_i,j_nearP_i) ) / ( (*mu1_p0)(i_nearP_i,j_nearP_i+1) - (*mu1_p0)(i_nearP_i,j_nearP_i) ) ;
1508 bP_i = (*mu2_p0)(i_nearP_i,j_nearP_i) - aP_i * (*mu1_p0)(i_nearP_i,j_nearP_i) ;
1509 double yP_i = (mu2- bP_i) /aP_i ;
1510
1511 if (yP_i < mu1)
1512 {
1513 Placei = 0; // two fluids
1514 }
1515 else
1516 {
1517 Placei = 2 ; // fluid 2 only
1518 }
1519
1520 }
1521
1522 } // end of the search of the table to be used.
1523
1524 /*
1525 cout << mu2* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << mu1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13) << endl ;
1526 cout << " Placei = " << Placei << endl ;
1527 */
1528
1529 // Interpolation in the different areas
1530 if (Placei == 3 ) // no fluids
1531 {
1532 press_i_near = 0. ;
1533 nbar1_i_near = 0. ;
1534 nbar2_i_near = 0. ;
1535 malpha_i_near = 0.;
1536 }
1537 else if (Placei == 1 ) // fluid 1 only
1538 {
1539 malpha_i_near = 0.;
1540 nbar2_i_near = 0. ;
1541 int i = 0;
1542 interpol_herm(*mu1_N, *press_N, *n_n_N, mu1, i, press_i_near , nbar1_i_near ) ;
1543 if (press_i_near < 0.) {
1544 cout << " INTERPOLATION FLUID N --> negative pressure " << endl ;
1545 abort();
1546 }
1547 if (nbar1_i_near < 0.) {
1548 cout << " INTERPOLATION FLUID N --> negative density " << endl ;
1549 abort();
1550 }
1551 }
1552 else if (Placei == 2 ) // fluid 2 only
1553 {
1554 malpha_i_near = 0.;
1555 nbar1_i_near = 0. ;
1556 int i =0;
1557 interpol_herm( *mu2_P, *press_P, *n_p_P, mu2, i, press_i_near, nbar2_i_near) ;
1558 if (press_i_near < 0.) {
1559 cout << " INTERPOLATION FLUID P --> negative pressure " << endl ;
1560 abort();
1561 }
1562 if (nbar2_i_near < 0.) {
1563 cout << " INTERPOLATION FLUID P --> negative density " << endl ;
1564 abort();
1565 }
1566 }
1567 else if (Placei == 0 ) { // two fluids
1568
1569 // interpolation of press, nbar1 and nbar2
1570
1571 Eos_bf_tabul::interpol_2d_Tbl3d(i_near, j_near, k_near, *mu1_tab, *mu2_tab,
1572 *press_tab, *n1_tab, *n2_tab, *d2psdmu1dmu2_tab ,
1573 mu1, mu2, press_i_near, nbar1_i_near , nbar2_i_near) ;
1574
1575 // interpolation of malpha
1576 double der1 = 0., der2 = 0. ;
1577
1578 Eos_bf_tabul::interpol_2d_Tbl3d(i_near, j_near, k_near, *mu1_tab, *mu2_tab,
1579 *dpsddelta_car_tab, *dn1sddelta_car_tab, *dn2sddelta_car_tab, *null_tab,
1580 mu1, mu2, malpha_i_near, der1 , der2) ;
1581
1582 if (press_i_near < 0.) {
1583 //cout << press_i_near << " --> negative pressure " << endl ;
1584 press_i_near = 0 ; // interpolation for very small values of press can lead to press <0...
1585 // to check the smallness of press -> abort();
1586 }
1587 if (nbar1_i_near < 0.) {
1588 //cout << nbar1_i_near << " --> negative density nbar1 " << endl ;
1589 nbar1_i_near = 0 ; // interpolation for very small values of nbar1 can lead to nbar1 <0...
1590 }
1591 if (nbar2_i_near < 0.) {
1592 //cout << nbar2_i_near << " --> negative density nbar2 " << endl ;
1593 nbar2_i_near = 0 ; // interpolation for very small values of nbar2 can lead to nbar2 <0...
1594 }
1595
1596 }
1597
1598
1599 /***********************************
1600 * 2D interpolation on Slice i + 1 *
1601 ***********************************/
1602
1603 // Looking for the table to be used (concerning protons only, neutrons only or both fluids).
1604 // -----------------------------------------------------------------------------------------
1605
1606 int Placei1 = 0 ; // 0 = two fluids, 1 = only neutrons (fluid 1), 2 = only protons (fluid 2), 3 = no fluids
1607
1608 int i_nearN_i1 = 0;
1609 int j_nearN_i1 = 0;
1610 int i_nearP_i1 = 0;
1611 int j_nearP_i1 = 0;
1612
1613 if ( mu1 > m_1 ) // both fluids are present or only the neutron fluid (fluid 1)
1614 {
1615 // to find if one or two fluid(s) is (are) present, we compare the position under consideration
1616 // with the curve n_p = 0 for the EoS which is used.
1617 // Note that the following procedure is adapted to DDH and DDHdelta EoSs (close to beta eq. and corotation)
1618 // but the curve n_p = 0 could be possibly much more complicated depending on the EoS...
1619 while ( ( (*delta_car_n0)(i_nearN_i1,0) <= (*delta_car_tab)(i_near+1, j_near, k_near) ) && ( ( n_deltaN-1 ) > i_nearN_i1 ) ) {
1620 i_nearN_i1++ ;
1621 }
1622 if (i_nearN_i1 != 0) {
1623 i_nearN_i1 -- ;
1624 }
1625 while ( ( (*mu1_n0)(i_nearN_i1,j_nearN_i1) <= mu1 ) && ( ( n_mu1N-1 ) > j_nearN_i1 ) ) {
1626 j_nearN_i1++ ;
1627 }
1628 if (j_nearN_i1 != 0) {
1629 j_nearN_i1 -- ;
1630 }
1631
1632 // some checks
1633 if ( ( (*delta_car_n0)(i_nearN_i1,0) > (*delta_car_tab)(i_near+1, j_near, k_near) ) || ((*delta_car_tab)(i_near+1, j_near, k_near) > (*delta_car_n0)(i_nearN_i1+1,0) ) )
1634 {
1635 cout << " bad location of delta_car_tab_i+1 in *delta_car_n0 (courbe limite np = 0) " << endl ;
1636 cout << (*delta_car_n0)(i_nearN_i1,0) << " " << (*delta_car_tab)(i_near+1, j_near, k_near) << " " << (*delta_car_n0)(i_nearN_i1+1,0) << endl;
1637 abort();
1638 }
1639 if ( ( (*mu1_n0)(i_nearN_i1,j_nearN_i1) > mu1 ) || (mu1 > (*mu1_n0)(i_nearN_i1,j_nearN_i1+1) ) )
1640 {
1641 cout << " bad location of mu_n in *mu1_n0 (limit curve np = 0) " << endl ;
1642 cout << (*mu1_n0)(i_nearN_i1,j_nearN_i1) << " " << mu1 << " " << (*mu1_n0)(i_nearN_i1,j_nearN_i1+1) << endl;
1643 abort();
1644 }
1645
1646 double aN_i1, bN_i1;
1647 aN_i1 = ((*mu2_n0)(i_nearN_i1,j_nearN_i1+1) - (*mu2_n0)(i_nearN_i1,j_nearN_i1) ) / ((*mu1_n0)(i_nearN_i1,j_nearN_i1+1) - (*mu1_n0)(i_nearN_i1,j_nearN_i1) ) ;
1648 bN_i1 = (*mu2_n0)(i_nearN_i1,j_nearN_i1) - aN_i1 * (*mu1_n0)(i_nearN_i1,j_nearN_i1) ;
1649 double zN_i1 = aN_i1 * mu1 + bN_i1 ;
1650
1651 if (zN_i1 < mu2)
1652 {
1653 Placei1 = 0 ; // two fluids
1654 }
1655 else {
1656 Placei1 = 1 ; // fluid 1 only
1657 }
1658
1659 }
1660 else // both fluids are present or only the proton fluid (fluid 2) or no fluids !
1661 {
1662 if ( mu2 <= m_2)
1663 {
1664 Placei = 3 ; // no fluids
1665 }
1666 else
1667 {
1668
1669 // to find if one or two fluid(s) is (are) present, we compare the position under consideration
1670 // with the curve n_n = 0 for the EoS which is used.
1671 // Note that the following procedure is adapted to DDH and DDHdelta EoSs (close to beta eq. and corotation)
1672 // but the curve n_n = 0 could be possibly much more complicated depending on the EoS...
1673 while ( ( (*delta_car_p0)(i_nearP_i1,0) <= (*delta_car_tab)(i_near+1, j_near, k_near)) && ( ( n_deltaP-1 ) > i_nearP_i1) ) {
1674 i_nearP_i1++ ;
1675 }
1676 if (i_nearP_i1 != 0) {
1677 i_nearP_i1 -- ;
1678 }
1679 while ( ( (*mu2_p0)(i_nearP_i1,j_nearP_i1) <= mu2 ) && ( ( n_mu2P-1 ) > j_nearP_i1 ) ) {
1680 j_nearP_i1++ ;
1681 }
1682 if (j_nearP_i1 != 0) {
1683 j_nearP_i1 -- ;
1684 }
1685
1686 // some checks
1687 if ( ( (*delta_car_p0)(i_nearP_i1,0) > (*delta_car_tab)(i_near+1, j_near, k_near) ) || ((*delta_car_tab)(i_near+1, j_near, k_near) > (*delta_car_p0)(i_nearP_i1+1,0) ) )
1688 {
1689 cout << " bad location of delta_car_tab_i+1 in *delta_car_p0 (courbe limite nn = 0) " << endl ;
1690 cout << (*delta_car_p0)(i_nearP_i1,0) << " " << (*delta_car_tab)(i_near+1, j_near, k_near) << " " << (*delta_car_p0)(i_nearP_i1+1,0) << endl;
1691 abort();
1692 }
1693 if ( ( (*mu2_p0)(i_nearP_i1,j_nearP_i1) > mu2 ) || (mu2 > (*mu2_p0)(i_nearP_i1,j_nearP_i1+1) ) ) {
1694 cout << " bad location of mu_p in *mu2_p0 (limit curve nn = 0) " << endl ;
1695 cout << (*mu2_p0)(i_nearP_i1,j_nearP_i1) << " " << mu2 << " " << (*mu2_p0)(i_nearP_i1,j_nearP_i1+1) << endl;
1696 abort();
1697 }
1698
1699 double aP_i1, bP_i1;
1700 aP_i1 = ( (*mu2_p0)(i_nearP_i1,j_nearP_i1+1) - (*mu2_p0)(i_nearP_i1,j_nearP_i1) ) / ( (*mu1_p0)(i_nearP_i1,j_nearP_i1+1) - (*mu1_p0)(i_nearP_i1,j_nearP_i1) ) ;
1701 bP_i1 = (*mu2_p0)(i_nearP_i1,j_nearP_i1) - aP_i1 * (*mu1_p0)(i_nearP_i1,j_nearP_i1) ;
1702 double yP_i1 = (mu2- bP_i1) /aP_i1 ;
1703
1704 if (yP_i1 < mu1)
1705 {
1706 Placei = 0; // two fluids
1707 }
1708 else
1709 {
1710 Placei = 2 ; // fluid 2 only
1711 }
1712
1713 }
1714
1715 } // end of the search of the table to be used.
1716
1717 /*
1718 cout << mu2* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << mu1 *2.99792458E+8 * 2.99792458E+8 * 1.66e-27 / (1.6021892E-13) << endl ;
1719 cout << " Placei = " << Placei << endl ;
1720 */
1721
1722 // Interpolation in the different areas
1723 if (Placei == 3 ) // no fluids
1724 {
1725 press_i1 = 0. ;
1726 nbar1_i1 = 0. ;
1727 nbar2_i1 = 0. ;
1728 malpha_i1 = 0.;
1729 }
1730 else if (Placei1 == 1 ) // fluid 1 only
1731 {
1732 malpha_i1 = 0. ;
1733 nbar2_i1 = 0. ;
1734 int i =0;
1735 interpol_herm(*mu1_N, *press_N, *n_n_N, mu1, i, press_i1 , nbar1_i1 ) ;
1736 if (press_i1 < 0.) {
1737 cout << " INTERPOLATION FLUID N i+1 --> negative pressure " << endl ;
1738 abort();
1739 }
1740 if (nbar1_i1 < 0.) {
1741 cout << " INTERPOLATION FLUID N i+1--> negative density " << endl ;
1742 abort();
1743 }
1744 }
1745 else if (Placei1 == 2 ) // fluid 2 only
1746 {
1747 malpha_i1 = 0.;
1748 nbar1_i1 = 0. ;
1749 int i =0;
1750 interpol_herm( *mu2_P, *press_P, *n_p_P, mu2, i, press_i1, nbar2_i1) ;
1751 if (press_i1 < 0.) {
1752 cout << " INTERPOLATION FLUID P i+1--> negative pressure " << endl ;
1753 abort();
1754 }
1755 if (nbar2_i1 < 0.) {
1756 cout << " INTERPOLATION FLUID P i+1 --> negative density " << endl ;
1757 abort();
1758 }
1759 }
1760 else if (Placei1 == 0 ) { // two fluids
1761
1762 // interpolation of press, nbar1 and nbar2
1763
1764 Eos_bf_tabul::interpol_2d_Tbl3d(i1, j_near, k_near, *mu1_tab, *mu2_tab,
1765 *press_tab, *n1_tab, *n2_tab, *d2psdmu1dmu2_tab ,
1766 mu1, mu2, press_i1, nbar1_i1 , nbar2_i1) ;
1767
1768 // interpolation of malpha
1769 double der1b = 0., der2b = 0.;
1770
1771 Eos_bf_tabul::interpol_2d_Tbl3d(i1, j_near, k_near, *mu1_tab, *mu2_tab,
1772 *dpsddelta_car_tab, *dn1sddelta_car_tab, *dn2sddelta_car_tab, *null_tab,
1773 mu1, mu2, malpha_i1, der1b , der2b ) ;
1774
1775
1776 if (press_i1 < 0.) {
1777 //cout << press_i1 << " --> negative pressure " << endl ;
1778 press_i1 = 0 ; // interpolation for very small values of press can lead to press <0...
1779 // to check the smallness of press -> abort();
1780 }
1781 if (nbar1_i1 < 0.) {
1782 //cout << nbar1_i1 << " --> negative density nbar1 " << endl ;
1783 nbar1_i1 = 0 ; // interpolation for very small values of nbar1 can lead to nbar1 <0...
1784 }
1785 if (nbar2_i1 < 0.) {
1786 //cout << nbar2_i1 << " --> negative density nbar2 " << endl ;
1787 nbar2_i1 = 0 ; // interpolation for very small values of nbar2 can lead to nbar2 <0...
1788 }
1789
1790 }
1791
1792
1793 /***********************************
1794 * linear interpolation in Delta^2 *
1795 ***********************************/
1796
1797 double x1 = (*delta_car_tab)(i_near, 0, 0) ;
1798 double x2 = (*delta_car_tab)(i1, 0, 0) ;
1799 double x12 = x1-x2 ;
1800
1801 //pressure
1802 double y1 = press_i_near;
1803 double y2 = press_i1;
1804 double a = (y1-y2)/x12 ;
1805 double b = (x1*y2-y1*x2)/x12 ;
1806 press = delta2*a+b ;
1807
1808
1809 //nbar1
1810 double y1_y = nbar1_i_near;
1811 double y2_y = nbar1_i1;
1812 double a_y = (y1_y-y2_y)/x12 ;
1813 double b_y = (x1*y2_y-y1_y*x2)/x12 ;
1814 nbar1 = delta2*a_y+b_y ;
1815
1816 //nbar2
1817 double y1_z = nbar2_i_near;
1818 double y2_z = nbar2_i1;
1819 double a_z = (y1_z-y2_z)/x12 ;
1820 double b_z = (x1*y2_z-y1_z*x2)/x12 ;
1821 nbar2 = delta2*a_z+b_z ;
1822
1823 // for alpha
1824 double y1_alpha = malpha_i_near;
1825 double y2_alpha = malpha_i1;
1826 double a_alpha = (y1_alpha-y2_alpha)/x12 ;
1827 double b_alpha = (x1*y2_alpha-y1_alpha*x2)/x12 ;
1828 malpha = delta2*a_alpha+b_alpha ;
1829 alpha = -malpha ;
1830
1831 delete null_tab;
1832 return;
1833
1834}
1835
1836// routine used in interpol_3d_bifluid to perform a 2D thermodynamically consistent interpolation
1837// on each slice of constant delta_car
1838//----------------------------------------------------------------------------------------------
1839void Eos_bf_tabul::interpol_2d_Tbl3d(const int i, const int j, const int k, const Tbl& ytab, const Tbl& ztab,
1840 const Tbl& ftab, const Tbl& dfdytab, const Tbl& dfdztab, const Tbl& d2fdydztab,
1841 const double y, const double z, double& f, double& dfdy, double& dfdz) const
1842{
1843
1844 assert(dfdytab.dim == ftab.dim ) ;
1845 assert(dfdztab.dim == ftab.dim ) ;
1846 assert(d2fdydztab.dim == ftab.dim ) ;
1847
1848 int j1 = j+1 ;
1849 int k1 = k+1 ;
1850
1851 double dy = ytab(i, j1, k) - ytab(i, j, k) ;
1852 double dz = ztab(i, j, k1) - ztab(i, j, k) ;
1853
1854 double u = (y - ytab(i, j, k)) / dy ;
1855 double v = (z - ztab(i, j, k)) / dz ;
1856
1857 double u2 = u*u ; double v2 = v*v ;
1858 double u3 = u2*u ; double v3 = v2*v ;
1859
1860 double psi0_u = 2.*u3 - 3.*u2 + 1. ;
1861 double psi0_1mu = -2.*u3 + 3.*u2 ;
1862 double psi1_u = u3 - 2.*u2 + u ;
1863 double psi1_1mu = -u3 + u2 ;
1864 double psi0_v = 2.*v3 - 3.*v2 + 1. ;
1865 double psi0_1mv = -2.*v3 + 3.*v2 ;
1866 double psi1_v = v3 - 2.*v2 + v ;
1867 double psi1_1mv = -v3 + v2 ;
1868
1869 f = ftab(i, j, k) * psi0_u * psi0_v
1870 + ftab(i, j1, k) * psi0_1mu * psi0_v
1871 + ftab(i, j, k1) * psi0_u * psi0_1mv
1872 + ftab(i, j1, k1) * psi0_1mu * psi0_1mv ;
1873
1874 f += (dfdytab(i, j, k) * psi1_u * psi0_v
1875 - dfdytab(i, j1, k) * psi1_1mu * psi0_v
1876 + dfdytab(i, j, k1) * psi1_u * psi0_1mv
1877 - dfdytab(i, j1, k1) * psi1_1mu * psi0_1mv) * dy ;
1878
1879 f += (dfdztab(i, j, k) * psi0_u * psi1_v
1880 + dfdztab(i, j1, k) * psi0_1mu * psi1_v
1881 - dfdztab(i, j, k1) * psi0_u * psi1_1mv
1882 - dfdztab(i, j1, k1) * psi0_1mu * psi1_1mv) * dz ;
1883
1884 f += (d2fdydztab(i, j, k) * psi1_u * psi1_v
1885 - d2fdydztab(i, j1, k) * psi1_1mu * psi1_v
1886 - d2fdydztab(i, j, k1) * psi1_u * psi1_1mv
1887 + d2fdydztab(i, j1, k1) * psi1_1mu * psi1_1mv) * dy * dz ;
1888
1889 double dpsi0_u = 6.*(u2 - u) ;
1890 double dpsi0_1mu = 6.*(u2 - u) ;
1891 double dpsi1_u = 3.*u2 - 4.*u + 1. ;
1892 double dpsi1_1mu = 3.*u2 - 2.*u ;
1893
1894 dfdy = (ftab(i, j, k) * dpsi0_u * psi0_v
1895 - ftab(i, j1, k) * dpsi0_1mu * psi0_v
1896 + ftab(i, j, k1) * dpsi0_u * psi0_1mv
1897 - ftab(i, j1, k1) * dpsi0_1mu * psi0_1mv ) / dy;
1898
1899 dfdy += (dfdytab(i, j, k) * dpsi1_u * psi0_v
1900 + dfdytab(i, j1, k) * dpsi1_1mu * psi0_v
1901 + dfdytab(i, j, k1) * dpsi1_u * psi0_1mv
1902 + dfdytab(i, j1, k1) * dpsi1_1mu * psi0_1mv) ;
1903
1904 dfdy += (dfdztab(i, j, k) * dpsi0_u* psi1_v
1905 - dfdztab(i, j1, k) * dpsi0_1mu * psi1_v
1906 - dfdztab(i, j, k1) * dpsi0_u * psi1_1mv
1907 + dfdztab(i, j1, k1) * dpsi0_1mu * psi1_1mv) * dz /dy ;
1908
1909 dfdy += (d2fdydztab(i, j, k) * dpsi1_u * psi1_v
1910 + d2fdydztab(i, j1, k) * dpsi1_1mu * psi1_v
1911 - d2fdydztab(i, j, k1) * dpsi1_u * psi1_1mv
1912 - d2fdydztab(i, j1, k1) * dpsi1_1mu * psi1_1mv) * dz ;
1913
1914 double dpsi0_v = 6.*(v2 - v) ;
1915 double dpsi0_1mv = 6.*(v2 - v) ;
1916 double dpsi1_v = 3.*v2 - 4.*v + 1. ;
1917 double dpsi1_1mv = 3.*v2 - 2.*v ;
1918
1919 dfdz = (ftab(i, j, k) * psi0_u * dpsi0_v
1920 + ftab(i, j1, k) * psi0_1mu * dpsi0_v
1921 - ftab(i, j, k1) * psi0_u * dpsi0_1mv
1922 - ftab(i, j1, k1) * psi0_1mu * dpsi0_1mv) / dz ;
1923
1924 dfdz += (dfdytab(i, j, k) * psi1_u * dpsi0_v
1925 - dfdytab(i, j1, k) * psi1_1mu * dpsi0_v
1926 - dfdytab(i, j, k1) * psi1_u * dpsi0_1mv
1927 + dfdytab(i, j1, k1) * psi1_1mu * dpsi0_1mv) * dy / dz ;
1928
1929 dfdz += (dfdztab(i, j, k) * psi0_u * dpsi1_v
1930 + dfdztab(i, j1, k) * psi0_1mu * dpsi1_v
1931 + dfdztab(i, j, k1) * psi0_u * dpsi1_1mv
1932 + dfdztab(i, j1, k1) * psi0_1mu * dpsi1_1mv) ;
1933
1934 dfdz += (d2fdydztab(i, j, k) * psi1_u* dpsi1_v
1935 - d2fdydztab(i, j1, k) * psi1_1mu * dpsi1_v
1936 + d2fdydztab(i, j, k1) * psi1_u* dpsi1_1mv
1937 - d2fdydztab(i, j1, k1) * psi1_1mu * dpsi1_1mv) * dy;
1938
1939 return ;
1940
1941}
1942
1943// Conversion functions
1944// ---------------------
1945
1946//This function is necessary for "Et_rot", which needs an eos of type "Eos"
1947//But this eos is not used in the code, except for the construction of the star
1948Eos* Eos_bf_tabul::trans2Eos() const {
1949
1950 Eos_poly* eos_simple = new Eos_poly(2.,0.016,1.008) ; // we can take whatever we want that makes sense as parameters
1951
1952 return eos_simple ;
1953}
1954
1955}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::allocate_all
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition cmp.C:326
Lorene::Cmp::get_etat
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:351
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Lorene::Cmp::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:292
Lorene::Cmp::get_mp
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Lorene::Eos_bf_tabul::n2_tab
Tbl * n2_tab
Table of .
Definition eos_bifluid.h:1482
Lorene::Eos_bf_tabul::mu2_max
double mu2_max
Upper boundary of the chemical-potential interval (fluid 2 = p).
Definition eos_bifluid.h:1464
Lorene::Eos_bf_tabul::press_ent_p
virtual double press_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes the pressure from the baryonic log-enthalpies and the relative velocity.
Definition eos_bf_tabul.C:1087
Lorene::Eos_bf_tabul::dpsddelta_car_tab
Tbl * dpsddelta_car_tab
Table of .
Definition eos_bifluid.h:1488
Lorene::Eos_bf_tabul::identify
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Definition eos_bf_file.C:105
Lorene::Eos_bf_tabul::operator=
void operator=(const Eos_bf_tabul &)
Assignment to another Eos_bf_tabul.
Definition eos_bf_tabul.C:177
Lorene::Eos_bf_tabul::dn2sddelta_car_tab
Tbl * dn2sddelta_car_tab
Table of .
Definition eos_bifluid.h:1494
Lorene::Eos_bf_tabul::nbar_ent_p
virtual bool nbar_ent_p(const double ent1, const double ent2, const double delta2, double &nbar1, double &nbar2) const
Computes both baryon densities from the log-enthalpies.
Definition eos_bf_tabul.C:943
Lorene::Eos_bf_tabul::delta_car_max
double delta_car_max
Upper boundary of the relative velocity interval.
Definition eos_bifluid.h:1452
Lorene::Eos_bf_tabul::authors
string authors
Authors.
Definition eos_bifluid.h:1446
Lorene::Eos_bf_tabul::operator!=
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference).
Definition eos_bf_tabul.C:267
Lorene::Eos_bf_tabul::calcule_interpol
void calcule_interpol(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, Cmp &K_nn, Cmp &K_np, Cmp &K_pp, Cmp &alpha_eos, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
Definition eos_bf_tabul.C:757
Lorene::Eos_bf_tabul::press_tab
Tbl * press_tab
Table of .
Definition eos_bifluid.h:1476
Lorene::Eos_bf_tabul::~Eos_bf_tabul
virtual ~Eos_bf_tabul()
Destructor.
Definition eos_bf_tabul.C:146
Lorene::Eos_bf_tabul::tablename
string tablename
Name of the file containing the tabulated data (be careful, Eos_bifluid uses char*).
Definition eos_bifluid.h:1443
Lorene::Eos_bf_tabul::mu2_min
double mu2_min
Lower boundary of the chemical-potential interval (fluid 2 = p).
Definition eos_bifluid.h:1461
Lorene::Eos_bf_tabul::mu1_min
double mu1_min
Lower boundary of the chemical-potential interval (fluid 1 = n).
Definition eos_bifluid.h:1455
Lorene::Eos_bf_tabul::alpha_ent_p
virtual double alpha_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes alpha, the derivative of the total energy density with respect to from the baryonic log-ent...
Definition eos_bf_tabul.C:1177
Lorene::Eos_bf_tabul::operator==
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality).
Definition eos_bf_tabul.C:244
Lorene::Eos_bf_tabul::delta_car_min
double delta_car_min
Lower boundary of the relative velocity interval.
Definition eos_bifluid.h:1449
Lorene::Eos_bf_tabul::mu1_max
double mu1_max
Upper boundary of the chemical-potential interval (fluid 1 = n).
Definition eos_bifluid.h:1458
Lorene::Eos_bf_tabul::interpol_2d_Tbl3d
void interpol_2d_Tbl3d(const int i, const int j, const int k, const Tbl &ytab, const Tbl &ztab, const Tbl &ftab, const Tbl &dfdytab, const Tbl &dfdztab, const Tbl &d2fdydztab, const double y, const double z, double &f, double &dfdy, double &dfdz) const
Routine used by interpol_3d_bifluid to perform the 2D interpolation on the chemical potentials on eac...
Definition eos_bf_tabul.C:1839
Lorene::Eos_bf_tabul::trans2Eos
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
Definition eos_bf_tabul.C:1948
Lorene::Eos_bf_tabul::nbar_ent_p2
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
Definition eos_bf_tabul.C:1034
Lorene::Eos_bf_tabul::ener_nbar_p
virtual double ener_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the total energy density from the baryonic densities and the relative velocity.
Definition eos_bf_tabul.C:1060
Lorene::Eos_bf_tabul::press_nbar_p
virtual double press_nbar_p(const double nbar1, const double nbar2, const double delta2) const
Computes the pressure from the baryonic densities and the relative velocity.
Definition eos_bf_tabul.C:1073
Lorene::Eos_bf_tabul::n1_tab
Tbl * n1_tab
Table of .
Definition eos_bifluid.h:1479
Lorene::Eos_bf_tabul::d2psdmu1dmu2_tab
Tbl * d2psdmu1dmu2_tab
Table of .
Definition eos_bifluid.h:1485
Lorene::Eos_bf_tabul::get_K22
virtual double get_K22(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy/(baryonic density 2) .
Definition eos_bf_tabul.C:1247
Lorene::Eos_bf_tabul::read_table
void read_table()
Reads the file containing the table and initializes the arrays mu1_tab, mu2_tab, delta_car_tab,...
Definition eos_bf_tabul.C:279
Lorene::Eos_bf_tabul::interpol_3d_bifluid
void interpol_3d_bifluid(const double delta2, const double mu1, const double mu2, double &press, double &nbar1, double &nbar2, double &alpha) const
General method computing the pressure, both baryon densities and alpha from the values of both chemic...
Definition eos_bf_tabul.C:1298
Lorene::Eos_bf_tabul::delta_car_tab
Tbl * delta_car_tab
Table of .
Definition eos_bifluid.h:1473
Lorene::Eos_bf_tabul::operator>>
virtual ostream & operator>>(ostream &) const
Operator >>.
Definition eos_bf_tabul.C:229
Lorene::Eos_bf_tabul::get_K12
virtual double get_K12(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to .
Definition eos_bf_tabul.C:1270
Lorene::Eos_bf_tabul::get_K11
virtual double get_K11(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to (baryonic density 1) .
Definition eos_bf_tabul.C:1223
Lorene::Eos_bf_tabul::mu2_tab
Tbl * mu2_tab
Table of where .
Definition eos_bifluid.h:1470
Lorene::Eos_bf_tabul::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos_bf_tabul.C:218
Lorene::Eos_bf_tabul::ener_ent_p
virtual double ener_ent_p(const double ent1, const double ent2, const double delta_car) const
Computes the total energy density from the baryonic log-enthalpies and the relative velocity.
Definition eos_bf_tabul.C:1142
Lorene::Eos_bf_tabul::nbar_ent_p1
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
Definition eos_bf_tabul.C:1009
Lorene::Eos_bf_tabul::dn1sddelta_car_tab
Tbl * dn1sddelta_car_tab
Table of .
Definition eos_bifluid.h:1491
Lorene::Eos_bf_tabul::Eos_bf_tabul
Eos_bf_tabul(const char *name_i, const char *table, const char *path, double mass1, double mass2)
Standard constructor.
Definition eos_bf_tabul.C:84
Lorene::Eos_bf_tabul::mu1_tab
Tbl * mu1_tab
Table of where .
Definition eos_bifluid.h:1467
Lorene::Eos_bifluid::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition eos_bifluid.C:250
Lorene::Eos_bifluid::m_1
double m_1
Individual particle mass [unit: ].
Definition eos_bifluid.h:191
Lorene::Eos_bifluid::name
string name
EOS name.
Definition eos_bifluid.h:186
Lorene::Eos_bifluid::identify
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
Lorene::Eos_bifluid::operator=
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
Definition eos_bifluid.C:237
Lorene::Eos_bifluid::m_2
double m_2
Individual particle mass [unit: ].
Definition eos_bifluid.h:196
Lorene::Eos_bifluid::Eos_bifluid
Eos_bifluid()
Standard constructor.
Definition eos_bifluid.C:154
Lorene::Eos_poly
Polytropic equation of state (relativistic case).
Definition eos.h:809
Lorene::Eos
Equation of state base class.
Definition eos.h:206
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::dim
Dim_tbl dim
Number of dimensions, size,...
Definition tbl.h:172
Lorene::Tbl::set_etat_zero
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition tbl.C:350
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::c_est_pas_fait
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Definition c_est_pas_fait.C:76
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::z
Coord z
z coordinate centered on the grid
Definition map.h:740
Lorene::y
Coord y
y coordinate centered on the grid
Definition map.h:739
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738
Unites
Standard units of space, time and mass.