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
30char eos_bf_tabul_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_tabul.C,v 1.3 2015/06/11 14:41:59 a_sourie Exp $" ;
31
32/*
33 * $Id: eos_bf_tabul.C,v 1.3 2015/06/11 14:41:59 a_sourie Exp $
34 * $Log: eos_bf_tabul.C,v $
35 * Revision 1.3 2015/06/11 14:41:59 a_sourie
36 * Corrected minor bug
37 *
38 * Revision 1.2 2015/06/11 13:50:19 j_novak
39 * Minor corrections
40 *
41 * Revision 1.1 2015/06/10 14:39:17 a_sourie
42 * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bf_tabul.C,v 1.3 2015/06/11 14:41:59 a_sourie Exp $
46 *
47 */
48
49// headers C // Are they all useful ?
50#include <cstdlib>
51#include <cstring>
52#include <cmath>
53#include <cstdlib>
54
55
56// Headers Lorene // Are they all useful ??
57#include "eos_bifluid.h"
58#include "param.h"
59#include "eos.h"
60#include "tbl.h"
61#include "utilitaires.h"
62#include "unites.h"
63#include "cmp.h"
64#include "nbr_spx.h"
65
66namespace Lorene {
67
68void interpol_linear(const Tbl& xtab, const Tbl& ytab,
69 double x, int& i, double& y) ;
70
71void interpol_herm(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
72 double x, int& i, double& y, double& dy) ;
73
74void interpol_mixed_3d(const Tbl&, const Tbl&, const Tbl&, const Tbl&,
75 const Tbl&, const Tbl&, const Tbl&,
76 double, double, double, double&, double&, double&) ;
77
78
79void interpol_mixed_3d_new(double m1, double m2, const Tbl& xtab, const Tbl& ytab, const Tbl& ztab,
80 const Tbl& ftab, const Tbl& dfdytab, const Tbl& dfdztab, const Tbl& d2fdydztab,
81 const Tbl& dlpsddelta_car, const Tbl& d2lpsdlent1ddelta_car, const Tbl& d2lpsdlent2ddelta_car,
82 const Tbl& mu2_P, const Tbl& n_p_P, const Tbl& press_P,
83 const Tbl& mu1_N, const Tbl& n_n_N, const Tbl& press_N,
84 const Tbl& delta_car_n0, const Tbl& mu1_n0, const Tbl& mu2_n0,
85 const Tbl& delta_car_p0, const Tbl& mu1_p0, const Tbl& mu2_p0,
86 double x, double y, double z, double& f, double& dfdy, double& dfdz, double& alpha) ;
87
88void interpolation_brutale(double x,double y, double z, // localisation
89 // sommet vitesse inferieure, mun inf, mup_inf
90 double delta111, double delta211,
91 double mu1_111, double mu1_121, double mu2_111, double mu2_112, double mu1_211, double mu1_221, double mu2_211, double mu2_212,
92 double p_111, double p_121, double p_112, double p_122, double p_211, double p_221, double p_212, double p_222,
93 double n1_111, double n1_121, double n1_112, double n1_122, double n1_211, double n1_221, double n1_212, double n1_222,
94 double n2_111, double n2_121, double n2_112, double n2_122, double n2_211, double n2_221, double n2_212, double n2_222,
95 double cross_111, double cross_121, double cross_112, double cross_122, double cross_211, double cross_221, double cross_212, double cross_222,
96 double& f, double& dfdy, double& dfdz) ;
97
98void interpolation_brutale_mod(double x,double y, double z, // localisation
99 // sommet vitesse inferieure, mun inf, mup_inf
100 double delta111, double delta211,
101 double mu1_111, double mu1_121, double mu2_111, double mu2_112, double mu1_211, double mu1_221, double mu2_211, double mu2_212,
102 double p_111, double p_121, double p_112, double p_122, double p_211, double p_221, double p_212, double p_222,
103 double n1_111, double n1_121, double n1_112, double n1_122, double n1_211, double n1_221, double n1_212, double n1_222,
104 double n2_111, double n2_121, double n2_112, double n2_122, double n2_211, double n2_221, double n2_212, double n2_222,
105 double& f, double& dfdy, double& dfdz) ;
106
107
108/*void interpol_mixed_3d_test3(const Tbl&, const Tbl&, const Tbl&, const Tbl&,
109 const Tbl&, const Tbl&, const Tbl&,
110 double, double, double, double&, double&, double&) ;
111
112void interpol_mixed_3d_test2(const Tbl&, const Tbl&, const Tbl&, const Tbl&,
113 const Tbl&, const Tbl&, const Tbl&, const Tbl&, const Tbl&,
114 double, double, double, double&, double&, double&) ;
115*/
116void interpol_mixed_3d_mod(const Tbl&, const Tbl&, const Tbl&, const Tbl&,
117 const Tbl&, const Tbl&,
118 double, double, double, double&, double&, double&) ;
119
120
121 //----------------------------//
122 // Constructors //
123 //----------------------------//
124
125
126// Standard constructor
127// --------------------
128Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* table,
129 const char* path, double mass1, double mass2) :
130Eos_bifluid(name_i, mass1, mass2)
131{
132 tablename = path ;
133 tablename += "/" ;
134 tablename += table ;
135
136 read_table() ;
137}
138
139// Standard constructor with full filename
140// ---------------------------------------
141Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* file_name,
142 double mass1, double mass2) :
143Eos_bifluid(name_i, mass1, mass2)
144{
145 tablename = file_name ;
146
147 read_table() ;
148}
149
150// Constructor from binary file
151// ----------------------------
153Eos_bifluid(fich)
154{
155 char tmp_string[160] ;
156 fread(tmp_string, sizeof(char), 160, fich) ;
157 tablename = tmp_string ;
158
159 read_table() ;
160}
161
162// // Constructor from a formatted file
163// // ---------------------------------
164 Eos_bf_tabul::Eos_bf_tabul(ifstream& fich, const char* table) :
165 Eos_bifluid(fich)
166 {
167 fich >> tablename ;
168 tablename += "/" ;
169 tablename += table ;
170 read_table() ;
171 }
172
173 Eos_bf_tabul::Eos_bf_tabul(ifstream& fich) :
174 Eos_bifluid(fich)
175 {
176 fich >> tablename ;
177 read_table() ;
178 }
179
180
181 //--------------//
182 // Destructor //
183 //--------------//
184
186
187 delete logent1 ;
188 delete logent2 ;
189 delete delta_car ;
190 delete logp ;
191 delete dlpsdlent1 ;
192 delete dlpsdlent2 ;
193 delete d2lpsdlent1dlent2 ;
194 delete d2lpsdlent1dlent1 ;
195 delete d2lpsdlent2dlent2 ;
196 delete dlpsddelta_car ;
197 delete d2lpsdlent1ddelta_car ;
198 delete d2lpsdlent2ddelta_car ;
199 delete d3lpsdlent1dlent2ddelta_car ; // if necessary
200 delete delta_car_n0 ;
201 delete mu1_n0 ;
202 delete mu2_n0 ;
203 delete delta_car_p0 ;
204 delete mu1_p0 ;
205 delete mu2_p0 ;
206 delete mu1_N ;
207 delete n_n_N ;
208 delete press_N ;
209 delete mu2_P ;
210 delete n_p_P ;
211 delete press_P ;
212
213}
214
215 //--------------//
216 // Assignment //
217 //--------------//
218
220
221 // Assignment of quantities common to all the derived classes of Eos_bifluid
223
224 tablename = eosi.tablename;
225 ent1_min = eosi.ent1_min ;
226 ent1_max = eosi.ent1_max ;
227 ent2_min = eosi.ent2_min ;
228 ent2_max = eosi.ent2_max ;
231 logent1 = eosi.logent1 ;
232 logent2 = eosi.logent2 ;
233 delta_car = eosi.delta_car ;
234 logp = eosi.logp ;
235 dlpsdlent1 = eosi.dlpsdlent1 ;
236 dlpsdlent2 = eosi.dlpsdlent2 ;
242 delta_car_n0= eosi.delta_car_n0 ;
243 mu1_n0= eosi.mu1_n0 ;
244 mu2_n0= eosi.mu2_n0 ;
245 delta_car_p0= eosi.delta_car_p0 ;
246 mu1_p0 = eosi.mu1_p0;
247 mu2_p0 = eosi.mu2_p0 ;
248 mu1_N = eosi.mu1_N ;
249 n_n_N = eosi.n_n_N;
250 press_N = eosi.press_N;
251 mu2_P = eosi.mu2_P ;
252 n_p_P = eosi.n_p_P;
253 press_P = eosi.press_P;
254
255}
256
257
258 //------------//
259 // Outputs //
260 //------------//
261
262void Eos_bf_tabul::sauve(FILE* fich) const {
263
264 Eos_bifluid::sauve(fich) ;
265
266 char tmp_string[160] ;
267 strcpy(tmp_string, tablename.c_str()) ;
268 fwrite(tmp_string, sizeof(char), 160, fich) ;
269
270}
271
272
273ostream& Eos_bf_tabul::operator>>(ostream & ost) const {
274
275 ost <<
276 "EOS of class Eos_bf_tabul : tabulated EOS depending on three variables (relative velocity and enthalpies of neutrons and protons)."
277 << '\n' ;
278 ost << "Table read from " << tablename << endl ;
279
280 return ost ;
281}
282
283 //------------------------//
284 // Comparison operators //
285 //------------------------//
286
287
288bool Eos_bf_tabul::operator==(const Eos_bifluid& eos_i) const {
289
290 bool resu = true ;
291
292 if ( eos_i.identify() != identify() ) {
293 cout << "The second EOS is not of type Eos_bf_tabul !" << endl ;
294 resu = false ;
295 }
296
297 else {
298 const Eos_bf_tabul& eos = dynamic_cast<const Eos_bf_tabul&>( eos_i ) ;
299
300 if (eos.tablename != tablename){
301 cout <<
302 "The two Eos_bf_tabul have different names and not refer to the same tables! " << endl ;
303 resu = false ;
304 }
305 }
306
307 return resu ;
308
309}
310
311
312bool Eos_bf_tabul::operator!=(const Eos_bifluid& eos_i) const {
313
314 return !(operator==(eos_i)) ;
315
316}
317
318
319 //------------------------//
320 // Reading of the tables //
321 //------------------------//
322
324
325 using namespace Unites ;
326 double m_b_si = 1.66E-27;
327
328 //-------------------------------------------
329 //---------- Table à deux fluides -----------
330 //-------------------------------------------
331
332 ifstream fich(tablename.data()) ;
333 //DEBUG : cout << tablename << endl;
334 if (!fich) {
335 cerr << "Eos_bf_tabul::read_table(): " << endl ;
336 cerr << "Problem in opening the EOS file!" << endl ;
337 cerr << "While trying to open " << tablename << endl ;
338 cerr << "Aborting..." << endl ;
339 abort() ;
340 }
341
342 fich.ignore(1000, '\n') ; // Jump over the first header
343 fich.ignore(2) ;
344 getline(fich, authors, '\n') ;
345 //cout << authors << endl ;
346 for (int i=0; i<5; i++) { // jump over the file
347 fich.ignore(1000, '\n') ; // header
348 } //
349
350 int nbp ;
351 fich >> nbp ; fich.ignore(1000, '\n') ; // number of data
352 //cout << nbp << endl ;
353
354 int n_delta, n_mu1, n_mu2; // number of different values for : delta_car, mu_n, mu_p
355 fich >> n_delta; fich.ignore(1000, '\n') ;
356 fich >> n_mu1; fich.ignore(1000, '\n') ;
357 fich >> n_mu2; fich.ignore(1000, '\n') ;
358 //cout << "n_delta : " << n_delta << endl ;
359 //cout << "n_mu1 : " << n_mu1 << endl ;
360 //cout << "n_mu2 : " <<n_mu2 << endl ;
361
362 if (nbp<=0) {
363 cerr << "Eos_bf_tabul::read_table(): " << endl ;
364 cerr << "Wrong value for the number of lines!" << endl ;
365 cerr << "nbp = " << nbp << endl ;
366 cerr << "Aborting..." << endl ;
367 abort() ;
368 }
369 if (n_delta<=0) {
370 cerr << "Eos_bf_tabul::read_table(): " << endl ;
371 cerr << "Wrong value for the number of values of delta_car!" << endl ;
372 cerr << "n_delta = " << n_delta << endl ;
373 cerr << "Aborting..." << endl ;
374 abort() ;
375 }
376 if (n_mu1<=0) {
377 cerr << "Eos_bf_tabul::read_table(): " << endl ;
378 cerr << "Wrong value for the number of values of mu_n!" << endl ;
379 cerr << "n_mu1 = " << n_mu1 << endl ;
380 cerr << "Aborting..." << endl ;
381 abort() ;
382 }
383 if (n_mu2<=0) {
384 cerr << "Eos_bf_tabul::read_table(): " << endl ;
385 cerr << "Wrong value for the number of values of mu_p!" << endl ;
386 cerr << "n_mu2 = " << n_mu2 << endl ;
387 cerr << "Aborting..." << endl ;
388 abort() ;
389 }
390 if (n_mu2*n_mu1*n_delta != nbp ) {
391 cerr << "Eos_bf_tabul::read_table(): " << endl ;
392 cerr << "Wrong value for the number of lines!" << endl ;
393 cerr << "Aborting..." << endl ;
394 abort() ;
395 }
396
397 for (int i=0; i<3; i++) { // jump over the table
398 fich.ignore(1000, '\n') ; // description
399 }
400
401 logent1 = new Tbl(n_delta, n_mu1, n_mu2) ;
402 logent2 = new Tbl(n_delta, n_mu1, n_mu2) ;
403 delta_car = new Tbl(n_delta, n_mu1, n_mu2) ;
404 logp = new Tbl(n_delta, n_mu1, n_mu2) ;
405 dlpsdlent1 = new Tbl(n_delta, n_mu1, n_mu2) ;
406 dlpsdlent2 = new Tbl(n_delta, n_mu1, n_mu2) ;
407 d2lpsdlent1dlent2 = new Tbl(n_delta, n_mu1, n_mu2) ;
408 dlpsddelta_car = new Tbl(n_delta, n_mu1, n_mu2) ;
409 d2lpsdlent1ddelta_car = new Tbl(n_delta, n_mu1, n_mu2) ;
410 d2lpsdlent2ddelta_car = new Tbl(n_delta, n_mu1, n_mu2) ;
411 d2lpsdlent1dlent1 = new Tbl(n_delta, n_mu1, n_mu2) ;
412 d2lpsdlent2dlent2 = new Tbl(n_delta, n_mu1, n_mu2) ;
413
417 logp->set_etat_qcq() ;
426
427 //------------------------------------------------------
428 // We have to choose the right unites (SI, cgs , LORENE)
429 //------------------------------------------------------
430 // sor far, all the calculations are done with Lorene units
431
432 //Grandeurs issues de la table à deux fluides (unités physiques adaptées)
433 double mu1_MeV, mu2_MeV, n1_fm3, n2_fm3;
434 double Knp_Mev_2, press_MeV_fm3;
435 double d2press_dmu1dmu2_MeV_fm3, dn1_ddelta_car_fm3, dn2_ddelta_car_fm3;
436
437 //Grandeurs stockées (unité Lorene)
438 double delta_car_adim, mu1_adim, mu2_adim, n1_01fm3, n2_01fm3, Knp_adim ;
439 double press_adim, dpress_ddelta_car_adim, dn1_ddelta_car_adim, dn2_ddelta_car_adim ;
440 double d2press_dmu1dmu2_adim;
441 double hbarc = 197.33 ; //en Mev*fm
442 double hbarc3 = hbarc * hbarc * hbarc ;
443
444
445 for (int j=0 ; j < n_delta ; j++) {
446 for (int k=0 ; k < n_mu1 ; k++) {
447 for (int l=0 ; l < n_mu2 ; l++) {
448
449 fich >> delta_car_adim ;
450 fich >> mu1_MeV ;
451 mu1_adim = mu1_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
452 fich >> mu2_MeV ;
453 mu2_adim = mu2_MeV * mev_si / (m_b_si * v_unit * v_unit ) ;
454 fich >> n1_fm3 ;
455 n1_01fm3 = 10. * n1_fm3 ;
456 fich >> n2_fm3 ;
457 n2_01fm3 = 10. * n2_fm3 ;
458 fich >> Knp_Mev_2 ;
459 Knp_adim = Knp_Mev_2 / ( m_b_si * v_unit * v_unit *10. )
460 * (mev_si *hbarc3 ) ;
461 dpress_ddelta_car_adim = - Knp_adim * n1_01fm3 * n2_01fm3
462 * pow( 1.-delta_car_adim, -1.5) / 2. ;
463 // cout << "neg" << dpress_ddelta_car_adim << endl ;
464 fich >> press_MeV_fm3 ;
465 press_adim = press_MeV_fm3 * mevpfm3 ;
466 fich >> d2press_dmu1dmu2_MeV_fm3 ;
467 d2press_dmu1dmu2_adim = d2press_dmu1dmu2_MeV_fm3
468 * (10. * m_b_si * v_unit * v_unit ) / mev_si ;
469 fich >> dn1_ddelta_car_fm3 ;
470 dn1_ddelta_car_adim = dn1_ddelta_car_fm3 * 10. ;
471 fich >> dn2_ddelta_car_fm3 ;
472 dn2_ddelta_car_adim = dn2_ddelta_car_fm3 * 10. ;
473
474 //To check the reading is well done :
475 //cout << j << " " << k << " " << l << " " << delta_car_adim << " "<< mu1_MeV << " "<< mu2_MeV << " " << Knp_Mev_2 << " " << dn2_ddelta_car_fm3 << endl;
476
477 fich.ignore(1000,'\n') ;
478
479 /* if ( (n1_fm3<0) || (n2_fm3<0) || (press_MeV_fm3 < 0)){
480 cout << "Eos_tabul::read_table(): " << endl ;
481 cout << "Negative value in table!" << endl ;
482 cout << "n_neutrons = " << n1_fm3 << ", n_protons " << n2_fm3 <<
483 ", p = " << press_MeV_fm3 << ", "<< endl ;
484 cout << "Aborting..." << endl ;
485 abort() ;
486 }*/ // il peut y avoir des densités négatives selon la table à deux fluides utilisés...
487
488
489
490 /* redefinition de la masse --> permet de redefinir le minimum en log enthalpie !
491 double ww1 = 1.;
492 double ww2 = 1.;
493 if ((k==0) && (l==0) ) {ww1 = mu1_adim; ww2 = mu2_adim;} ;
494 cout << setprecision(16) ;
495 cout << ww1 ;
496 mu1_adim = mu1_adim/ww1 ; //+1e-14 ;
497 mu2_adim = mu2_adim/ww2 ;
498 */
499
500 // on enregistre les grandeurs dans des tableaux sans passer par des logarithmes !!! --> si on garde il faudra choisir des meilleurs noms de tableaux !
501 logent1->set(j,k,l) = mu1_adim ;
502 logent2->set(j,k,l) = mu2_adim ;
503 delta_car->set(j,k,l) = delta_car_adim ;
504 // logp->set(j,k,l) = log(press_adim) ;
505 if (press_adim<=0) {press_adim = 0. ;}
506 logp->set(j,k,l) = press_adim ; // non en log
507 // dlpsdlent1->set(j,k,l) = n1_01fm3 /press_adim ;
508 dlpsdlent1->set(j,k,l) = n1_01fm3 ; // non en log
509 // dlpsdlent2->set(j,k,l) = n2_01fm3 /press_adim ;
510 dlpsdlent2->set(j,k,l) = n2_01fm3 ; // non en log
511
512 if ((n1_01fm3 < 1e-16) && (n2_01fm3 <1e-16 )) {
513 d2press_dmu1dmu2_adim = 0. ;
514 dpress_ddelta_car_adim = 0. ;
515 dn1_ddelta_car_adim =0. ;
516 dn2_ddelta_car_adim = 0. ;
517 }
518 // d2lpsdlent1dlent2->set(j,k,l) = d2press_dmu1dmu2_adim
519 // / press_adim - n1_01fm3 /press_adim * n2_01fm3 /press_adim ;
520 d2lpsdlent1dlent2->set(j,k,l) = d2press_dmu1dmu2_adim ; // non en log
521 dlpsddelta_car->set(j,k,l) = dpress_ddelta_car_adim ;
522 d2lpsdlent1ddelta_car->set(j,k,l) = dn1_ddelta_car_adim ;
523 d2lpsdlent2ddelta_car->set(j,k,l) = dn2_ddelta_car_adim ;
524 }
525
526 fich.ignore(1000, '\n') ; // BIEN VERIFIER QU ON SAUTE UNE LIGNE DANS TOUTES LES TABLES MAIS NORMALEMENT OUI !
527
528 }
529
530 fich.ignore(1000, '\n') ;
531
532 }
533 //abort();
534 delta_car_min = (*delta_car)(0, 0, 0) ;
535 delta_car_max = (*delta_car)(n_delta-1, 0, 0) ;
536
537 // si on passe par des log
538 /*ent1_min = pow( double(10), (*logent1)(0, 0, 0)) ;
539 ent1_max = pow( double(10), (*logent1)(0, n_mu1-1, 0)) ;
540 ent2_min = pow( double(10), (*logent2)(0, 0, 0)) ;
541 ent2_max = pow( double(10), (*logent2)(0, 0, n_mu2-1));
542 */
543
544 /*ent1_min = log((*logent1)(0, 0, 0)) ;
545 ent1_max = log((*logent1)(0, n_mu1-1, n_mu2-1)) ;
546 ent2_min = log((*logent2)(0, 0, 0)) ;
547 ent2_max = log((*logent2)(0, n_mu1-1, n_mu2-1));
548 */
549
550 // sans passer par les log !
551 ent1_min = (*logent1)(0, 0, 0) ; // en fait il s'agit de mu
552 ent1_max = (*logent1)(0, n_mu1-1, 0) ;
553 ent2_min = (*logent2)(0, 0, 0) ;
554 ent2_max = (*logent2)(0, 0, n_mu2-1);
555
556 //cout <<"DEBUG : ent1_min "<< ent1_min << endl ;
557 //cout <<"DEBUG : ent1_max "<< ent1_max << endl ;
558 //cout <<"DEBUG : ent2_min "<< ent2_min << endl ;
559 //cout <<"DEBUG : ent2_max "<< ent2_max << endl ;
560
561
562
563
564 //cout <<"DEBUG : ent1_min en MeV "<< ent1_min / mev_si * (m_b_si * v_unit * v_unit ) << endl ;
565 //cout <<"DEBUG : ent1_max en MeV "<< ent1_max / mev_si * (m_b_si * v_unit * v_unit ) << endl ;
566 //cout <<"DEBUG : ent2_min en MeV "<< ent2_min / mev_si * (m_b_si * v_unit * v_unit ) << endl ;
567 //cout <<"DEBUG : ent2_max en MeV "<< ent2_max / mev_si * (m_b_si * v_unit * v_unit )<< endl ;
568 //abort();
569 fich.close();
570
571
572 //---------------------------------------------------
573 //---------- Table à un fluide : fluide N -----------
574 //---------------------------------------------------
575
576 ifstream fichN ;
577 fichN.open("eos_bf_test_1_fluide_N.d") ;
578
579 if (!fichN) {
580 cerr << "Eos_bf_tabul::read_table(): " << endl ;
581 cerr << "Problem in opening the EOS file!" << endl ;
582 cerr << "While trying to open " << tablename << endl ;
583 cerr << "Aborting..." << endl ;
584 abort() ;
585 }
586
587 fichN.ignore(1000, '\n') ; // Jump over the first header
588 fichN.ignore(2) ;
589 fichN.ignore(1000, '\n') ;
590
591 for (int i=0; i<5; i++) {
592 fichN.ignore(1000, '\n') ;
593 }
594
595 int nbp_N ;
596 fichN >> nbp_N ; fichN.ignore(1000, '\n') ; // number of data
597 cout<< "nbp_N = " << nbp_N ;
598 int n_mu1_N; // number of different values for : delta_car, n_n, n_p // SEULS n_n sert le reste est inutile et pourra etre enlever dans une version finale
599 fichN >> n_mu1_N;fichN.ignore(1000, '\n') ;
600
601 //cout << n_mu1_N endl ;
602
603 if (nbp_N<=0) {
604 cerr << "Eos_bf_tabul::read_table(): " << endl ;
605 cerr << "Wrong value for the number of lines!" << endl ;
606 cerr << "nbp = " << nbp << endl ;
607 cerr << "Aborting..." << endl ;
608 abort() ;
609 }
610 if (n_mu1_N<=0) {
611 cerr << "Eos_bf_tabul::read_table(): " << endl ;
612 cerr << "Wrong value for the number of values of mu_n!" << endl ;
613 cerr << "n_mu1 = " << n_mu1 << endl ;
614 cerr << "Aborting..." << endl ;
615 abort() ;
616 }
617 for (int i=0; i<3; i++) { // jump over the table
618 fichN.ignore(1000, '\n') ; // description
619 }
620
621 mu1_N = new Tbl(n_mu1_N) ;
622 n_n_N = new Tbl(n_mu1_N) ;
623 press_N = new Tbl(n_mu1_N) ;
624
625 mu1_N ->set_etat_qcq() ;
626 n_n_N->set_etat_qcq() ;
627 press_N ->set_etat_qcq() ;
628
629 //Grandeurs issues de la table à un fluide (unités physiques adaptées)
630 double mu1_MeV_N, n1_fm3_N,press_MeV_fm3_N;
631
632 //Grandeurs stockées (unité Lorene)
633 double mu1_adimN, n1_01fm3N,press_adimN;
634
635 for (int k=0 ; k < n_mu1_N ; k++) {
636
637 fichN >> mu1_MeV_N ;
638 mu1_adimN = mu1_MeV_N * mev_si / (m_b_si * v_unit * v_unit ) ;
639 fichN >> n1_fm3_N ;
640 n1_01fm3N = 10. * n1_fm3_N ;
641 fichN >> press_MeV_fm3_N;
642 press_adimN = press_MeV_fm3_N * mevpfm3 ;
643 fichN.ignore(1000,'\n') ;
644
645 //cout << k << " " << press_MeV_fm3_N << endl;
646
647 if ( (n1_01fm3N<0) || (press_adimN < 0)){
648 cout << "Eos_tabul::read_table(): " << endl ;
649 cout << "Negative value in table!" << endl ;
650 cout << "n_neutrons = " << n1_01fm3N <<
651 ", p = " << press_adimN << ", "<< endl ;
652 cout << "Aborting..." << endl ;
653 abort() ;
654 }
655
656 mu1_N ->set(k) = mu1_adimN ;
657 n_n_N->set(k) = n1_01fm3N ;
658 press_N ->set(k) = press_adimN ;
659
660 }
661 //abort();
662 //cout << setprecision(16) ;
663 //cout << endl ;
664 //cout << m_b_si * v_unit * v_unit / mev_si << endl ;
665 //cout << m_b_si << endl ;
666 //cout << v_unit << endl ;
667 //cout << mev_si << endl ;
668 //cout << rhonuc_si << endl ;
669 //abort() ;
670 fichN.close();
671
672
673
674 //---------------------------------------------------
675 //---------- Table à un fluide : fluide P -----------
676 //---------------------------------------------------
677
678
679 ifstream fichP ;
680 fichP.open("eos_bf_test_1_fluide_P.d") ;
681
682 if (!fichP) {
683 cerr << "Eos_bf_tabul::read_table(): " << endl ;
684 cerr << "Problem in opening the EOS file!" << endl ;
685 cerr << "While trying to open " << tablename << endl ;
686 cerr << "Aborting..." << endl ;
687 abort() ;
688 }
689
690 fichP.ignore(1000, '\n') ; // Jump over the first header
691 fichP.ignore(2) ;
692 fichP.ignore(1000, '\n') ;
693
694 for (int i=0; i<5; i++) { // jump over the file
695 fichP.ignore(1000, '\n') ; // header
696 } //
697
698 int nbp_P ;
699 fichP >> nbp_P ; fichP.ignore(1000, '\n') ; // number of data
700
701 int n_mu2_P; // number of different values for : delta_car, n_n, n_p // SEUL n_mu2_P a un sens !!!!!!!!!!!! ENLEVER LES AUTRES DANS UNE FUTURE VERSION
702 fichP >> n_mu2_P;fichP.ignore(1000, '\n') ;
703
704 if (nbp_P<=0) {
705 cerr << "Eos_bf_tabul::read_table(): " << endl ;
706 cerr << "Wrong value for the number of lines!" << endl ;
707 cerr << "nbp = " << nbp << endl ;
708 cerr << "Aborting..." << endl ;
709 abort() ;
710 }
711 if (n_mu2_P<=0) {
712 cerr << "Eos_bf_tabul::read_table(): " << endl ;
713 cerr << "Wrong value for the number of values of mu_p!" << endl ;
714 cerr << "n_mu2 = " << n_mu2 << endl ;
715 cerr << "Aborting..." << endl ;
716 abort() ;
717 }
718
719 for (int i=0; i<3; i++) { // jump over the table
720 fichP.ignore(1000, '\n') ; // description
721 }
722
723 mu2_P = new Tbl(n_mu2_P) ;
724 n_p_P = new Tbl(n_mu2_P) ;
725 press_P = new Tbl(n_mu2_P) ;
726
727 mu2_P ->set_etat_qcq() ;
728 n_p_P->set_etat_qcq() ;
729 press_P ->set_etat_qcq() ;
730
731 //Grandeurs issues de la table à un fluide (unités physiques adaptées)
732 double mu2_MeV_P, n2_fm3_P,press_MeV_fm3_P;
733
734 //Grandeurs stockées (unité Lorene)
735 double mu2_adimP, n2_01fm3P, press_adimP;
736
737 for (int l=0 ; l < n_mu2_P ; l++) {
738
739 fichP >> mu2_MeV_P ;
740 mu2_adimP = mu2_MeV_P * mev_si / (m_b_si * v_unit * v_unit ) ;
741 fichP >> n2_fm3_P ;
742 n2_01fm3P = 10. * n2_fm3_P ;
743 fichP >> press_MeV_fm3_P;
744 press_adimP = press_MeV_fm3_P * mevpfm3 ;
745
746 //cout << l << " " << press_MeV_fm3_P << endl;
747
748 fichP.ignore(1000,'\n') ;
749 if ( (n2_01fm3P<0) || (press_adimP < 0)){
750 cout << "Eos_tabul::read_table(): " << endl ;
751 cout << "Pegative value in table!" << endl ;
752 cout <<", n_protons " << n2_01fm3P <<
753 ", p = " << press_adimP << ", "<< endl ;
754 cout << "Aborting..." << endl ;
755 abort() ;
756 }
757
758 mu2_P ->set(l) = mu2_adimP ;
759 n_p_P->set(l) = n2_01fm3P ;
760 press_P ->set(l) = press_adimP ;
761
762 }
763 //abort();
764 fichP.close();
765
766
767 //------------------------------------------------------------------
768 //---------- COURBE LIMITE TABLE A DEUX FLUIDES : np = 0 -----------
769 //------------------------------------------------------------------
770
771 ifstream fich1 ;
772 fich1.open("np=0.dat") ;
773 // table classee en mun croissant !
774 if (!fich1) {
775 cerr << "Eos_bf_tabul::read_table(): " << endl ;
776 cerr << "Problem in opening the EOS file!" << endl ;
777 cerr << "While trying to open " << tablename << endl ;
778 cerr << "Aborting..." << endl ;
779 abort() ;
780 }
781 int n_delta_n0, n_mu1_n0;
782 fich1 >> n_delta_n0;fich1.ignore(1000, '\n') ;
783 fich1 >> n_mu1_n0;fich1.ignore(1000, '\n') ;
784 fich1.ignore(1000, '\n') ; // Jump over the first header
785 //cout << n_delta_n0 << " " << n_mu1_n0 << endl;
786 delta_car_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
787 mu1_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
788 mu2_n0 = new Tbl(n_delta_n0, n_mu1_n0) ;
789
790 delta_car_n0 ->set_etat_qcq() ;
791 mu1_n0->set_etat_qcq() ;
792 mu2_n0->set_etat_qcq() ;
793
794 double delta_car_nn0, mu1_MeV_nn0, mu2_MeV_nn0;
795
796 for (int o = 0; o < n_delta_n0 ; o++ ) {
797 for (int p = 0 ; p < n_mu1_n0 ; p++) {
798
799 fich1 >> delta_car_nn0 ;
800 fich1 >> mu1_MeV_nn0 ;
801 fich1 >> mu2_MeV_nn0 ;
802
803 fich1.ignore(1000,'\n') ;
804 // cout << o << " " << p << " " << delta_car_nn0 << " " << mu1_MeV_nn0 << " " << mu2_MeV_nn0<< endl ;
805 delta_car_n0 ->set(o,p) = delta_car_nn0;
806 mu1_n0 ->set(o,p) = mu1_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
807 mu2_n0 ->set(o,p) = mu2_MeV_nn0 * mev_si / (m_b_si * v_unit * v_unit ) ;
808
809 }
810 fich1.ignore(1000,'\n') ;
811 }
812 //abort();
813 fich1.close();
814
815 //------------------------------------------------------------------
816 //---------- COURBE LIMITE TABLE A DEUX FLUIDES : nn = 0 -----------
817 //------------------------------------------------------------------
818
819 ifstream fich2 ;
820 fich2.open("nn=0.dat") ;
821 // table classe en mup croissant ---------------> attention
822
823 if (!fich2) {
824 cerr << "Eos_bf_tabul::read_table(): " << endl ;
825 cerr << "Problem in opening the EOS file!" << endl ;
826 cerr << "While trying to open " << tablename << endl ;
827 cerr << "Aborting..." << endl ;
828 abort() ;
829 }
830 int n_delta_p0, n_mu2_p0;
831 fich2 >> n_delta_p0;fich2.ignore(1000, '\n') ;
832 fich2 >> n_mu2_p0;fich2.ignore(1000, '\n') ;
833 fich2.ignore(1000, '\n') ; // Jump over the first header
834 //cout << n_delta_p0 << " " << n_mu2_p0 << endl;
835 delta_car_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
836 mu1_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
837 mu2_p0 = new Tbl(n_delta_p0, n_mu2_p0) ;
838
839 delta_car_p0 ->set_etat_qcq() ;
840 mu1_p0->set_etat_qcq() ;
841 mu2_p0 ->set_etat_qcq() ;
842
843 double delta_car_np0, mu1_MeV_np0, mu2_MeV_np0;
844
845 for (int o = 0; o < n_delta_p0 ; o++ ) {
846 for (int p = 0 ; p < n_mu2_p0 ; p++) {
847
848 fich2 >> delta_car_np0 ;
849 fich2 >> mu1_MeV_np0 ;
850 fich2 >> mu2_MeV_np0 ;
851
852 fich2.ignore(1000,'\n') ;
853 // cout << o << " " << p << " " << delta_car_np0 << " " << mu1_MeV_np0 << " " << mu2_MeV_np0<< endl ;
854 delta_car_p0 ->set(o,p) = delta_car_np0;
855 mu1_p0 ->set(o,p) = mu1_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
856 mu2_p0 ->set(o,p) = mu2_MeV_np0 * mev_si / (m_b_si * v_unit * v_unit ) ;
857
858 //cout << o << " " << p << " " << mu2_MeV_np0 << endl;
859 }
860 fich2.ignore(1000,'\n') ;
861 }
862 //abort();
863 fich2.close();
864
865}
866
867
868 //------------------------------//
869 // Computational routines //
870 //------------------------------//
871
872
873// Complete computational routine giving all thermo variables
874//-----------------------------------------------------------
875
876void Eos_bf_tabul::calcule_interpol(const Cmp& ent1, const Cmp& ent2,
877 const Cmp& delta2, Cmp& nbar1, Cmp& nbar2,
878 Cmp& ener, Cmp& press, Cmp& K_nn, Cmp& K_np, Cmp& K_pp,
879 int nzet, int l_min) const {
880
881 assert(ent1.get_etat() != ETATNONDEF) ;
882 assert(ent2.get_etat() != ETATNONDEF) ;
883 assert(delta2.get_etat() != ETATNONDEF) ;
884
885 const Map* mp = ent1.get_mp() ; // Mapping
886 assert(mp == ent2.get_mp()) ;
887 assert(mp == delta2.get_mp()) ;
888 assert(mp == ener.get_mp()) ;
889
890 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
891 nbar1.set_etat_zero() ;
892 nbar2.set_etat_zero() ;
893 ener.set_etat_zero() ;
894 press.set_etat_zero() ;
895 K_nn.set_etat_zero() ;
896 K_np.set_etat_zero() ;
897 K_pp.set_etat_zero() ;
898 return ;
899 }
900 nbar1.allocate_all() ;
901 nbar2.allocate_all() ;
902 ener.allocate_all() ;
903 press.allocate_all() ;
904 K_nn.allocate_all() ;
905 K_np.allocate_all() ;
906 K_pp.allocate_all() ;
907 const Mg3d* mg = mp->get_mg() ; // Multi-grid
908
909 int nz = mg->get_nzone() ; // total number of domains
910
911 // resu is set to zero in the other domains :
912
913 if (l_min > 0) {
914 nbar1.annule(0, l_min-1) ;
915 nbar2.annule(0, l_min-1) ;
916 ener.annule(0, l_min-1) ;
917 press.annule(0, l_min-1) ;
918 K_nn.annule(0, l_min-1) ;
919 K_np.annule(0, l_min-1) ;
920 K_pp.annule(0, l_min-1) ;
921 }
922
923 if (l_min + nzet < nz) {
924 nbar1.annule(l_min + nzet, nz - 1) ;
925 nbar2.annule(l_min + nzet, nz - 1) ;
926 ener.annule(l_min + nzet, nz - 1) ;
927 press.annule(l_min + nzet, nz - 1) ;
928 K_nn.annule(l_min + nzet, nz - 1) ;
929 K_np.annule(l_min + nzet, nz - 1) ;
930 K_pp.annule(l_min + nzet, nz - 1) ;
931 }
932
933 int np0 = mp->get_mg()->get_np(0) ;
934 int nt0 = mp->get_mg()->get_nt(0) ;
935 for (int l=l_min; l<l_min+nzet; l++) {
936 assert(mp->get_mg()->get_np(l) == np0) ;
937 assert(mp->get_mg()->get_nt(l) == nt0) ;
938 }
939
940
941 for (int k=0; k<np0; k++) {
942 for (int j=0; j<nt0; j++) {
943
944 for (int l=l_min; l< l_min + nzet; l++) {
945 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
946
947 double xx, xx1, xx2; // xx1 = H1 = ln(mu1/m1)
948 xx1 = ent1(l,k,j,i) ;
949 xx2 = ent2(l,k,j,i) ;
950 xx = delta2(l,k,j,i) ;
951
952
953
954/*
955 * AVANT MODIFICATION
956 *
957 // ZONE a zero fluide ----------------------------------------------------------------------------------------------
958 if ((exp(xx1) <= 1. ) && (exp(xx2) <= 1. ) ) {
959 nbar1.set(l,k,j,i) = 0. ;
960 nbar2.set(l,k,j,i) = 0. ;
961 press.set(l,k,j,i) = 0. ;
962 ener.set(l,k,j,i) = 0.;
963 K_nn.set(l,k,j,i) = 0.;
964 K_np.set(l,k,j,i) = 0.;
965 K_pp.set(l,k,j,i) = 0.;
966 }
967
968 // ZONE a deux fluides ou fluide seul -----------------------------------------------------------------------------
969 else {
970
971 double n1 = 0. ;
972 double n2 = 0.;
973 double p = 0. ;
974
975 inside = (!nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
976 nbar1.set(l,k,j,i) = n1 ;
977 nbar2.set(l,k,j,i) = n2 ;
978
979 //---si n1<0 -------------------------
980// if (n1 <0. )
981// {
982 // soit on met a zero
983// n1 = 0.;
984 // soit on fait une interpolation lineaire (au choix)
985 //int i_tri1 = 0;
986 //interpol_linear(*mu1_N, *n_1_N, mu1_adim, i_tri1 , n1) ;
987// }
988 //---si n1<0 -------------------------
989 //---si n2<0 -------------------------
990// if (n2 <0. )
991// {
992 // soit on met a zero
993// n2 = 0.;
994 // soit on fait une interpolation lineaire (au choix)
995 //int i_tri2 = 0;
996 //interpol_linear(*mu2_P, *n_2_P, mu2_adim, i_tri2 , n2) ;
997// }
998 //---si n2<0 -------------------------
999
1000 //---si p<0 -------------------------
1001// if (p <0. )
1002// {
1003 // soit on met a zero
1004// p = 0.;
1005 // soit on fait une interpolation lineaire (au choix) // pas terrible
1006 //int i_trip = 0;
1007 //interpol_linear(*mu2_P, *press_P, mu2_adim, i_trip , p) ;
1008// }
1009 //---si p<0 -------------------------
1010
1011 p = press_ent_p(xx1, xx2, xx) ;
1012 press.set(l,k,j,i) = p ;
1013 double mu1 = m_1 * exp (xx1 );
1014 double mu2 = m_2 * exp( xx2) ;
1015 ener.set(l,k,j,i) = mu1 * n1 + mu2 * n2 - p ;//ener_ent_p(xx1, xx2, xx) ;
1016 K_nn.set(l,k,j,i) = get_K11(xx, xx1, xx2);
1017 K_np.set(l,k,j,i) = get_K12(xx, xx1, xx2);
1018 K_pp.set(l,k,j,i) = get_K22(xx, xx1, xx2);
1019 }
1020
1021*
1022* AVANT MODIFICATION
1023*/
1024
1025 //---------------------------------------------------------------------------------------------------------------
1026 // plus d'appel à plein de fonctions, tout est fait d'un coup!
1027 //-----------------------------------------------------------------------------------------------------------------
1028
1029 if (xx < delta_car_min) {
1030 cout << "Eos_bf_tabul::calcule_tout : delta2 < delta_car_min !"
1031 << endl ;
1032 abort() ;
1033 }
1034 if (xx > delta_car_max) {
1035 cout << "Eos_bf_tabul::calcule_tout : delta2 > delta_car_max !"
1036 << endl ;
1037 abort() ;
1038 }
1039 if (m_1 * exp(xx1) > ent1_max) {
1040 cout << "Eos_bf_tabul::calcule_tout : ent1 > ent1_max !" << endl ;
1041 abort() ;
1042 }
1043 if (m_2 *exp(xx2) > ent2_max) {
1044 cout << "Eos_bf_tabul::calcule_tout : ent2 > ent2_max !" << endl ;
1045 abort() ;
1046 }
1047
1048 double n1 = 0 ;
1049 double n2 = 0 ;
1050 double pressure = 0 ;
1051 double alpha = 0 ;
1052 double K11 = 0 ;
1053 double K12 = 0 ;
1054 double K22 = 0 ;
1055
1056 double mu1 = m_1 * exp(xx1);
1057 double mu2 = m_2 * exp(xx2);
1058
1059 if ( (exp(xx1) < 1.) && (exp(xx2) < 1.) ) {
1060 n1 = 0. ;
1061 n2 = 0. ;
1062 pressure = 0.;
1063 alpha = 0 ;
1064 K11 = 0 ;
1065 K12 = 0 ;
1066 K22 = 0 ;
1067 }
1068 else {
1069
1070
1071
1072 double p_interpol ;
1073 double dpsdent1_interpol ;
1074 double dpsdent2_interpol ;
1075 double alpha_interpol ;
1076
1077 // cout << "DEBUG:" << xx << " " << mu1 << " " << mu2 << endl;
1078
1079 interpol_mixed_3d_new(m_1, m_2, *delta_car, *logent1, *logent2,
1082 *mu2_P, *n_p_P, *press_P,
1083 *mu1_N, *n_n_N, *press_N,
1084 *delta_car_n0, *mu1_n0, *mu2_n0,
1085 *delta_car_p0, *mu1_p0, *mu2_p0,
1086 xx, mu1, mu2,
1087 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha_interpol) ;
1088
1089 // en log :
1090// pressure = exp(p_interpol) ;
1091// n1 = dpsdent1_interpol * pressure ;
1092// n2 = dpsdent2_interpol * pressure ;
1093 n1 = dpsdent1_interpol ;
1094 n2 = dpsdent2_interpol ;
1095 pressure = p_interpol;
1096 alpha = alpha_interpol ;
1097
1098 }
1099
1100 if (n1 < 0 ) {
1101 cout << "Eos_bf_tabul::nbar_ent_p : nbar1<0 !" << endl ;
1102// abort() ;
1103// n1 = 0 ;
1104 }
1105 if (n2 < 0 ) {
1106 cout << "Eos_bf_tabul::nbar_ent_p : nbar2<0 !" << endl ;
1107// // abort() ;
1108// n2 = 0 ;
1109 }
1110 //cout << "nbar1 " << nbar1 << endl ;
1111 //cout << "nbar2 " << nbar2 << endl ;
1112
1113 if (pressure < 0 ) {
1114 cout << "Eos_bf_tabul::press_ent_p : pressure < 0 !" << endl ;
1115// abort() ;
1116// pressure = 0 ;
1117 }
1118
1119
1120 if (n1 >0.) {
1121
1122 K11 = mu1 / n1 - double(2) * alpha * ( 1. - xx) / ( n1 * n1 ) ; //OK
1123
1124// K11 = mu1 *n1 - double(2) * alpha * ( 1. - xx) ; // KNN----> KNN * Nn^2
1125
1126 }
1127 if (n2 >0.) {
1128
1129 K22 = mu2 / n2 - double(2) * alpha * ( 1. - xx) / ( n2 * n2 ) ; //OK
1130
1131// K22 = mu2 * n2 - double(2) * alpha * ( 1. - xx) ; // KNN----> KNN * Nn^2
1132
1133 }
1134 if ((n1 <= 0.) || (n2 <= 0.) ) {
1135
1136 K12 = 0. ;
1137
1138 }
1139 else {
1140
1141 K12 = double(2) * alpha * pow(1.-xx, 1.5)/ ( n1 * n2 ); //OK
1142// K12 = double(2) * alpha * pow(1.-xx, 1.5) ;// KNN----> KNN * Nn^2
1143
1144 }
1145
1146 nbar1.set(l,k,j,i) = n1 ;
1147 nbar2.set(l,k,j,i) = n2 ;
1148 press.set(l,k,j,i) = pressure ;
1149 ener.set(l,k,j,i) = mu1 * n1 + mu2 * n2 - pressure ;//ener_ent_p(xx1, xx2, xx) ;
1150 K_nn.set(l,k,j,i) = K11 ;
1151 K_np.set(l,k,j,i) = K12;
1152 K_pp.set(l,k,j,i) = K22 ;
1153
1154 //--------------------------------------------------------------------------------------------------------------------
1155
1156
1157
1158 //cout<< setprecision(16) ;
1159 //cout << "DEBUG-------> delta2 = " << xx << " test de la fonction pow : pow(1-delta2,1.5) = " << pow(1.-xx,1.5) << endl;
1160
1161 //Comparaison entre interpolation et calcul direct
1162 // cout<< setprecision(16) ;
1163 // cout << m_1 * exp (xx1 )* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " " << m_2 * exp (xx2 )* 2.99792458E+8 * 2.99792458E+8 * 1.66e-27 /(1.6021892E-13) << " :" << endl;
1164 // cout << "H1 = " << xx1 << " H2 = " << xx2 << " n1 = "<< nbar1(l,k,j,i) << " n2 = " << nbar2(l,k,j,i) << " P = " << press(l,k,j,i) << " E = " << ener(l,k,j,i) << endl;
1165
1166 //abort();
1167
1168 }
1169 }
1170 //abort();
1171
1172 }
1173 }
1174
1175}
1176
1177
1178// Baryon densities from enthalpies
1179//---------------------------------
1180
1181bool Eos_bf_tabul::nbar_ent_p(const double ent1, const double ent2,
1182 const double delta2, double& nbar1,
1183 double& nbar2) const {
1184
1185 bool one_fluid = false;
1186
1187 if (delta2 < delta_car_min) {
1188 cout << "Eos_bf_tabul::nbar_ent_p : delta2 < delta_car_min !" << endl ;
1189 abort() ;
1190 }
1191 if (delta2 > delta_car_max) {
1192 cout << "Eos_bf_tabul::nbar_ent_p : delta2 > delta_car_max !" << endl ;
1193 abort() ;
1194 }
1195 if (m_1 * exp(ent1) > ent1_max) {
1196 cout << "Eos_bf_tabul::nbar_ent_p : ent1 > ent1_max !" << endl ;
1197 abort() ;
1198 }
1199 if (m_2 *exp(ent2) > ent2_max) {
1200 cout << "Eos_bf_tabul::nbar_ent_p : ent2 > ent2_max !" << endl ;
1201 abort() ;
1202 }
1203
1204 if ( (exp(ent1) < 1.) && (exp(ent2) < 1.) ) {
1205 nbar1 = 0. ;
1206 nbar2 = 0. ;
1207 }
1208 else {
1209
1210 double mu1 = m_1 * exp(ent1);
1211 double mu2 = m_2 * exp(ent2);
1212
1213 double p_interpol ;
1214 double dpsdent1_interpol ;
1215 double dpsdent2_interpol ;
1216 double alpha ;
1217
1218 interpol_mixed_3d_new(m_1, m_2, *delta_car, *logent1, *logent2,
1221 *mu2_P, *n_p_P, *press_P,
1222 *mu1_N, *n_n_N, *press_N,
1223 *delta_car_n0, *mu1_n0, *mu2_n0,
1224 *delta_car_p0, *mu1_p0, *mu2_p0,
1225 delta2, mu1, mu2,
1226 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1227
1228 nbar1 = dpsdent1_interpol ;
1229 nbar2 = dpsdent2_interpol ;
1230
1231 }
1232
1233 if (nbar1 < 0 ) {
1234 cout << "Eos_bf_tabul::nbar_ent_p : nbar1<0 !" << endl ;
1235// abort() ;
1236 nbar1 = 0 ;
1237 }
1238 if (nbar2 < 0 ) {
1239 cout << "Eos_bf_tabul::nbar_ent_p : nbar2<0 !" << endl ;
1240// abort() ;
1241 nbar2 = 0 ;
1242 }
1243 //cout << "nbar1 " << nbar1 << endl ;
1244 //cout << "nbar2 " << nbar2 << endl ;
1245 one_fluid = ((nbar1 <= 0.)||(nbar2 <= 0.)) ;
1246
1247 return one_fluid;
1248}
1249
1250// One fluid sub-EOSs
1251//-------------------
1252
1253//PLUS APPELE ACTUELLEMENT!!!!!!!!!!!!!!!!!!!
1254double Eos_bf_tabul::nbar_ent_p1(const double ent1) const {
1255 //cout << " appel a nbar_ent_p1" << endl;
1256
1257 double pressN_interpol ;
1258 double n_n_N_interpol ;
1259 double mu1 = m_1 * exp(ent1);
1260 int i =0;
1261
1262 if (exp(ent1) < 1. ) {
1263 n_n_N_interpol = 0. ;
1264 }
1265 else {
1266 interpol_herm( *mu1_N, *press_N, *n_n_N,
1267 mu1, i, pressN_interpol , n_n_N_interpol) ;
1268 }
1269 return n_n_N_interpol ;
1270
1271}
1272
1273//PLUS APPELE ACTUELLEMENT!!!!!!!!!!!!!!!!!!!
1274 double Eos_bf_tabul::nbar_ent_p2(const double ent2) const {
1275 //cout << " appel a nbar_ent_p2" << endl;
1276 double pressP_interpol ;
1277 double n_p_P_interpol ;
1278 double mu2 = m_2 * exp(ent2);
1279 int i =0;
1280 if (exp(ent2) < 1. ) {
1281 n_p_P_interpol = 0. ;
1282 }
1283 else {
1284 interpol_herm( *mu2_P, *press_P, *n_p_P,
1285 mu2, i, pressP_interpol , n_p_P_interpol) ;
1286 }
1287 return n_p_P_interpol ;
1288
1289}
1290
1291 // Energy density from baryonic densities
1292//-----------------------------------------
1293double Eos_bf_tabul::ener_nbar_p(const double nbar1, const double nbar2,
1294 const double delta2) const{
1295
1296 c_est_pas_fait("Eos_bf_tabul::ener_nbar_p" ) ;
1297
1298 return nbar1 + nbar2 + delta2;
1299
1300 }
1301
1302// Pressure from baryonic densities
1303//----------------------------------
1304double Eos_bf_tabul::press_nbar_p(const double nbar1, const double nbar2,
1305 const double delta2) const {
1306
1307 c_est_pas_fait("Eos_bf_tabul::press_nbar_p" ) ;
1308
1309 return nbar1 + nbar2 + delta2;
1310
1311}
1312
1313
1314// Pressure from baryonic log-enthalpies
1315//--------------------------------------
1316 double Eos_bf_tabul::press_ent_p(const double ent1, const double ent2, const double delta2) const {
1317
1318 if (delta2 < delta_car_min) {
1319 cout << "Eos_bf_tabul::press_ent_p : delta2 < delta_car_min !" << endl ;
1320 abort() ;
1321 }
1322 if (delta2 > delta_car_max) {
1323 cout << "Eos_bf_tabul::press_ent_p : delta2 > delta_car_max !" << endl ;
1324 abort() ;
1325 }
1326 if (m_1 * exp(ent1) > ent1_max) {
1327 cout << "Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1328 abort() ;
1329 }
1330 if (m_2 * exp(ent2) > ent2_max) {
1331 cout << "Eos_bf_tabul::press_ent_p : ent2 > ent2_max !" << endl ;
1332 abort() ;
1333 }
1334
1335 double pressure ;
1336
1337 if ( (exp(ent1) < 1.) && (exp(ent2) < 1.)) {
1338 //abort();
1339 pressure = 0. ;
1340 }
1341
1342 else {
1343
1344 double mu1 = m_1 * exp(ent1);
1345 double mu2 = m_2 * exp(ent2);
1346
1347 double p_interpol ;
1348 double dpsdent1_interpol ;
1349 double dpsdent2_interpol ;
1350 double alpha ;
1351
1352 interpol_mixed_3d_new(m_1, m_2, *delta_car, *logent1, *logent2,
1355 *mu2_P, *n_p_P, *press_P,
1356 *mu1_N, *n_n_N, *press_N,
1357 *delta_car_n0, *mu1_n0, *mu2_n0,
1358 *delta_car_p0, *mu1_p0, *mu2_p0,
1359 delta2, mu1, mu2,
1360 p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1361
1362 pressure = p_interpol;
1363
1364 }
1365
1366 if (pressure < 0 ) {
1367 cout << "Eos_bf_tabul::press_ent_p : pressure < 0 !" << endl ;
1368// abort() ;
1369 pressure = 0 ;
1370 }
1371 return pressure ;
1372 }
1373
1374// PLUS APPELE MAINTEANTN
1375 double Eos_bf_tabul::press_ent_p1(const double ent1) const {
1376
1377 double pressure ;
1378 if (m_1 * exp(ent1) > ent1_max) {
1379 cout << "Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1380 abort() ;
1381 }
1382
1383 else {
1384
1385 double pressN_interpol ;
1386 double n_n_N_interpol ;
1387
1388 double mu1 = m_1 * exp(ent1);
1389 int i =0;
1390 interpol_herm( *mu1_N, *press_N, *n_n_N,
1391 mu1, i, pressN_interpol , n_n_N_interpol) ;
1392
1393 pressure = pressN_interpol ;
1394
1395 }
1396
1397 return pressure ;
1398 }
1399
1400 //PLUS APPELE MAINTENANT
1401double Eos_bf_tabul::press_ent_p2(const double ent2) const {
1402 double pressure ;
1403
1404 if (m_2*exp(ent2) > ent2_max) {
1405 cout << "Eos_bf_tabul::press_ent_p : ent1 > ent1_max !" << endl ;
1406 abort() ;
1407 }
1408 else {
1409
1410 double pressP_interpol ;
1411 double n_p_P_interpol ;
1412 double mu2 = m_2* exp(ent2);
1413 int i =0;
1414 interpol_herm( *mu2_P, *press_P, *n_p_P,
1415 mu2, i, pressP_interpol , n_p_P_interpol) ;
1416
1417 pressure = pressP_interpol ;
1418
1419 }
1420 return pressure ;
1421}
1422
1423
1424// Energy from baryonic log - enthalpies
1425//--------------------------------------
1426
1427 double Eos_bf_tabul::ener_ent_p(const double ent1, const double ent2,
1428 const double delta2) const {
1429 double energy= 0.;
1430
1431 if ( (exp(ent1) < 1.) && ( exp(ent2) < 1.)) {
1432 energy = 0. ;
1433 }
1434 else {
1435 /*Eos_bf_tabul::nbar_ent_p( ent1, ent2, delta2, nbar1, nbar2);
1436
1437 if ( (nbar1 > 0) && (nbar2 > 0) ) {
1438
1439 pressure = Eos_bf_tabul::press_ent_p(ent1,ent2,delta2);
1440 double mu_1 = m_1 * exp ( ent1 );
1441 double mu_2 = m_2 * exp ( ent2 );
1442 energy = nbar1 * mu_1 + nbar2 * mu_2 - pressure;
1443
1444 }
1445 else if ( (nbar1 <= 0) && (nbar2 > 0) ) {
1446
1447 //nbar1 = 0.;
1448 //nbar2 = Eos_bf_tabul::nbar_ent_p2(ent2);
1449 pressure = Eos_bf_tabul::press_ent_p(ent1,ent2,delta2);
1450 //pressure = Eos_bf_tabul::press_ent_p2(ent2);
1451 double mu_2 = m_2 * exp ( ent2 );
1452 energy = nbar2 * mu_2 - pressure;
1453 //cout << nbar1 << " " << pressure << " " << energy << endl ;
1454
1455 }
1456 else if ((nbar2 <= 0) && (nbar1 > 0)) {
1457 //nbar2 = 0.;
1458 //nbar1 = Eos_bf_tabul:: nbar_ent_p1(ent1);
1459 //pressure = Eos_bf_tabul::press_ent_p1(ent1);
1460 pressure = Eos_bf_tabul::press_ent_p(ent1,ent2,delta2);
1461 double mu_1 = m_1 * exp (ent1);
1462 energy = nbar1 * mu_1 - pressure;
1463 //cout << nbar1 << " " << pressure << " " << energy << endl ;
1464
1465 }
1466 else if ((nbar2 <= 0) && (nbar1 <= 0)) {
1467
1468 energy = 0.;*/
1469
1470 double mu1 = m_1 * exp(ent1);
1471 double mu2 = m_2 * exp(ent2);
1472
1473 double p_interpol ;
1474 double dpsdent1_interpol ;
1475 double dpsdent2_interpol ;
1476 double alpha ;
1477
1478 interpol_mixed_3d_new(m_1, m_2, *delta_car, *logent1, *logent2,
1481 *mu2_P, *n_p_P, *press_P,
1482 *mu1_N, *n_n_N, *press_N,
1483 *delta_car_n0, *mu1_n0, *mu2_n0,
1484 *delta_car_p0, *mu1_p0, *mu2_p0,
1485 delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1486
1487 energy = mu1 * dpsdent1_interpol + mu2 * dpsdent2_interpol - p_interpol ;
1488
1489 }
1490
1491 if (energy < 0 ) {
1492 cout << "Eos_bf_tabul::ener_ent_p : energy < 0 !" << endl ;
1493// abort() ;
1494 }
1495 return energy;
1496
1497}
1498
1499
1500// Alpha from baryonic log - enthalpies
1501//---------------------------------------
1502
1503double Eos_bf_tabul::alpha_ent_p(const double ent1, const double ent2,
1504 const double delta2) const {
1505
1506 if (delta2 < delta_car_min) {
1507 cout << "Eos_bf_tabul::alpha_ent_p : delta2 < delta_car_min !"
1508 << endl ;
1509 abort() ;
1510 }
1511 if (delta2 > delta_car_max) {
1512 cout << "Eos_bf_tabul::alpha_ent_p : delta2 > delta_car_max !"
1513 << endl ;
1514 abort() ;
1515 }
1516 if (m_1 * exp(ent1) > ent1_max) {
1517 cout << "Eos_bf_tabul::alpha_ent_p : ent1 > ent1_max !" << endl ;
1518 abort() ;
1519 }
1520 if (m_2 * exp(ent2) > ent2_max) {
1521 cout << "Eos_bf_tabul::alpha_ent_p : ent2 > ent2_max !" << endl ;
1522 abort() ;
1523 }
1524
1525 double alpha;
1526
1527 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ) {
1528 alpha = 0. ;
1529 }
1530
1531 else {
1532
1533 double mu1 = m_1 * exp(ent1);
1534 double mu2 = m_2 * exp(ent2);
1535
1536 double p_interpol ;
1537 double dpsdent1_interpol ;
1538 double dpsdent2_interpol ;
1539
1540 interpol_mixed_3d_new(m_1, m_2, *delta_car, *logent1, *logent2,
1543 *mu2_P, *n_p_P, *press_P,
1544 *mu1_N, *n_n_N, *press_N,
1545 *delta_car_n0, *mu1_n0, *mu2_n0,
1546 *delta_car_p0, *mu1_p0, *mu2_p0,
1547 delta2, mu1, mu2, p_interpol, dpsdent1_interpol, dpsdent2_interpol, alpha) ;
1548//cout << dpsdent1_interpol << " " << dpsdent2_interpol << " " << alpha <<endl;
1549 //cout << " alpha " << alpha << endl ;
1550 }
1551 //alpha = -alpha ; // NON NON NON LA CONVERSION EN -alpha est déjà faite dans l'interpolation
1552 //cout << ent1 << " " << ent2 << " " << alpha << endl;
1553 return alpha;
1554
1555
1556 }
1557
1558
1559
1560// Derivatives of energy
1561//----------------------
1562
1563double Eos_bf_tabul::get_K11(const double delta2, const double ent1, const double ent2) const {
1564
1565 double xx=0.;
1566 double mu_1 = m_1 * exp(ent1);
1567 double nbar1;
1568 double nbar2;
1569
1570 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1571 xx = 0. ;
1572 }
1573 else {
1574 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1, nbar2) ;
1575
1576 /*if (nbar1 <= 0.) {
1577 xx = 0. ;
1578 }
1579 else if ( (nbar1 > 0.) && ( nbar2 > 0. )) {
1580 //double mu_1 = m_1 * exp ( ent1 );
1581 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1582 xx = mu_1 / nbar1 - double(2) * alpha / ( nbar1 * nbar1 ) * ( 1. - delta2) ;
1583 }
1584 else if ( (nbar1 > 0.) && ( nbar2 <= 0. )) {
1585 //double mu_1 = m_1 * exp ( ent1 );
1586 //nbar1 = Eos_bf_tabul::nbar_ent_p1(ent1) ;
1587
1588 if (nbar1 !=0.) {
1589 xx = mu_1 / nbar1 ;
1590 }
1591 else { xx = 0.; }
1592 }*/
1593 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1594 if (nbar1 >0.) {
1595
1596 xx = mu_1 / nbar1 - double(2) * alpha * ( 1. - delta2) / ( nbar1 * nbar1 ) ;
1597 //cout << alpha << " " << mu_1 / nbar1 << " " << double(2) * alpha / ( nbar1 * nbar1 ) * ( 1. - delta2) << endl;
1598 }
1599 //cout << " test : " << " nbar = " << nbar1 << " " << mu_1 / nbar1 << " " << double(2) * alpha * ( 1. - delta2) / ( nbar1 * nbar1 ) << " " << " Knn= " <<
1600// xx << endl;
1601 //cout << "get_K11 : ent1 =" << ent1 << " ent2 = " << ent2 << " delta2 = " << delta2 << " nbar1 = " << nbar1 << " nbar2 = " << nbar2 << " alpha = " << alpha << " Knn = " << xx << endl;
1602 }
1603
1604
1605
1606 return xx;
1607
1608}
1609
1610double Eos_bf_tabul::get_K22( const double delta2, const double ent1, const double ent2) const {
1611
1612 double xx=0.;
1613 double mu_2 = m_2 * exp (ent2) ;
1614 double nbar1;
1615 double nbar2;
1616
1617 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1618 xx = 0. ;
1619 }
1620 else {
1621
1622 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1,nbar2) ;
1623
1624 /*if (nbar2 <= 0.) {
1625 xx = 0. ;
1626 }
1627 else if ( (nbar2 > 0.) && ( nbar1 > 0. )) {
1628 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1629 xx = mu_2 / nbar2 - double(2) * alpha / ( nbar2 * nbar2 ) * ( 1. - delta2);
1630 }
1631 else if ( (nbar2 > 0.) && ( nbar1 <= 0. )) {
1632
1633 //double mu_2 = m_2 * exp ( ent2);
1634 //nbar2 = Eos_bf_tabul::nbar_ent_p2(ent2) ;
1635 if (nbar2 !=0.) {
1636 xx = mu_2 / nbar2 ;
1637 }
1638 else { xx = 0. ;}
1639 }*/
1640 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1641 if (nbar2 >0.) {
1642
1643 xx = mu_2 / nbar2 - double(2) * alpha * ( 1. - delta2) / ( nbar2 * nbar2 ) ;
1644
1645 //cout << alpha << " " << mu_2 / nbar2 << " " << double(2) * alpha / ( nbar2 * nbar2 ) * ( 1. - delta2) << endl;
1646 }
1647// cout << " test : " << " nbar2 = " << nbar2 << " " << mu_2 / nbar2 << " " << double(2) * alpha * ( 1. - delta2) / ( nbar2 * nbar2 ) << " " << " Kpp= " <<
1648// xx << endl;
1649 //cout << "get_K22 : ent1 =" << ent1 << " ent2 = " << ent2 << " delta2 = " << delta2 << " nbar1 = " << nbar1 << " nbar2 = " << nbar2 << " alpha = " << alpha << " Kpp = " << xx << endl;
1650 }
1651
1652 return xx;
1653
1654}
1655
1656double Eos_bf_tabul::get_K12(const double delta2, const double ent1, const double ent2) const {
1657 double xx =0.;
1658 double nbar1;
1659 double nbar2;
1660
1661 if ((exp(ent1) <= 1.) && (exp(ent2) <= 1.) ){
1662 xx = 0. ;
1663 }
1664 else {
1665
1666 Eos_bf_tabul::nbar_ent_p(ent1,ent2, delta2, nbar1,nbar2) ;
1667
1668 double alpha = Eos_bf_tabul::alpha_ent_p(ent1,ent2,delta2) ;
1669 if ((nbar1 <= 0.) || (nbar2 <= 0.) ) {
1670 xx = 0. ;
1671 }
1672 else {
1673
1674
1675 xx = double(2) * alpha * pow(1.-delta2, 1.5)/ ( nbar1 * nbar2 );
1676 //cout << alpha << " " << double(2) * alpha / ( nbar1 * nbar2 ) * pow(1.-delta2, 1.5)<< endl;
1677 }
1678 //cout << "get_K12 : ent1 =" << ent1 << " ent2 = " << ent2 << " delta2 = " << delta2 << " nbar1 = " << nbar1 << " nbar2 = " << nbar2 << " alpha = " << alpha << " Knp = " << xx << endl;
1679 }
1680 // cout << ent1 << " " << ent2 << " " << xx << endl;
1681 return xx;
1682}
1683
1684
1685
1686// Conversion functions
1687// ---------------------
1688
1689//This function is necessary for "Et_rot", which needs an eos of type "Eos"
1690//But this eos is not used in the code, except for the construction of the star
1692
1693 Eos_poly* eos_simple = new Eos_poly(2.,0.016,1.008) ; // we can take whatever we want that makes sense as parameters
1694
1695 return eos_simple ;
1696}
1697
1698
1699
1700}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition cmp.C:323
int get_etat() const
Returns the logical state.
Definition cmp.h:899
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Class for a two-fluid (tabulated) equation of state.
Tbl * d3lpsdlent1dlent2ddelta_car
if necessary for the interpolation to find alpha (derivee seconde croisee) ie, if it's possible to ca...
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.
double ent1_max
Upper boundary of the log-enthalpy interval (fluid 1 = n)
virtual int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Definition eos_bf_file.C:99
void operator=(const Eos_bf_tabul &)
Assignment to another Eos_bf_tabul.
Tbl * logp
Table of .
double ent2_max
Upper boundary of the log-enthalpy interval (fluid 2 = p)
double ent2_min
Lower boundary of the log-enthalpy interval (fluid 2 = 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.
Tbl * d2lpsdlent2ddelta_car
Table of .
double delta_car_max
Upper boundary of the relative velocity interval --> 1 ?
Tbl * logent1
Table of where .
string authors
Authors.
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
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, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
Tbl * dlpsdlent2
Table of .
virtual ~Eos_bf_tabul()
Destructor.
string tablename
Name of the file containing the tabulated data (be careful, Eos_bifluid uses char*)
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...
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Tbl * delta_car
Table of .
double delta_car_min
Lower boundary of the relative velocity interval --> 0 ?
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
Tbl * dlpsdlent1
Table of .
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
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.
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.
virtual double get_K22(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy/(baryonic density 2) .
double ent1_min
Lower boundary of the log-enthalpy interval (fluid 1 = n)
void read_table()
Reads the file containing the table and initializes the arrays logent1, logent2, delta_car,...
virtual ostream & operator>>(ostream &) const
Operator >>
Tbl * d2lpsdlent2dlent2
Table of .
Tbl * logent2
Table of where .
Tbl * d2lpsdlent1dlent1
Table of .
Tbl * d2lpsdlent1ddelta_car
Table of .
virtual double press_ent_p1(const double ent1) const
Computes the pressure from the baryonic log-enthalpies asuming that only fluid 1 is present.
virtual double get_K12(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to .
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) .
virtual void sauve(FILE *) const
Save in a file.
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.
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
Tbl * d2lpsdlent1dlent2
Table of .
Eos_bf_tabul(const char *name_i, const char *table, const char *path, double mass1, double mass2)
Standard constructor.
Tbl * dlpsddelta_car
Table of
virtual double press_ent_p2(const double ent2) const
Computes the pressure from the baryonic log-enthalpies assuming that only fluid 2 is present.
2-fluids equation of state base class.
virtual void sauve(FILE *) const
Save in a file.
double m_1
Individual particle mass [unit: ].
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
double m_2
Individual particle mass [unit: ].
Polytropic equation of state (relativistic case).
Definition eos.h:757
Equation of state base class.
Definition eos.h:190
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.