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 char 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 
66 namespace Lorene {
67 
68 void interpol_linear(const Tbl& xtab, const Tbl& ytab,
69  double x, int& i, double& y) ;
70 
71 void interpol_herm(const Tbl& xtab, const Tbl& ytab, const Tbl& dytab,
72  double x, int& i, double& y, double& dy) ;
73 
74 void 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 
79 void 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 
88 void 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 
98 void 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 
112 void 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 */
116 void 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 // --------------------
128 Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* table,
129  const char* path, double mass1, double mass2) :
130 Eos_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 // ---------------------------------------
141 Eos_bf_tabul::Eos_bf_tabul(const char* name_i, const char* file_name,
142  double mass1, double mass2) :
143 Eos_bifluid(name_i, mass1, mass2)
144 {
145  tablename = file_name ;
146 
147  read_table() ;
148 }
149 
150 // Constructor from binary file
151 // ----------------------------
153 Eos_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
222  Eos_bifluid::operator=(eosi) ;
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 
262 void 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 
273 ostream& 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 
288 bool 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 
312 bool 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 
414  logent1->set_etat_qcq() ;
415  logent2->set_etat_qcq() ;
416  delta_car->set_etat_qcq() ;
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 
876 void 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 
1181 bool 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!!!!!!!!!!!!!!!!!!!
1254 double 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 //-----------------------------------------
1293 double 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 //----------------------------------
1304 double 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
1401 double 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 
1503 double 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 
1563 double 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 
1610 double 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 
1656 double 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 }
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_bifluid.C:245
string authors
Authors.
Definition: eos_bifluid.h:1544
double m_1
Individual particle mass [unit: ].
Definition: eos_bifluid.h:184
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Tbl * d3lpsdlent1dlent2ddelta_car
if necessary for the interpolation to find alpha (derivee seconde croisee) ie, if it&#39;s possible to ca...
Definition: eos_bifluid.h:1611
Tbl * logent2
Table of where .
Definition: eos_bifluid.h:1576
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
Tbl * d2lpsdlent2ddelta_car
Table of .
Definition: eos_bifluid.h:1606
Equation of state base class.
Definition: eos.h:190
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to. ...
virtual ostream & operator>>(ostream &) const
Operator >>
Definition: eos_bf_tabul.C:273
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Base class for coordinate mappings.
Definition: map.h:670
double ent1_min
Lower boundary of the log-enthalpy interval (fluid 1 = n)
Definition: eos_bifluid.h:1555
virtual double nbar_ent_p1(const double ent1) const
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present.
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 int identify() const
Returns a number to identify the sub-classe of Eos the object belongs to.
Definition: eos_bf_file.C:99
virtual double press_ent_p2(const double ent2) const
Computes the pressure from the baryonic log-enthalpies assuming that only fluid 2 is present...
Tbl * dlpsdlent1
Table of .
Definition: eos_bifluid.h:1585
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
Definition: eos_bifluid.C:232
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 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...
Tbl * logent1
Table of where .
Definition: eos_bifluid.h:1570
Tbl * d2lpsdlent1ddelta_car
Table of .
Definition: eos_bifluid.h:1603
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
double m_2
Individual particle mass [unit: ].
Definition: eos_bifluid.h:189
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Tbl * dlpsddelta_car
Table of .
Definition: eos_bifluid.h:1600
virtual double nbar_ent_p2(const double ent2) const
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present.
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 &#39;s, it computes both baryon densities, energy and pressure profi...
Definition: eos_bf_tabul.C:876
double ent2_max
Upper boundary of the log-enthalpy interval (fluid 2 = p)
Definition: eos_bifluid.h:1564
double delta_car_min
Lower boundary of the relative velocity interval –> 0 ?
Definition: eos_bifluid.h:1547
Tbl * d2lpsdlent1dlent1
Table of .
Definition: eos_bifluid.h:1594
virtual bool operator==(const Eos_bifluid &) const
Comparison operator (egality)
Definition: eos_bf_tabul.C:288
Polytropic equation of state (relativistic case).
Definition: eos.h:757
virtual double press_ent_p1(const double ent1) const
Computes the pressure from the baryonic log-enthalpies asuming that only fluid 1 is present...
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
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_max
Upper boundary of the log-enthalpy interval (fluid 1 = n)
Definition: eos_bifluid.h:1558
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_bf_tabul.C:262
virtual ~Eos_bf_tabul()
Destructor.
Definition: eos_bf_tabul.C:185
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void read_table()
Reads the file containing the table and initializes the arrays logent1, logent2, delta_car, logp, dlpsdlent1, dlpsdlent2, d2lpsdlent1dlent2, dlpsddelta_car, d2lpsdlent1ddelta_car, d2lpsdlent2ddelta_car, d3lpsdlent1dlent2ddelta_car.
Definition: eos_bf_tabul.C:323
virtual bool operator!=(const Eos_bifluid &) const
Comparison operator (difference)
Definition: eos_bf_tabul.C:312
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.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: cmp.C:323
Multi-domain grid.
Definition: grilles.h:273
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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_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. ...
void c_est_pas_fait(const char *)
Helpful function to say something is not implemented yet.
Tbl * dlpsdlent2
Table of .
Definition: eos_bifluid.h:1588
Tbl * delta_car
Table of .
Definition: eos_bifluid.h:1579
Tbl * d2lpsdlent2dlent2
Table of .
Definition: eos_bifluid.h:1596
Class for a two-fluid (tabulated) equation of state.
Definition: eos_bifluid.h:1534
virtual double get_K12(const double delta2, const double ent1, const double ent2) const
Computes the derivative of the energy with respect to .
Tbl * d2lpsdlent1dlent2
Table of .
Definition: eos_bifluid.h:1591
double delta_car_max
Upper boundary of the relative velocity interval –> 1 ?
Definition: eos_bifluid.h:1551
Basic array class.
Definition: tbl.h:161
void operator=(const Eos_bf_tabul &)
Assignment to another Eos_bf_tabul.
Definition: eos_bf_tabul.C:219
double ent2_min
Lower boundary of the log-enthalpy interval (fluid 2 = p)
Definition: eos_bifluid.h:1561
virtual Eos * trans2Eos() const
Makes a translation from Eos_bifluid to Eos .
string tablename
Name of the file containing the tabulated data (be careful, Eos_bifluid uses char*) ...
Definition: eos_bifluid.h:1541
Eos_bf_tabul(const char *name_i, const char *table, const char *path, double mass1, double mass2)
Standard constructor.
Definition: eos_bf_tabul.C:128
Tbl * logp
Table of .
Definition: eos_bifluid.h:1582
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.