LORENE
et_rot_bifluid.C
1 /*
2  * Methods for two fluids rotating relativistic stars.
3  *
4  * See the file et_rot_bifluid.h for documentation
5  *
6  */
7 
8 /*
9  * Copyright (c) 2001 Jerome Novak
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 et_rot_bifluid_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_bifluid.C,v 1.15 2015/06/11 13:50:19 j_novak Exp $" ;
31 
32 /*
33  * $Id: et_rot_bifluid.C,v 1.15 2015/06/11 13:50:19 j_novak Exp $
34  * $Log: et_rot_bifluid.C,v $
35  * Revision 1.15 2015/06/11 13:50:19 j_novak
36  * Minor corrections
37  *
38  * Revision 1.14 2015/06/10 14:39:17 a_sourie
39  * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
40  *
41  * Revision 1.13 2014/10/13 08:52:57 j_novak
42  * Lorene classes and functions now belong to the namespace Lorene.
43  *
44  * Revision 1.12 2004/03/25 10:29:04 j_novak
45  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
46  *
47  * Revision 1.11 2003/12/04 14:28:26 r_prix
48  * allow for the case of "slow-rot-style" EOS inversion, in which we need to adapt
49  * the inner domain to n_outer=0 instead of mu_outer=0 ...
50  * (this should only be used for comparison to analytic slow-rot solution!)
51  *
52  * Revision 1.10 2003/11/20 14:01:26 r_prix
53  * changed member names to better conform to Lorene coding standards:
54  * J_euler -> j_euler, EpS_euler -> enerps_euler, Delta_car -> delta_car
55  *
56  * Revision 1.9 2003/11/18 18:38:11 r_prix
57  * use of new member EpS_euler: matter sources in equilibrium() and global quantities
58  * no longer distinguish Newtonian/relativistic, as all terms should have the right limit...
59  *
60  * Revision 1.8 2003/11/17 13:49:43 r_prix
61  * - moved superluminal check into hydro_euler()
62  * - removed some warnings
63  *
64  * Revision 1.7 2003/11/13 12:07:57 r_prix
65  * *) changed xxx2 -> Delta_car
66  * *) added (non 2-fluid specific!) members sphph_euler J_euler
67  * *) more or less rewritten hydro_euler() to see if I understand it ;)
68  * - somewhat simplified and more adapted to the notation used in our notes/paper.
69  * - Main difference: u_euler is no longer used!!, the "output" instead
70  * consists of ener_euler, s_euler, sphph_euler and J_euler, which are
71  * the general 3+1 components for Tmunu.
72  *
73  * Revision 1.6 2003/09/17 08:27:50 j_novak
74  * New methods: mass_b1() and mass_b2().
75  *
76  * Revision 1.5 2002/10/18 08:42:58 j_novak
77  * Take into account the sign for uuu and uuu2
78  *
79  * Revision 1.4 2002/01/16 15:03:28 j_novak
80  * *** empty log message ***
81  *
82  * Revision 1.3 2002/01/11 14:09:34 j_novak
83  * Added newtonian version for 2-fluid stars
84  *
85  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
86  *
87  * All writing/reading to a binary file are now performed according to
88  * the big endian convention, whatever the system is big endian or
89  * small endian, thanks to the functions fwrite_be and fread_be
90  *
91  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
92  * LORENE
93  *
94  * Revision 1.3 2001/08/28 16:04:22 novak
95  * Use of new definition of relative velocity and new declarations for EOS
96  *
97  * Revision 1.2 2001/08/27 09:58:43 novak
98  * *** empty log message ***
99  *
100  * Revision 1.1 2001/06/22 15:39:17 novak
101  * Initial revision
102  *
103  *
104  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_bifluid.C,v 1.15 2015/06/11 13:50:19 j_novak Exp $
105  *
106  */
107 // Headers C
108 #include "math.h"
109 
110 // Headers Lorene
111 #include "et_rot_bifluid.h"
112 #include "nbr_spx.h"
113 #include "utilitaires.h"
114 #include "unites.h"
115 
116  //--------------//
117  // Constructors //
118  //--------------//
119 // Standard constructor
120 // --------------------
121 namespace Lorene {
122 Et_rot_bifluid::Et_rot_bifluid(Map& mpi, int nzet_i, bool relat, const Eos_bifluid& eos_i):
123  Etoile_rot(mpi, nzet_i, relat, *eos_i.trans2Eos()),
124  eos(eos_i),
125  ent2(mpi),
126  nbar2(mpi),
127  K_nn(mpi),
128  K_np(mpi),
129  K_pp(mpi),
130  sphph_euler(mpi),
131  j_euler(mpi, 1, CON, mp.get_bvect_cart()),
132  j_euler1 (mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments cinétiques de chaque fluide
133  j_euler2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments cinétiques de chaque fluide
134  j_euler11_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
135  j_euler12_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
136  j_euler21_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
137  j_euler22_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
138 
139  j_euler11_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
140  j_euler12_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
141  j_euler21_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
142  j_euler22_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
143 
144 
145  enerps_euler(mpi),
146  uuu2(mpi),
147  gam_euler2(mpi),
148  delta_car(mpi)
149 {
150  // All the matter quantities are initialized to zero :
151  nbar2 = 0 ;
152  ent2 = 0 ;
153  K_nn = 0 ;
154  K_np = 0 ;
155  K_pp = 0 ;
156  sphph_euler = 0;
157  j_euler = 0;
158  j_euler1 = 0 ;//Rajouté pour le calcul des moments cinétiques de chaque fluide
159  j_euler2 = 0;//Rajouté pour le calcul des moments cinétiques de chaque fluide
160  j_euler11_1 = 0 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
161  j_euler12_1 = 0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
162  j_euler21_1 = 0 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
163  j_euler22_1 = 0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
164  j_euler11_2 = 0 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
165  j_euler12_2 = 0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
166  j_euler21_2 = 0 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
167  j_euler22_2 = 0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
168 
169 
170  enerps_euler = 0;
172 
173  // Initialization to a static state :
174  omega2 = 0 ;
175  uuu2 = 0 ;
176  gam_euler2 = 1 ;
177  delta_car = 0 ;
178 
179  // Pointers of derived quantities initialized to zero :
180  set_der_0x0() ;
181 
182 }
183 
184 // Copy constructor
185 // ----------------
186 
188  Etoile_rot(et),
189  eos(et.eos),
190  ent2(et.ent2),
191  nbar2(et.nbar2),
192  K_nn(et.K_nn),
193  K_np(et.K_np),
194  K_pp(et.K_pp),
195  sphph_euler(et.sphph_euler),
196  j_euler(et.j_euler),
197  j_euler1(et.j_euler1),//Rajouté pour le calcul des moments cinétiques de chaque fluide
198  j_euler2(et.j_euler2),//Rajouté pour le calcul des moments cinétiques de chaque fluide
199  j_euler11_1(et.j_euler11_1),//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
200  j_euler12_1(et.j_euler12_1),//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
201  j_euler21_1(et.j_euler21_1),//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
202  j_euler22_1(et.j_euler22_1),//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)*
203  j_euler11_2(et.j_euler11_2),//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
204  j_euler12_2(et.j_euler12_2),//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
205  j_euler21_2(et.j_euler21_2),//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
206  j_euler22_2(et.j_euler22_2),//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
207  enerps_euler(et.enerps_euler),
208  uuu2(et.uuu2),
209  gam_euler2(et.gam_euler2),
210  delta_car(et.delta_car)
211 {
212  omega2 = et.omega2 ;
213 
214  // Pointers of derived quantities initialized to zero :
215  set_der_0x0() ;
216 }
217 
218 
219 // Constructor from a file
220 // ------------------------
221 Et_rot_bifluid::Et_rot_bifluid(Map& mpi, const Eos_bifluid& eos_i, FILE* fich):
222  Etoile_rot(mpi, *eos_i.trans2Eos(), fich),
223  eos(eos_i),
224  ent2(mpi),
225  nbar2(mpi),
226  K_nn(mpi),
227  K_np(mpi),
228  K_pp(mpi),
229  sphph_euler(mpi),
230  j_euler(mpi, 1, CON, mp.get_bvect_cart()),
231  j_euler1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments cinétiques de chaque fluide
232  j_euler2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments cinétiques de chaque fluide
233  j_euler11_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
234  j_euler12_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
235  j_euler21_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
236  j_euler22_1(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
237  j_euler11_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
238  j_euler12_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
239  j_euler21_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
240  j_euler22_2(mpi, 1, CON, mp.get_bvect_cart()), //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
241 
242 
243  enerps_euler(mpi),
244  uuu2(mpi),
245  gam_euler2(mpi),
246  delta_car(mpi)
247 {
248 
249  // Etoile parameters
250  // -----------------
251  // omega2 is read in the file:
252  fread_be(&omega2, sizeof(double), 1, fich) ;
253 
254 
255  // Read of the saved fields:
256  // ------------------------
257 
258  Tenseur ent2_file(mp, fich) ;
259  ent2 = ent2_file ;
260 
261  // All other fields are initialized to zero :
262  // ----------------------------------------
263  uuu2 = 0 ;
264  delta_car = 0 ;
265 
266  // Pointers of derived quantities initialized to zero
267  // --------------------------------------------------
268  set_der_0x0() ;
269 
270 }
271 
272  //------------//
273  // Destructor //
274  //------------//
275 
277 
278  del_deriv() ;
279 
280 }
281 
282  //----------------------------------//
283  // Management of derived quantities //
284  //----------------------------------//
285 
287 
289 
290  if (p_ray_eq2 != 0x0) delete p_ray_eq2 ;
291  if (p_ray_eq2_pis2 != 0x0) delete p_ray_eq2_pis2 ;
292  if (p_ray_eq2_pi != 0x0) delete p_ray_eq2_pi ;
293  if (p_ray_pole2 != 0x0) delete p_ray_pole2 ;
294  if (p_l_surf2 != 0x0) delete p_l_surf2 ;
295  if (p_xi_surf2 != 0x0) delete p_xi_surf2 ;
296  if (p_r_circ2 != 0x0) delete p_r_circ2 ;
297  if (p_area2 != 0x0) delete p_area2 ;
298  if (p_aplat2 != 0x0) delete p_aplat2 ;
299  if (p_mass_b1 != 0x0) delete p_mass_b1 ;
300  if (p_mass_b2 != 0x0) delete p_mass_b2 ;
301  if (p_angu_mom_1 != 0x0) delete p_angu_mom_1 ;//Rajouté pour le calcul des moments cinétiques de chaque fluide
302  if (p_angu_mom_2 != 0x0) delete p_angu_mom_2 ;//Rajouté pour le calcul des moments cinétiques de chaque fluide
303  if (p_angu_mom_1_part1_1 != 0x0) delete p_angu_mom_1_part1_1 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
304  if (p_angu_mom_2_part1_1 != 0x0) delete p_angu_mom_2_part1_1 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
305  if (p_angu_mom_1_part2_1 != 0x0) delete p_angu_mom_1_part2_1 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
306  if (p_angu_mom_2_part2_1 != 0x0) delete p_angu_mom_2_part2_1 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
307  if (p_angu_mom_1_part1_2 != 0x0) delete p_angu_mom_1_part1_2 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
308  if (p_angu_mom_2_part1_2 != 0x0) delete p_angu_mom_2_part1_2 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
309  if (p_angu_mom_1_part2_2 != 0x0) delete p_angu_mom_1_part2_2 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
310  if (p_angu_mom_2_part2_2 != 0x0) delete p_angu_mom_2_part2_2 ;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
311 
312 
313  set_der_0x0() ;
314 }
315 
316 
317 
318 
320 
322 
323  p_ray_eq2 = 0x0 ;
324  p_ray_eq2_pis2 = 0x0 ;
325  p_ray_eq2_pi = 0x0 ;
326  p_ray_pole2 = 0x0 ;
327  p_l_surf2 = 0x0 ;
328  p_xi_surf2 = 0x0 ;
329  p_r_circ2 = 0x0 ;
330  p_area2 = 0x0 ;
331  p_aplat2 = 0x0 ;
332  p_mass_b1 = 0x0;
333  p_mass_b2 = 0x0;
334  p_angu_mom_1 = 0x0;//Rajouté pour le calcul des moments cinétiques de chaque fluide
335  p_angu_mom_2 = 0x0;//Rajouté pour le calcul des moments cinétiques de chaque fluide
336  p_angu_mom_1_part1_1 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
337  p_angu_mom_2_part1_1 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
338  p_angu_mom_1_part2_1 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
339  p_angu_mom_2_part2_1 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
340  p_angu_mom_1_part1_2 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
341  p_angu_mom_2_part1_2 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
342  p_angu_mom_1_part2_2 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
343  p_angu_mom_2_part2_2 = 0x0;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
344 }
345 
347 
351  j_euler1.set_etat_nondef();//Rajouté pour le calcul des moments cinétiques de chaque fluide
352  j_euler2.set_etat_nondef();//Rajouté pour le calcul des moments cinétiques de chaque fluide
353  j_euler11_1.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
354  j_euler12_1.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
355  j_euler21_1.set_etat_nondef(); //Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
356  j_euler22_1.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
357  j_euler11_2.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
358  j_euler12_2.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
359  j_euler21_2.set_etat_nondef(); //Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
360  j_euler22_2.set_etat_nondef();//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
365  K_nn.set_etat_nondef(); //tester avec ou sans ces trois lignes !
368  del_deriv() ;
369 }
370 
371 // Assignment of the enthalpy field
372 // --------------------------------
373 
374 void Et_rot_bifluid::set_enthalpies(const Cmp& ent_i, const Cmp& ent2_i) {
375 
376  ent = ent_i ;
377  ent2 = ent2_i ;
378 
379  // Update of (nbar, ener, press) :
380  equation_of_state() ;
381 
382  // The derived quantities are obsolete:
383  del_deriv() ;
384 
385 }
386 
387 
388  //--------------//
389  // Assignment //
390  //--------------//
391 
392 // Assignment to another Et_rot_bifluid
393 // --------------------------------
395 
396  // Assignment of quantities common to all the derived classes of Etoile
397  Etoile_rot::operator=(et) ;
398 
399  assert( &(et.eos) == &eos ) ; // Same EOS
400  // Assignement of proper quantities of class Et_rot_bifluid
401  omega2 = et.omega2 ;
402 
403  ent2 = et.ent2 ;
404  nbar2 = et.nbar2 ;
405  K_nn = et.K_nn ;
406  K_np = et.K_np ;
407  K_pp = et.K_pp ;
409  j_euler = et.j_euler;
410  j_euler1 = et.j_euler1;//Rajouté pour le calcul des moments cinétiques de chaque fluide
411  j_euler2 = et.j_euler2;//Rajouté pour le calcul des moments cinétiques de chaque fluide
412  j_euler11_1 = et.j_euler11_1;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
413  j_euler12_1 = et.j_euler12_1;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
414  j_euler21_1 = et.j_euler21_1;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
415  j_euler22_1 = et.j_euler22_1;//Rajouté pour le calcul des moments d'inertie de chaque fluide (1ère version)
416  j_euler11_2 = et.j_euler11_2;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
417  j_euler12_2 = et.j_euler12_2;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
418  j_euler21_2 = et.j_euler21_2;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
419  j_euler22_2 = et.j_euler22_2;//Rajouté pour le calcul des moments d'inertie de chaque fluide (2ème version)
421  uuu2 = et.uuu2 ;
422  gam_euler2 = et.gam_euler2 ;
423  delta_car = et.delta_car ;
424 
425 
426  del_deriv() ; // Deletes all derived quantities
427 
428 }
429 
430  //--------------//
431  // Outputs //
432  //--------------//
433 
434 // Save in a file
435 // --------------
436 void Et_rot_bifluid::sauve(FILE* fich) const {
437 
438  Etoile_rot::sauve(fich) ;
439 
440  fwrite_be(&omega2, sizeof(double), 1, fich) ;
441 
442  ent2.sauve(fich) ;
443 
444 
445 }
446 
447 // Printing
448 // --------
449 
450 ostream& Et_rot_bifluid::operator>>(ostream& ost) const {
451 
452  using namespace Unites ;
453 
454  Etoile::operator>>(ost) ;
455 
456  ost << endl ;
457  ost << "Bifluid rotating star" << endl ;
458  ost << "-------------" << endl ;
459 
460  double freq = omega / (2.*M_PI) ;
461  ost << "Omega1 : " << omega * f_unit
462  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
463  ost << "Rotation period 1: " << 1000. / (freq * f_unit) << " ms"
464  << endl ;
465 
466  double freq2 = omega2 / (2.*M_PI) ;
467  ost << "Omega2 : " << omega2 * f_unit
468  << " rad/s f : " << freq2 * f_unit << " Hz" << endl ;
469  ost << "Rotation period 2: " << 1000. / (freq2 * f_unit) << " ms"
470  << endl ;
471 
472 //c----------------------------- angular momentum of fluid 1 and fluid 2 -------------------------------------------------
473  ost << "Total angular momentum J : "
474  << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
475  << endl ;
476  ost << "c J / (G M^2) : "
477  << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;
478 
479  double mom_iner = fabs(angu_mom() / omega2) ;
480  double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
481  ost << "Total moment of inertia I = J/Omega2 : "
482  << mom_iner_38si << " 10^38 kg m^2"
483  << endl ;
484 
485  ost << "Angular momentum J of fluid 1 : "
486  << angu_mom_1()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
487  << endl ;
488  ost << "Angular momentum J of fluid 2 : "
489  << angu_mom_2()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
490  << endl ;
491 
492 
493 //c----------------------------- momentum of inertia of fluid 1 and fluid 2 FIRST VERSION ----------------------------------------------
494  double mom_iner_1 = 0.;
495  double En_1=0.;
496  double mom_iner2_1 = 0. ;
497  double Ep_1=0.;
498  ost << "First way to define In, Ip, En and Ep : " << endl;
499  if (omega != __infinity) {
500  mom_iner_1 = fabs(angu_mom_1_part1_1() / omega) ;
501  double mom_iner_38si_1 = mom_iner_1 * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
502  ost << "Moment of inertia of fluid 1 : " << mom_iner_38si_1 << " 10^38 kg m^2"
503  << endl ;
504  if (omega != omega2) {
505  En_1 = angu_mom_1_part2_1() / ( (omega2 - omega) * mom_iner_1 ) ;
506  }
507  ost << "Coupling constante of fluid 1 (En) : " << En_1 << endl ;
508  }
509  if (omega2 != __infinity) {
510  mom_iner2_1 = fabs(angu_mom_2_part1_1() / omega2) ;
511  double mom_iner2_38si_1 = mom_iner2_1 * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
512  ost << "Moment of inertia of fluid 2 : " << mom_iner2_38si_1 << " 10^38 kg m^2"
513  << endl;
514 
515  if (omega != omega2) {
516  Ep_1 = angu_mom_2_part2_1() / ( (omega - omega2) * mom_iner2_1 ) ;
517  }
518  ost << "Coupling constante of fluid 2 (Ep) : " << Ep_1 << endl ;
519  }
520 
521  if (mom_iner2_1 * Ep_1 != 0. ) {
522  double X_1 = mom_iner_1 * En_1 / (mom_iner2_1 * Ep_1 );
523  ost << "Ratio X = (InEn)/(IpEp) : " << X_1 << endl ;
524  }
525 
526 //c----------------------------- momentum of inertia of fluid 1 and fluid 2 SECOND VERSION ----------------------------------------------
527  double mom_iner_2 = 0.;
528  double En_2=0.;
529  double mom_iner2_2 = 0. ;
530  double Ep_2=0.;
531  ost << "Second way to define In, Ip, En and Ep : " << endl;
532  if (omega != __infinity) {
533  mom_iner_2 = fabs(angu_mom_1_part1_2() / omega) ;
534  double mom_iner_38si_2 = mom_iner_2 * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
535  ost << "Moment of inertia of fluid 1 : " << mom_iner_38si_2 << " 10^38 kg m^2"
536  << endl ;
537  if (omega != omega2) {
538  En_2 = angu_mom_1_part2_2() / ( (omega2 - omega) * mom_iner_2 ) ;
539  }
540  ost << "Coupling constante of fluid 1 (En) : " << En_2 << endl ;
541  }
542  if (omega2 != __infinity) {
543  mom_iner2_2 = fabs(angu_mom_2_part1_2() / omega2) ;
544  double mom_iner2_38si_2 = mom_iner2_2 * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
545  ost << "Moment of inertia of fluid 2 : " << mom_iner2_38si_2 << " 10^38 kg m^2"
546  << endl ;
547 
548  if (omega != omega2) {
549  Ep_2 = angu_mom_2_part2_2() / ( (omega - omega2) * mom_iner2_2 ) ;
550  }
551  ost << "Coupling constante of fluid 2 (Ep) : " << Ep_2 << endl ;
552  }
553 
554  if (mom_iner2_2 * Ep_2 != 0. ) {
555  double X_2 = mom_iner_2 * En_2 / (mom_iner2_2 * Ep_2 );
556  ost << "Ratio X = (InEn)/(IpEp) : " << X_2 << endl ;
557  }
558 
559  ost.precision(16) ;
560  ost << "TEST ---- Jn + Jp = " <<angu_mom_2() + angu_mom_1() << " Jtot = " << angu_mom() << endl ;
561  ost << "TEST (vers. 1) ---- In + Ip = " <<mom_iner2_1 + mom_iner_1 << " I = " << mom_iner << endl ;
562  ost << "TEST (vers. 2) ---- In + Ip = " <<mom_iner2_2 + mom_iner_2 << " I = " << mom_iner << endl ;
563  ost << "TEST (vers. 1) ---- InOmegan + IpOmegap = " <<mom_iner2_1*omega2 + mom_iner_1*omega << " J = " << angu_mom() << endl ;
564  ost << "TEST (vers. 2) ---- InOmegan + IpOmegap = " <<mom_iner2_2*omega2 + mom_iner_2*omega << " J = " << angu_mom() << endl ;
565  ost << "TEST (vers. 1) ---- InOmegan + EnIn(Omegap-Omegan) = " <<mom_iner_1*omega + mom_iner_1*En_1*(omega2-omega) << " Jn = " << angu_mom_1() << endl ;
566  ost << "TEST (vers. 1) ---- IpOmegap + EpIp(Omegan-Omegap) = " <<mom_iner2_1*omega2 + mom_iner2_1*Ep_1*(omega-omega2) << " Jp = " << angu_mom_2() << endl ;
567  ost << "TEST (vers. 2) ---- InOmegan + EnIn(Omegap-Omegan) = " <<mom_iner_2*omega + mom_iner_2*En_2*(omega2-omega) << " Jn = " << angu_mom_1() << endl ;
568  ost << "TEST (vers. 2) ---- IpOmegap + EpIp(Omegan-Omegap) = " <<mom_iner2_2*omega2 + mom_iner2_2*Ep_2*(omega-omega2) << " Jp = " << angu_mom_2() << endl ;
569 
570 
571 // pour verifier
572 
573 //ost << "Verification : " << "Jn = "<< angu_mom_1() << " In * ... " << mom_iner * omega + mom_iner * En * (omega2 - omega) << endl ;
574 //ost << "Verification : " << "Jp = "<< angu_mom_2() << " Ip * ... " << mom_iner2* omega2 + mom_iner2 * Ep * (omega - omega2) << endl ;
575 
576 //c------------------------------------------------------------------------------------------------------------------------------------
577 
578  double nphi_c = nphi()(0, 0, 0, 0) ;
579  if ( (omega==0) && (nphi_c==0) ) {
580  ost << "Central N^phi : " << nphi_c << endl ;
581  }
582  else{
583  ost << "Central N^phi/Omega : " << nphi_c / omega << endl ;
584  }
585  if (omega2!=0)
586  ost << "Central N^phi/Omega2 : " << nphi_c / omega2 << endl ;
587 
588  ost << "Error on the virial identity GRV2 : " << endl ;
589  ost << "GRV2 = " << grv2() << endl ;
590  ost << "Error on the virial identity GRV3 : " << endl ;
591  double xgrv3 = grv3(&ost) ;
592  ost << "GRV3 = " << xgrv3 << endl ;
593 
594  ost << "Circumferential equatorial radius R_circ : "
595  << r_circ()/km << " km" << endl ;
596  ost << "Mean radius R_mean : "
597  << mean_radius()/km << " km" << endl ;
598  ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
599  << endl ;
600  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
601  ost << "Circumferential equatorial radius R_circ2 : "
602  << r_circ2()/km << " km" << endl ;
603  ost << "Mean radius R_mean2 : "
604  << mean_radius2()/km << " km" << endl ;
605  ost << "Coordinate equatorial radius r_eq2 : " << ray_eq2()/km << " km"
606  << endl ;
607  ost << "Flattening r_pole2/r_eq2 : " << aplat2() << endl ;
608 
609  int lsurf = nzet - 1;
610  int nt = mp.get_mg()->get_nt(lsurf) ;
611  int nr = mp.get_mg()->get_nr(lsurf) ;
612  ost << "Equatorial value of the velocity U: "
613  << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
614  ost << "Equatorial value of the velocity U2: "
615  << uuu2()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
616  ost << "Redshift at the equator (forward) : " << z_eqf() << endl ;
617  ost << "Redshift at the equator (backward): " << z_eqb() << endl ;
618  ost << "Redshift at the pole : " << z_pole() << endl ;
619 
620 
621  ost << "Central value of log(N) : "
622  << logn()(0, 0, 0, 0) << endl ;
623 
624  ost << "Central value of dzeta=log(AN) : "
625  << dzeta()(0, 0, 0, 0) << endl ;
626 
627  ost << "Central value of B^2 : " << b_car()(0,0,0,0) << endl ;
628 
629  Tbl diff_a_b = diffrel( a_car(), b_car() ) ;
630  ost <<
631  "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
632  << endl ;
633  for (int l=0; l<diff_a_b.get_taille(); l++) {
634  ost << diff_a_b(l) << " " ;
635  }
636  ost << endl;
637  ost << "Quadrupole moment : " << mom_quad() << endl ;
638  double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.)) / double(1.e38) ) ;
639  ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
640  << endl ;
641  ost << "Old quadrupole moment : " << mom_quad_old() << endl ;
642  ost << "Coefficient b : " << mom_quad_Bo() / pow(mass_g(), 2.) << endl ;
643  ost << "q = c^4 Q / (G^2 M^3) : "
644  << mom_quad() / ( ggrav * ggrav * pow(mass_g(), 3.) ) << endl ;
645  ost << "j = c J / (G M^2) : "
646  << angu_mom()/( ggrav * pow(mass_g(), 2.) ) << endl ;
647 
648  // Pour vérifier la loi de Laarakkers et Poisson, 1999
649  ost << " Mb (Msol) Mg (Msol) omegan=omega (rad/s) q q_old b j j2 " << endl;
650  ost << mass_b() / msol << " "
651  << mass_g() / msol << " "
652  << omega * f_unit << " "
653  << mom_quad() / ( ggrav * ggrav * pow(mass_g(), 3.) ) << " "
654  << mom_quad_old() / ( ggrav * ggrav * pow(mass_g(), 3.) ) << " "
655  << mom_quad_Bo() / pow(mass_g(), 2.) << " "
656  << angu_mom()/( ggrav * pow(mass_g(), 2.) ) << " "
657  << pow(angu_mom()/( ggrav * pow(mass_g(), 2.) ) , 2.) << " "
658  << endl ;
659  ost << endl ;
660 
661  return ost ;
662 
663 }
664 
665 
666 void Et_rot_bifluid::partial_display(ostream& ost) const {
667 
668  using namespace Unites ;
669 
670  double freq = omega / (2.*M_PI) ;
671  ost << "Omega : " << omega * f_unit
672  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
673  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
674  << endl ;
675  ost << endl << "Central enthalpy : " << ent()(0,0,0,0) << " c^2" << endl ;
676  ost << "Central proper baryon density : " << nbar()(0,0,0,0)
677  << " x 0.1 fm^-3" << endl ;
678  double freq2 = omega2 / (2.*M_PI) ;
679  ost << "Omega2 : " << omega2 * f_unit
680  << " rad/s f : " << freq2 * f_unit << " Hz" << endl ;
681  ost << "Rotation period 2: " << 1000. / (freq2 * f_unit) << " ms"
682  << endl ;
683  ost << endl << "Central enthalpy 2: " << ent2()(0,0,0,0) << " c^2" << endl ;
684  ost << "Central proper baryon density 2: " << nbar2()(0,0,0,0)
685  << " x 0.1 fm^-3" << endl ;
686  ost << "Central proper energy density : " << ener()(0,0,0,0)
687  << " rho_nuc c^2" << endl ;
688  ost << "Central pressure : " << press()(0,0,0,0)
689  << " rho_nuc c^2" << endl ;
690 
691  ost << "Central value of log(N) : "
692  << logn()(0, 0, 0, 0) << endl ;
693  ost << "Central lapse N : " << nnn()(0,0,0,0) << endl ;
694  ost << "Central value of dzeta=log(AN) : "
695  << dzeta()(0, 0, 0, 0) << endl ;
696  ost << "Central value of A^2 : " << a_car()(0,0,0,0) << endl ;
697  ost << "Central value of B^2 : " << b_car()(0,0,0,0) << endl ;
698 
699  double nphi_c = nphi()(0, 0, 0, 0) ;
700  if ( (omega==0) && (nphi_c==0) ) {
701  ost << "Central N^phi : " << nphi_c << endl ;
702  }
703  else{
704  ost << "Central N^phi/Omega : " << nphi_c / omega << endl ;
705  }
706 
707 
708  int lsurf = nzet - 1;
709  int nt = mp.get_mg()->get_nt(lsurf) ;
710  int nr = mp.get_mg()->get_nr(lsurf) ;
711  ost << "Equatorial value of the velocity U: "
712  << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
713 
714  ost << "Equatorial value of the velocity U2: "
715  << uuu2()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
716 
717  ost << endl
718  << "Coordinate equatorial radius r_eq = "
719  << ray_eq()/km << " km" << endl ;
720  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
721  ost << endl
722  << "Coordinate equatorial radius r_eq2 = "
723  << ray_eq2()/km << " km" << endl ;
724  ost << "Flattening r_pole2/r_eq2 : " << aplat2() << endl ;
725 
726 }
727 
728 //
729 // Equation of state
730 //
731 
733 
734  Cmp ent_eos = ent() ;
735  Cmp ent2_eos = ent2() ;
736  Tenseur rel_vel(delta_car) ;
737 
738  if (nzet > 1) {
739  // Slight rescale of the enthalpy field in case of 2 domains inside the
740  // star
741 
742  if (nzet > 2) {
743  cout << "Et_rot_bifluid::equation_of_state: not ready yet for nzet > 2 !" << endl ;
744  }
745 
746  double epsilon = 1.e-12 ;
747 
748  const Mg3d* mg = mp.get_mg() ;
749  int nz = mg->get_nzone() ;
750 
751  Mtbl xi(mg) ;
752  xi.set_etat_qcq() ;
753  for (int l=0; l<nz; l++) {
754  xi.t[l]->set_etat_qcq() ;
755  for (int k=0; k<mg->get_np(l); k++) {
756  for (int j=0; j<mg->get_nt(l); j++) {
757  for (int i=0; i<mg->get_nr(l); i++) {
758  xi.set(l,k,j,i) =
759  mg->get_grille3d(l)->x[i] ;
760  }
761  }
762  }
763 
764  }
765 
766  Cmp fact_ent(mp) ;
767  fact_ent.allocate_all() ;
768 
769  fact_ent.set(0) = 1 + epsilon * xi(0) * xi(0) ;
770  fact_ent.set(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
771 
772  for (int l=nzet; l<nz; l++) {
773  fact_ent.set(l) = 1 ;
774  }
775 
776  ent_eos = fact_ent * ent_eos ;
777  ent2_eos = fact_ent * ent2_eos ;
778  ent_eos.std_base_scal() ;
779  ent2_eos.std_base_scal() ;
780 
781  } // if nzet > 1
782 
783 
784  // Call to the EOS
785  nbar.set_etat_qcq() ;
786  nbar2.set_etat_qcq() ;
787  ener.set_etat_qcq() ;
788  press.set_etat_qcq() ;
789  K_nn.set_etat_qcq() ;
790  K_np.set_etat_qcq() ;
791  K_pp.set_etat_qcq() ;
792 
793 
794  const Eos_bf_tabul* eos_t = dynamic_cast<const Eos_bf_tabul*>(&eos) ;
795 
796  if (eos_t != 0x0) { // if the equation of state is tabulated
797  eos_t->calcule_interpol(ent_eos, ent2_eos, rel_vel(), nbar.set(), nbar2.set(),
798  ener.set(), press.set(), K_nn.set(), K_np.set(), K_pp.set(),
799  nzet) ;
800  }
801  else { // if the equation of state is not tabulated
802 
803  eos.calcule_tout(ent_eos, ent2_eos, rel_vel(), nbar.set(), nbar2.set(),
804  ener.set(), press.set(), nzet) ;
805  // est-ce que Knn, Knp, et Kpp bien calculée ?
806 
807  }
808  // Set the bases for spectral expansion
809  nbar.set_std_base() ;
810  nbar2.set_std_base() ;
811  ener.set_std_base() ;
812  press.set_std_base() ;
813  K_pp.set_std_base() ;
814  K_nn.set_std_base() ;
815  K_np.set_std_base() ;
816 
817  // The derived quantities are obsolete
818  del_deriv() ;
819 
820 }
821 
822 //
823 // Computation of hydro quantities
824 //
825 
827 
828  const Mg3d* mg = mp.get_mg();
829  int nz = mg->get_nzone() ;
830  int nzm1 = nz - 1 ;
831 
832  // RP: I prefer to use the 3-vector j_euler instead of u_euler
833  // for better physical "encapsulation"
834  // (i.e. --> use same form of Poisson-equations for all etoile sub-classes!)
835  u_euler.set_etat_nondef(); // make sure it's not used
836 
837  // (Norm of) Euler-velocity of the first fluid
838  //------------------------------
839  uuu.set_etat_qcq();
840 
841  uuu.set() = bbb() * (omega - nphi() ) / nnn();
842  uuu.annule(nzm1) ;
843 
844  // gosh, we have to exclude the thing being zero here... :(
845  if( uuu.get_etat() != ETATZERO )
846  {
847  (uuu.set()).std_base_scal() ;
848  (uuu.set()).mult_rsint();
849  }
850  uuu.set_std_base();
851 
852 
853  // (Norm of) Euler-velocity of the second fluid
854  //----------------------------------------
855  uuu2.set_etat_qcq();
856 
857  uuu2.set() = bbb() * (omega2 - nphi() ) / nnn();
858  uuu2.annule(nzm1) ;
859 
860  if( uuu2.get_etat() != ETATZERO )
861  {
862  (uuu2.set()).std_base_scal();
863  (uuu2.set()).mult_rsint();
864  }
865  uuu2.set_std_base();
866 
867  // Sanity check:
868  // Is one of the new velocities larger than c in the equatorial plane ?
869  //----------------------------------------
870 
871  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
872  int l_b = nzet - 1 ;
873  int j_b = mg->get_nt(l_b) - 1 ;
874 
875  bool superlum = false ;
876  if (relativistic) {
877  for (int l=0; l<nzet; l++) {
878  for (int i=0; i<mg->get_nr(l); i++) {
879 
880  double u1 = uuu()(l, 0, j_b, i) ;
881  double u2 = uuu2()(l, 0, j_b, i) ;
882  if ((u1 >= 1.) || (u2>=1.)) { // superluminal velocity
883  superlum = true ;
884  cout << "U > c for l, i : " << l << " " << i
885  << " U1 = " << u1 << endl ;
886  cout << " U2 = " << u2 << endl ;
887  }
888  }
889  }
890  if ( superlum ) {
891  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
892  abort() ;
893  }
894  }
895 
896 
897  Tenseur uuu_car = uuu * uuu;
898  Tenseur uuu2_car = uuu2 * uuu2;
899 
900 
901  // Lorentz factors
902  // --------------
903  Tenseur gam_car = 1.0 / (1.0 - unsurc2 * uuu_car) ;
904  gam_euler = sqrt(gam_car) ;
905  gam_euler.set_std_base() ; // sets the standard spectral bases for a scalar field
906 
907 
908  Tenseur gam2_car = 1.0 / (1.0 - unsurc2 * uuu2_car) ;
909  gam_euler2 = sqrt(gam2_car) ;
911 
912  // Update of "relative velocity" squared: $\Delta^2$
913  // ---------------------------
914 
915  delta_car = (uuu - uuu2)*(uuu - uuu2) / ( (1 - unsurc2* uuu*uuu2) *(1 - unsurc2* uuu*uuu2 ) ) ;
917 
918 //------------------------------------------------------------------------------------------------------------------------------
919  Tenseur Ann(ent) ;
920  Tenseur Anp(ent) ;
921  Tenseur App(ent) ;
922 
923  if (eos.identify() == 3) {
924 
925 // Ann = gam_car() * nbar() * nbar() * eos.get_Knn_ent(ent(), ent2(), delta_car(), nzet);
926 // Anp = gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
927 // App = gam2_car() * nbar2() * nbar2() * eos.get_Kpp_ent(ent(), ent2(), delta_car(), nzet);
928 
929  Ann = gam_car() * nbar() * nbar() * K_nn() ; // OK
930  Anp = gam_euler() * gam_euler2() * nbar() * nbar2() * K_np(); //OK
931  App = gam2_car() * nbar2() * nbar2() * K_pp(); //OK
932 
933 // Ann = gam_car() * K_nn() ; // Knn->Knn * nn * nn
934 // Anp = gam_euler() * gam_euler2() * K_np(); // Knn->Knn * nn * nn
935 // App = gam2_car() * K_pp(); // Knn->Knn * nn * nn
936  }
937 
938  else {
939 
940  Ann = gam_car() * nbar() * nbar() * eos.get_Knn(nbar(), nbar2(), delta_car(), nzet);
941  Anp = gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
942  App = gam2_car() * nbar2() * nbar2() * eos.get_Kpp(nbar(), nbar2(), delta_car(), nzet);
943 
944 
945  }
946 
947 
948  // Energy density E with respect to the Eulerian observer
949  //------------------------------------
950  // use of ener_euler is deprecated, because it's useless in Newtonian limit!
951  ener_euler.set_etat_nondef(); // make sure, it's not used
952 
953  Tenseur E_euler(mp);
954  E_euler = Ann + 2. * Anp + App - press ;
955  E_euler.set_std_base() ;
956 
957 
958  // S^phi_phi component of stress-tensor S^i_j
959  //------------------------------------
960  sphph_euler = press() + Ann() * uuu_car() + 2. * Anp() * uuu() * uuu2() + App() * uuu2_car();
962 
963 
964  // Trace of the stress tensor with respect to the Eulerian observer
965  //------------------------------------
966  s_euler = 2. * press() + sphph_euler();
967  s_euler.set_std_base() ;
968 
969  // The combination enerps_euler := (E + S_i^i) which has Newtonian limit -> rho
970  if (relativistic)
971  enerps_euler = E_euler + s_euler;
972  else
973  enerps_euler = eos.get_m1() * nbar() + eos.get_m2() * nbar2();
974 
975 
976  // the (flat-space) angular-momentum 3-vector j_euler^i
977  //-----------------------------------
978  Tenseur Jph(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
979  Jph = Ann*uuu + Anp*(uuu + uuu2) + App*uuu2 ;
980 
982 
983  j_euler.set(0) = 0; // r tetrad component
984  j_euler.set(1) = 0; // theta tetrad component
985  j_euler.set(2) = Jph()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
988 
989  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
990  j_euler.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
991 
992  if( (j_euler(0).get_etat() == ETATZERO)&&(j_euler(1).get_etat() == ETATZERO)&&(j_euler(2).get_etat()==ETATZERO))
993  j_euler = 0;
994 
995  // Rajouté pour le calcul des moments cinétiques de chaque fluide
996  // the (flat-space) angular-momentum 3-vector j_euler^i for fluid 1
997  //-----------------------------------
998  Tenseur Jph1(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
999  Jph1 = Ann*uuu + Anp*uuu2 ;
1000 
1002 
1003  j_euler1.set(0) = 0; // r tetrad component
1004  j_euler1.set(1) = 0; // theta tetrad component
1005  j_euler1.set(2) = Jph1()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1008 
1009  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1010  j_euler1.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1011 
1012  if( (j_euler1(0).get_etat() == ETATZERO)&&(j_euler1(1).get_etat() == ETATZERO)&&(j_euler1(2).get_etat()==ETATZERO))
1013  j_euler1 = 0;
1014 
1015  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1016  // the (flat-space) angular-momentum 3-vector j_euler^i for fluid 2
1017  //-----------------------------------
1018  Tenseur Jph2(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1019  Jph2 = Anp*uuu + App*uuu2 ;
1020 
1022 
1023  j_euler2.set(0) = 0; // r tetrad component
1024  j_euler2.set(1) = 0; // theta tetrad component
1025  j_euler2.set(2) = Jph2()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1028 
1029  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1030  j_euler2.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1031 
1032  if( (j_euler2(0).get_etat() == ETATZERO)&&(j_euler2(1).get_etat() == ETATZERO)&&(j_euler2(2).get_etat()==ETATZERO))
1033  j_euler2 = 0;
1034 
1035 
1036 
1037  // Computation of elements used for the first version of the definition of In, Ip, En and Ep...
1038  // the (flat-space) first part angular-momentum 3-vector j_euler^i for fluid 1
1039  //-----------------------------------
1040  Tenseur Jph11_1(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1041 
1042  Tenseur A11_1(ent) ;
1043 
1044  if (eos.identify() == 3) {
1045 
1046 // A11_1 = gam_car() * nbar() * nbar() * eos.get_Knn_ent(ent(), ent2(), delta_car(), nzet) +
1047 // gam_car() * nbar() * nbar2() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1048 
1049  A11_1 = gam_car() * nbar() * nbar() * K_nn() + gam_car() * nbar() * nbar2() * K_np() * pow (1.-delta_car(), -0.5) ; // OK
1050 
1051 // A11_1 = gam_car() * K_nn() + gam_car() * K_np() * pow (1.-delta_car(), -0.5) ; // Knn->Knn * nn * nn
1052 
1053  //A11_1 = gam_car() * nbar() * eos.get_m1() * exp(ent()); // erreur légèrement augmentée par cette méthode
1054  }
1055 
1056  else {
1057 
1058  A11_1 = gam_car() * nbar() * nbar() * eos.get_Knn(nbar(), nbar2(), delta_car(), nzet) +
1059  gam_car() * nbar() * nbar2() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1060  }
1061 
1062  Jph11_1 = A11_1 * uuu ;
1063 
1065 
1066  j_euler11_1.set(0) = 0; // r tetrad component
1067  j_euler11_1.set(1) = 0; // theta tetrad component
1068  j_euler11_1.set(2) = Jph11_1()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1071 
1072  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1073  j_euler11_1.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1074 
1075  if( (j_euler11_1(0).get_etat() == ETATZERO)&&(j_euler11_1(1).get_etat() == ETATZERO)&&(j_euler11_1(2).get_etat()==ETATZERO))
1076  j_euler11_1 = 0;
1077 
1078  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1079  // the (flat-space) second part angular-momentum 3-vector j_euler^i for fluid 1
1080  //-----------------------------------
1081  Tenseur Jph12_1(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1082 
1083  Tenseur A12_1(ent) ;
1084  Tenseur B12_1(ent) ;
1085 
1086  if (eos.identify() == 3) {
1087 // A12_1 = - gam_car() * nbar() * nbar2() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1088 // B12_1 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1089 
1090  A12_1 = - gam_car() * nbar() * nbar2() * K_np() * pow (1.-delta_car(), -0.5) ; //OK
1091  B12_1 = + gam_euler() * gam_euler2() * nbar() * nbar2() * K_np(); //OK
1092 
1093 // A12_1 = - gam_car() * K_np() * pow (1.-delta_car(), -0.5) ; // Knn->Knn * nn * nn
1094 // B12_1 = + gam_euler() * gam_euler2() * K_np(); // Knn->Knn * nn * nn
1095  }
1096 
1097  else {
1098  A12_1 = - gam_car() * nbar() * nbar2() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1099  B12_1 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
1100  }
1101 
1102  Jph12_1 = A12_1 * uuu + B12_1 * uuu2;
1103 
1105 
1106  j_euler12_1.set(0) = 0; // r tetrad component
1107  j_euler12_1.set(1) = 0; // theta tetrad component
1108  j_euler12_1.set(2) = Jph12_1()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1111 
1112  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1113  j_euler12_1.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1114 
1115  if( (j_euler12_1(0).get_etat() == ETATZERO)&&(j_euler12_1(1).get_etat() == ETATZERO)&&(j_euler12_1(2).get_etat()==ETATZERO))
1116  j_euler12_1 = 0;
1117 
1118  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1119  // the (flat-space) first part angular-momentum 3-vector j_euler^i for fluid 2
1120  //-----------------------------------
1121  Tenseur Jph21_1(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1122 
1123  Tenseur A21_1(ent) ;
1124 
1125  if (eos.identify() == 3) {
1126 
1127 // A21_1 = gam2_car() * nbar2() * nbar2() * eos.get_Kpp_ent(ent(), ent2(), delta_car(), nzet) +
1128 // gam2_car() * nbar2() * nbar() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1129 
1130  A21_1 = gam2_car() * nbar2() * nbar2() * K_pp() + gam2_car() * nbar2() * nbar() * K_np() * pow (1.-delta_car(), -0.5) ; //OK
1131 
1132 // A21_1 = gam2_car() * K_pp() + gam2_car() * K_np() * pow (1.-delta_car(), -0.5) ; // Knn->Knn * nn * nn
1133 
1134 
1135 // A21_1 = gam2_car() * nbar2() * eos.get_m2() * exp(ent2()); // erreures légèrement augmentées avec cette méthode
1136  }
1137 
1138  else {
1139 
1140  A21_1 = gam2_car() * nbar2() * nbar2() * eos.get_Kpp(nbar(), nbar2(), delta_car(), nzet) +
1141  gam2_car() * nbar2() * nbar() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1142  }
1143 
1144  Jph21_1 = A21_1 * uuu2 ;
1145 
1146  j_euler21_1.set_etat_qcq();
1147 
1148  j_euler21_1.set(0) = 0; // r tetrad component
1149  j_euler21_1.set(1) = 0; // theta tetrad component
1150  j_euler21_1.set(2) = Jph21_1()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1151  j_euler21_1.set_triad (mp.get_bvect_spher());
1152  j_euler21_1.set_std_base();
1153 
1154  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1155  j_euler21_1.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1156 
1157  if( (j_euler21_1(0).get_etat() == ETATZERO)&&(j_euler21_1(1).get_etat() == ETATZERO)&&(j_euler21_1(2).get_etat()==ETATZERO))
1158  j_euler21_1 = 0;
1159 
1160  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1161  // the (flat-space) second part angular-momentum 3-vector j_euler^i for fluid 1
1162  //-----------------------------------
1163  Tenseur Jph22_1(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1164 
1165  Tenseur A22_1(ent) ;
1166  Tenseur B22_1(ent) ;
1167 
1168  if (eos.identify() == 3) {
1169 // A22_1 = - gam2_car() * nbar2() * nbar() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1170 // B22_1 = + gam_euler() * gam_euler2() * nbar()*nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1171 
1172  A22_1 = - gam2_car() * nbar2() * nbar() * K_np() * pow (1.-delta_car(), -0.5) ; //OK
1173  B22_1 = + gam_euler() * gam_euler2() * nbar()*nbar2() * K_np(); //OK
1174 
1175 // A22_1 = - gam2_car() * K_np() * pow (1.-delta_car(), -0.5) ; // Knn->Knn * nn * nn
1176 // B22_1 = + gam_euler() * gam_euler2() * K_np(); // Knn->Knn * nn * nn
1177  }
1178 
1179  else {
1180  // some auxiliary EOS variables
1181  A22_1 = - gam2_car() * nbar2() * nbar() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) * pow (1.-delta_car(), -0.5) ;
1182  B22_1 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
1183  }
1184 
1185  Jph22_1 = A22_1 * uuu2 + B22_1 * uuu;
1186 
1187  j_euler22_1.set_etat_qcq();
1188 
1189  j_euler22_1.set(0) = 0; // r tetrad component
1190  j_euler22_1.set(1) = 0; // theta tetrad component
1191  j_euler22_1.set(2) = Jph22_1()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1192  j_euler22_1.set_triad (mp.get_bvect_spher());
1193  j_euler22_1.set_std_base();
1194 
1195  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1196  j_euler22_1.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1197 
1198  if( (j_euler22_1(0).get_etat() == ETATZERO)&&(j_euler22_1(1).get_etat() == ETATZERO)&&(j_euler22_1(2).get_etat()==ETATZERO))
1199  j_euler22_1 = 0;
1200 
1201 
1202 
1203  // Computation of elements used for the second version of the definition of In, Ip, En and Ep...
1204  // the (flat-space) first part angular-momentum 3-vector j_euler^i for fluid 1
1205  //-----------------------------------
1206  Tenseur Jph11_2(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1207 
1208  Tenseur A11_2(ent) ;
1209 
1210  if (eos.identify() == 3) {
1211 
1212 // A11_2 = gam_car() * nbar() * nbar() * eos.get_Knn_ent(ent(), ent2(), delta_car(), nzet) +
1213 // gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) ;
1214 
1215  A11_2 = gam_car() * nbar() * nbar() * K_nn() + gam_euler() * gam_euler2() * nbar() * nbar2() * K_np() ; //OK
1216 
1217 // A11_2 = gam_car() * K_nn() + gam_euler() * gam_euler2() * K_np() ; // Knn->Knn * nn * nn
1218 
1219  }
1220 
1221  else {
1222 
1223  A11_2 = gam_car() * nbar() * nbar() * eos.get_Knn(nbar(), nbar2(), delta_car(), nzet) +
1224  gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) ;
1225  }
1226 
1227  Jph11_2 = A11_2 * uuu ;
1228 
1229  j_euler11_2.set_etat_qcq();
1230 
1231  j_euler11_2.set(0) = 0; // r tetrad component
1232  j_euler11_2.set(1) = 0; // theta tetrad component
1233  j_euler11_2.set(2) = Jph11_2()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1234  j_euler11_2.set_triad (mp.get_bvect_spher());
1235  j_euler11_2.set_std_base();
1236 
1237  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1238  j_euler11_2.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1239 
1240  if( (j_euler11_2(0).get_etat() == ETATZERO)&&(j_euler11_2(1).get_etat() == ETATZERO)&&(j_euler11_2(2).get_etat()==ETATZERO))
1241  j_euler11_2 = 0;
1242 
1243  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1244  // the (flat-space) second part angular-momentum 3-vector j_euler^i for fluid 1
1245  //-----------------------------------
1246  Tenseur Jph12_2(mp);
1247 
1248  Tenseur A12_2(ent) ;
1249  Tenseur B12_2(ent) ;
1250 
1251  if (eos.identify() == 3) {
1252 // A12_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1253 // B12_2 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1254 
1255  A12_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * K_np(); //OK
1256  B12_2 = + gam_euler() * gam_euler2() * nbar() * nbar2() * K_np(); //OK
1257 
1258 // A12_2 = - gam_euler() * gam_euler2() * K_np(); // Knn --> Knn *nn *nn
1259 // B12_2 = + gam_euler() * gam_euler2() * K_np(); // Knn --> Knn *nn *nn
1260  }
1261 
1262  else {
1263  A12_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
1264  B12_2 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
1265  }
1266 
1267  Jph12_2 = A12_2 * uuu + B12_2 * uuu2;
1268 
1269  j_euler12_2.set_etat_qcq();
1270 
1271  j_euler12_2.set(0) = 0; // r tetrad component
1272  j_euler12_2.set(1) = 0; // theta tetrad component
1273  j_euler12_2.set(2) = Jph12_2()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1274  j_euler12_2.set_triad (mp.get_bvect_spher());
1275  j_euler12_2.set_std_base();
1276 
1277  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1278  j_euler12_2.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1279 
1280  if( (j_euler12_2(0).get_etat() == ETATZERO)&&(j_euler12_2(1).get_etat() == ETATZERO)&&(j_euler12_2(2).get_etat()==ETATZERO))
1281  j_euler12_2 = 0;
1282 
1283  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1284  // the (flat-space) first part angular-momentum 3-vector j_euler^i for fluid 2
1285  //-----------------------------------
1286  Tenseur Jph21_2(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1287 
1288  Tenseur A21_2(ent) ;
1289 
1290  if (eos.identify() == 3) {
1291 
1292 // A21_2 = gam2_car() * nbar2() * nbar2() * eos.get_Kpp_ent(ent(), ent2(), delta_car(), nzet) +
1293 // gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(), ent2(), delta_car(), nzet) ;
1294 
1295  A21_2 = gam2_car() * nbar2() * nbar2() * K_pp() + gam_euler() * gam_euler2() * nbar() * nbar2() * K_np() ; //OK
1296 
1297 // A21_2 = gam2_car() * K_pp() + gam_euler() * gam_euler2() * K_np() ; //Knn--> Knn *nn *nn
1298  }
1299 
1300  else {
1301 
1302  A21_2 = gam2_car() * nbar2() * nbar2() * eos.get_Kpp(nbar(), nbar2(), delta_car(), nzet) +
1303  gam_euler() * gam_euler2() * nbar2() * nbar() * eos.get_Knp(nbar(), nbar2(), delta_car(), nzet) ;
1304  }
1305 
1306  Jph21_2 = A21_2 * uuu2 ;
1307 
1308  j_euler21_2.set_etat_qcq();
1309 
1310  j_euler21_2.set(0) = 0; // r tetrad component
1311  j_euler21_2.set(1) = 0; // theta tetrad component
1312  j_euler21_2.set(2) = Jph21_2()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1313  j_euler21_2.set_triad (mp.get_bvect_spher());
1314  j_euler21_2.set_std_base();
1315 
1316  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1317  j_euler21_2.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1318 
1319  if( (j_euler21_2(0).get_etat() == ETATZERO)&&(j_euler21_2(1).get_etat() == ETATZERO)&&(j_euler21_2(2).get_etat()==ETATZERO))
1320  j_euler21_2 = 0;
1321 
1322  // Rajouté pour le calcul des moments cinétiques de chaque fluide
1323  // the (flat-space) second part angular-momentum 3-vector j_euler^i for fluid 1
1324  //-----------------------------------
1325  Tenseur Jph22_2(mp); // the normalized phi-component of J^i: Sqrt[g_phiphi]*J^phi
1326 
1327  Tenseur A22_2(ent) ;
1328  Tenseur B22_2(ent) ;
1329 
1330  if (eos.identify() == 3) {
1331 // A22_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1332 // B22_2 = + gam_euler() * gam_euler2() * nbar()*nbar2() * eos.get_Knp_ent(ent(),ent2(),delta_car(),nzet);
1333 
1334  A22_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * K_np(); //OK
1335  B22_2 = + gam_euler() * gam_euler2() * nbar()*nbar2() * K_np(); //OK
1336 
1337 // A22_2 = - gam_euler() * gam_euler2() * K_np(); //Knn --> Knn * nn * nn
1338 // B22_2 = + gam_euler() * gam_euler2() * K_np(); //Knn --> Knn * nn * nn
1339 
1340  }
1341 
1342  else {
1343  // some auxiliary EOS variables
1344  A22_2 = - gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet) ;
1345  B22_2 = + gam_euler() * gam_euler2() * nbar() * nbar2() * eos.get_Knp(nbar(),nbar2(),delta_car(),nzet);
1346  }
1347 
1348  Jph22_2 = A22_2 * uuu2 + B22_2 * uuu;
1349 
1350  j_euler22_2.set_etat_qcq();
1351 
1352  j_euler22_2.set(0) = 0; // r tetrad component
1353  j_euler22_2.set(1) = 0; // theta tetrad component
1354  j_euler22_2.set(2) = Jph22_2()/ bbb(); // phi tetrad component ... = J^phi r sin(th)
1355  j_euler22_2.set_triad (mp.get_bvect_spher());
1356  j_euler22_2.set_std_base();
1357 
1358  // RP: it seems that j_euler _HAS_ to have cartesian triad set on exit from here...!!
1359  j_euler22_2.change_triad( mp.get_bvect_cart() ) ; // Triad = Cartesian triad
1360 
1361  if( (j_euler22_2(0).get_etat() == ETATZERO)&&(j_euler22_2(1).get_etat() == ETATZERO)&&(j_euler22_2(2).get_etat()==ETATZERO))
1362  j_euler22_2 = 0;
1363 
1364 
1365 
1366 
1367  // The derived quantities are obsolete
1368  // -----------------------------------
1369  del_deriv() ;
1370 
1371 } // hydro_euler()
1372 
1373 }
1374 
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
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
Class for a two-fluid (tabulated) equation of state.
Definition: eos_bifluid.h:1534
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...
Definition: eos_bf_tabul.C:876
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
virtual void calcule_tout(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
Definition: eos_bifluid.C:281
double get_m1() const
Return the individual particule mass
Definition: eos_bifluid.h:256
double get_m2() const
Return the individual particule mass
Definition: eos_bifluid.h:262
Cmp get_Knn(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 1) .
Definition: eos_bifluid.C:745
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
Cmp get_Kpp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 2) .
Definition: eos_bifluid.C:767
Cmp get_Knp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy with respect to .
Definition: eos_bifluid.C:756
Class for two-fluid rotating relativistic stars.
Et_rot_bifluid(Map &mp_i, int nzet_i, bool relat, const Eos_bifluid &eos_i)
Standard constructor.
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
double * p_angu_mom_1_part2_2
To compute Xn (2nd version)
Tenseur gam_euler2
Lorentz factor between the fluid 2 and Eulerian observers
Tenseur enerps_euler
the combination : useful because in the Newtonian limit .
Tenseur K_np
Coefficient Knp.
virtual ~Et_rot_bifluid()
Destructor.
double * p_mass_b1
Baryon mass of fluid 1.
virtual double angu_mom_1_part1_1() const
To compute In (1st version)
double * p_mass_b2
Baryon mass of fluid 2.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Tenseur j_euler12_1
To compute Ip (1st version)
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
void set_enthalpies(const Cmp &, const Cmp &)
Sets both enthalpy profiles.
double * p_angu_mom_1_part1_2
To compute In (2nd version)
Tenseur j_euler2
To compute Jp.
Tenseur sphph_euler
The component of the stress tensor .
double * p_angu_mom_2_part2_1
To compute Xp (1st version)
Tenseur K_pp
Coefficient Kpp.
void operator=(const Et_rot_bifluid &)
Assignment to another Et_rot_bifluid.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
virtual double angu_mom_2_part1_2() const
To compute Ip (2nd version)
Tbl * p_xi_surf2
Description of the surface of fluid 2: 2-D Tbl containing the values of the radial coordinate on the...
virtual double angu_mom_2_part1_1() const
To compute Ip (1st version)
double * p_angu_mom_1
Angular momentum of fluid 1.
virtual double angu_mom_2() const
Angular momentum of fluid 2.
Tenseur j_euler11_1
To compute In (1st version)
virtual void del_deriv() const
Deletes all the derived quantities.
virtual double angu_mom_2_part2_2() const
To compute Xp (2nd version)
double * p_ray_eq2
Coordinate radius at , .
virtual void sauve(FILE *) const
Save in a file.
double * p_angu_mom_1_part2_1
To compute Xn (1st version)
virtual double mass_b() const
Total Baryon mass.
Tenseur j_euler1
To compute Jn.
Tenseur j_euler
Total angular momentum (flat-space!) 3-vector , which is related to of the "3+1" decomposition,...
virtual void equation_of_state()
Computes the proper baryon and energy densities, as well as pressure and the coefficients Knn,...
double * p_ray_pole2
Coordinate radius at .
Tenseur K_nn
Coefficient Knn.
virtual double mass_g() const
Gravitational mass.
virtual double angu_mom_2_part2_1() const
To compute Xp (1st version)
double * p_aplat2
Flatening r_pole/r_eq of fluid no.2.
double * p_angu_mom_2_part2_2
To compute Xp (2nd version)
double * p_area2
Surface area of fluid no.2.
double * p_angu_mom_2
Angular momentum of fluid 2.
Tenseur delta_car
The "relative velocity" (squared) of the two fluids.
double omega2
Rotation angular velocity for fluid 2 ([f_unit] )
const Eos_bifluid & eos
Equation of state for two-fluids model.
double * p_angu_mom_2_part1_2
To compute Ip (2nd version)
virtual double angu_mom_1_part1_2() const
To compute In (2nd version)
virtual double angu_mom_1_part2_1() const
To compute Xn (1st version)
Tenseur uuu2
Norm of the (fluid no.2) 3-velocity with respect to the eulerian observer.
double ray_eq2() const
Coordinate radius for fluid 2 at , [r_unit].
virtual double r_circ2() const
Circumferential radius for fluid 2.
virtual double mom_quad() const
Quadrupole moment.
virtual double aplat2() const
Flatening r_pole/r_eq for fluid 2.
virtual double grv2() const
Error on the virial identity GRV2.
double * p_angu_mom_1_part1_1
To compute In (1st version)
double * p_r_circ2
Circumferential radius of fluid no.2.
Tenseur ent2
Log-enthalpy for the second fluid.
double * p_ray_eq2_pi
Coordinate radius at , .
virtual double angu_mom_1_part2_2() const
To compute Xn (2nd version)
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
double * p_angu_mom_2_part1_1
To compute Ip (1st version)
Itbl * p_l_surf2
Description of the surface of fluid 2: 2-D Itbl containing the values of the domain index l on the su...
Tenseur nbar2
Baryon density in the fluid frame, for fluid 2.
virtual double angu_mom_1() const
Angular momentum of fluid 1.
virtual double angu_mom() const
Angular momentum.
double * p_ray_eq2_pis2
Coordinate radius at , .
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
virtual double mean_radius2() const
Mean radius for fluid 2.
virtual double mom_quad_old() const
Part of the quadrupole moment.
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition: etoile.h:1496
Tenseur uuu
Norm of u_euler.
Definition: etoile.h:1518
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile_rot.C:335
double omega
Rotation angular velocity ([f_unit] )
Definition: etoile.h:1501
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition: etoile_rot.C:404
virtual double r_circ() const
Circumferential radius.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: etoile_rot.C:364
Tenseur & logn
Metric potential = logn_auto.
Definition: etoile.h:1521
virtual double aplat() const
Flatening r_pole/r_eq.
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
Tenseur & dzeta
Metric potential = beta_auto.
Definition: etoile.h:1534
virtual double mean_radius() const
Mean radius.
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition: etoile_rot.C:389
virtual double z_eqf() const
Forward redshift factor at equator.
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
virtual double z_eqb() const
Backward redshift factor at equator.
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_rot.C:441
virtual double z_pole() const
Redshift factor at North pole.
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:460
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile.C:511
Tenseur press
Fluid pressure.
Definition: etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:465
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:468
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:457
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: etoile.h:442
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:209
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
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
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
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:500
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
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain array.
Definition: mtbl.h:118
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition: mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition: mtbl.h:221
Basic array class.
Definition: tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void sauve(FILE *) const
Save in a file.
Definition: tenseur.C:1325
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:900
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Definition: tenseur.C:650
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:70
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.