LORENE
binaire_global.C
1 /*
2  * Methods of class Binaire to compute global quantities
3  *
4  * (see file binaire.h for documentation)
5  */
6 
7 /*
8  * Copyright (c) 2000-2001 Eric Gourgoulhon
9  * Copyright (c) 2000-2001 Keisuke Taniguchi
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: binaire_global.C,v 1.8 2016/12/05 16:17:47 j_novak Exp $
34  * $Log: binaire_global.C,v $
35  * Revision 1.8 2016/12/05 16:17:47 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.7 2014/10/13 08:52:44 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.6 2004/03/25 10:28:59 j_novak
42  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
43  *
44  * Revision 1.5 2003/09/08 09:32:40 e_gourgoulhon
45  * Corrected a problem of spectral basis initialisation in virial_gb() and
46  * virial_fus(): introduced the new variable a1.
47  *
48  * Revision 1.4 2001/12/20 14:18:40 k_taniguchi
49  * Addition of the Komar mass, the virial error by Gourgoulhon and Bonazzola, and the virial error by Friedman, Uryu, and Shibata.
50  *
51  * Revision 1.3 2001/12/16 16:21:38 e_gourgoulhon
52  * #include "unites.h" is now local to Binaire::mass_adm(), in order
53  * not to make Lorene's units global variables.
54  *
55  * Revision 1.2 2001/12/14 09:45:14 k_taniguchi
56  * Correction of missing 16 pi G factor in the ADM mass
57  *
58  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
59  * LORENE
60  *
61  * Revision 2.3 2000/03/08 12:26:33 eric
62  * Ajout de l'appel a std_base_scal() sur le Cmp source dans le cas
63  * relativiste (masse ADM).
64  *
65  * Revision 2.2 2000/02/23 11:26:00 keisuke
66  * Changement de "virial relation".
67  *
68  * Revision 2.1 2000/02/18 15:48:55 eric
69  * *** empty log message ***
70  *
71  * Revision 2.0 2000/02/18 14:53:09 eric
72  * *** empty log message ***
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_global.C,v 1.8 2016/12/05 16:17:47 j_novak Exp $
76  *
77  */
78 
79 // Headers C
80 #include "math.h"
81 
82 // Headers Lorene
83 #include "binaire.h"
84 #include "unites.h"
85 
86  //---------------------------------//
87  // ADM mass //
88  //---------------------------------//
89 
90 namespace Lorene {
91 double Binaire::mass_adm() const {
92  using namespace Unites ;
93 
94  if (p_mass_adm == 0x0) { // a new computation is requireed
95 
96  p_mass_adm = new double ;
97 
98  if (star1.is_relativistic()) { // Relativistic case
99  // -----------------
100  assert( star2.is_relativistic() ) ;
101 
102  *p_mass_adm = 0 ;
103 
104  for (int i=0; i<=1; i++) { // loop on the stars
105 
106  const Cmp& a2 = (et[i]->get_a_car())() ;
107  const Cmp& ee = (et[i]->get_ener_euler())() ;
108  const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
109  const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
110 
111  Cmp source = pow(a2, 1.25) * ee
112  + pow(a2, 0.25) * (ak2_auto + ak2_comp) / (4.*qpig) ;
113 
114  source.std_base_scal() ;
115 
116  *p_mass_adm += source.integrale() ;
117 
118  }
119 
120  }
121  else { // Newtonian case
122  // --------------
123 
124  *p_mass_adm = star1.mass_b() + star2.mass_b() ;
125 
126  }
127 
128  } // End of the case where a new computation was necessary
129 
130  return *p_mass_adm ;
131 
132 }
133 
134 
135  //---------------------------------//
136  // Komar mass //
137  //---------------------------------//
138 
139 double Binaire::mass_kom() const {
140 
141  using namespace Unites ;
142 
143  if (p_mass_kom == 0x0) { // a new computation is requireed
144 
145  p_mass_kom = new double ;
146 
147  if (star1.is_relativistic()) { // Relativistic case
148  // -----------------
149  assert( star2.is_relativistic() ) ;
150 
151  *p_mass_kom = 0 ;
152 
153  for (int i=0; i<=1; i++) { // loop on the stars
154 
155  const Cmp& lapse = (et[i]->get_nnn())() ;
156  const Cmp& a2 = (et[i]->get_a_car())() ;
157  const Cmp& ee = (et[i]->get_ener_euler())() ;
158  const Cmp& se = (et[i]->get_s_euler())() ;
159  const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
160  const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
161 
162  const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
163  const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
164  const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
165  const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
166 
167  Tenseur dndb = flat_scalar_prod(dnu_auto,
168  dbe_auto + dbe_comp) ;
169  Tenseur dndn = flat_scalar_prod(dnu_auto,
170  dnu_auto + dnu_comp) ;
171 
172  Cmp source = lapse * ( a2 * (ee + se)
173  + (ak2_auto + ak2_comp)/qpig
174  - dndb()/qpig + dndn()/qpig ) ;
175 
176  source.std_base_scal() ;
177 
178  *p_mass_kom += source.integrale() ;
179 
180  }
181 
182  }
183  else { // Newtonian case
184  // --------------
185 
186  *p_mass_kom = star1.mass_b() + star2.mass_b() ;
187 
188  }
189 
190  } // End of the case where a new computation was necessary
191 
192  return *p_mass_kom ;
193 
194 }
195 
196 
197  //---------------------------------//
198  // Total angular momentum //
199  //---------------------------------//
200 
201 const Tbl& Binaire::angu_mom() const {
202 
203  if (p_angu_mom == 0x0) { // a new computation is requireed
204 
205  p_angu_mom = new Tbl(3) ;
206 
207  p_angu_mom->annule_hard() ; // fills the double array with zeros
208 
209  for (int i=0; i<=1; i++) { // loop on the stars
210 
211  const Map& mp = et[i]->get_mp() ;
212 
213  Cmp xx(mp) ;
214  Cmp yy(mp) ;
215  Cmp zz(mp) ;
216 
217  xx = mp.xa ;
218  yy = mp.ya ;
219  zz = mp.za ;
220 
221  const Cmp& vx = (et[i]->get_u_euler())(0) ;
222  const Cmp& vy = (et[i]->get_u_euler())(1) ;
223  const Cmp& vz = (et[i]->get_u_euler())(2) ;
224 
225  Cmp rho(mp) ;
226 
227  if ( et[i]->is_relativistic() ) {
228  const Cmp& a2 = (et[i]->get_a_car())() ;
229  const Cmp& ee = (et[i]->get_ener_euler())() ;
230  const Cmp& pp = (et[i]->get_press())() ;
231  rho = pow(a2, 2.5) * (ee + pp) ;
232  }
233  else {
234  rho = (et[i]->get_nbar())() ;
235  }
236 
237 
238  Base_val** base = (et[i]->get_mp()).get_mg()->std_base_vect_cart() ;
239 
240  // X component
241  // -----------
242 
243  Cmp source = rho * ( yy * vz - zz * vy ) ;
244 
245  (source.va).set_base( *(base[2]) ) ; // same basis as V^z
246 
247 //## p_angu_mom->set(0) += source.integrale() ;
248 
249  p_angu_mom->set(0) += 0 ;
250 
251  // y component
252  // -----------
253 
254  source = rho * ( zz * vx - xx * vz ) ;
255 
256  (source.va).set_base( *(base[2]) ) ; // same basis as V^z
257 
258 //## p_angu_mom->set(1) += source.integrale() ;
259  p_angu_mom->set(1) += 0 ;
260 
261 
262  // Z component
263  // -----------
264 
265  source = rho * ( xx * vy - yy * vx ) ;
266 
267  source.std_base_scal() ; // same basis as V^x (standard scalar
268  // field)
269 
270  p_angu_mom->set(2) += source.integrale() ;
271 
272  delete base[0] ;
273  delete base[1] ;
274  delete base[2] ;
275  delete [] base ;
276 
277  } // End of the loop on the stars
278 
279  } // End of the case where a new computation was necessary
280 
281  return *p_angu_mom ;
282 
283 }
284 
285 
286 
287 
288  //---------------------------------//
289  // Total energy //
290  //---------------------------------//
291 
292 double Binaire::total_ener() const {
293 
294  if (p_total_ener == 0x0) { // a new computation is requireed
295 
296  p_total_ener = new double ;
297 
298  if (star1.is_relativistic()) { // Relativistic case
299  // -----------------
300 
301  assert( star2.is_relativistic() ) ;
302 
303  *p_total_ener = mass_adm() - star1.mass_b() - star2.mass_b() ;
304 
305  }
306  else { // Newtonian case
307  // --------------
308 
309  *p_total_ener = 0 ;
310 
311  for (int i=0; i<=1; i++) { // loop on the stars
312 
313  const Cmp e_int = (et[i]->get_ener())()
314  - (et[i]->get_nbar())() ;
315 
316  const Cmp& rho = (et[i]->get_nbar())() ;
317 
318  // Fluid velocity with respect to the inertial frame
319  const Tenseur& vit = et[i]->get_u_euler() ;
320 
321  Cmp vit2 = flat_scalar_prod(vit, vit)() ;
322 
323  // Gravitational potential
324  const Cmp nu = (et[i]->get_logn_auto())()
325  + (et[i]->get_logn_comp())() ;
326 
327  Cmp source = e_int + .5 * rho * vit2 + .5 * rho * nu ;
328 
329  *p_total_ener += source.integrale() ;
330 
331 
332  } // End of the loop on the stars
333 
334  } // End of Newtonian case
335 
336  } // End of the case where a new computation was necessary
337 
338 
339  return *p_total_ener ;
340 
341 }
342 
343 
344  //---------------------------------//
345  // Error on the virial theorem //
346  //---------------------------------//
347 
348 double Binaire::virial() const {
349 
350  if (p_virial == 0x0) { // a new computation is requireed
351 
352  p_virial = new double ;
353 
354  if (star1.is_relativistic()) { // Relativistic case
355  // -----------------
356 
357  assert( star2.is_relativistic() ) ;
358 
359  *p_virial = 1. - mass_kom() / mass_adm() ;
360 
361  }
362  else { // Newtonian case
363  // --------------
364 
365  *p_virial = 0 ;
366 
367 
368  double vir_mat = 0 ;
369  double vir_grav = 0 ;
370 
371  for (int i=0; i<=1; i++) { // loop on the stars
372 
373  const Cmp& pp = (et[i]->get_press())() ;
374 
375  const Cmp& rho = (et[i]->get_nbar())() ;
376 
377  // Fluid velocity with respect to the inertial frame
378  const Tenseur& vit = et[i]->get_u_euler() ;
379 
380  Cmp vit2 = flat_scalar_prod(vit, vit)() ;
381 
382  // Gravitational potential
383  const Cmp nu = (et[i]->get_logn_auto())()
384  + (et[i]->get_logn_comp())() ;
385 
386  Cmp source = 3*pp + rho * vit2 ;
387 
388  vir_mat += source.integrale() ;
389 
390  source = .5 * rho * nu ;
391 
392  vir_grav += source.integrale() ;
393 
394  } // End of the loop on the stars
395 
396  *p_virial = ( vir_mat + vir_grav ) / fabs(vir_grav) ;
397 
398  } // End of the Newtonian case
399 
400  } // End of the case where a new computation was necessary
401 
402  return *p_virial ;
403 
404 }
405 
406 
407  //----------------------------------------------//
408  // Virial error by Gourgoulhon and Bonazzola //
409  //----------------------------------------------//
410 
411 double Binaire::virial_gb() const {
412 
413  using namespace Unites ;
414 
415  if (p_virial_gb == 0x0) { // a new computation is requireed
416 
417  p_virial_gb = new double ;
418 
419  if (star1.is_relativistic()) { // Relativistic case
420  // -----------------
421 
422  assert( star2.is_relativistic() ) ;
423 
424  *p_virial_gb = 0 ;
425 
426  double vir_pres = 0. ;
427  double vir_extr = 0. ;
428  double vir_grav = 0. ;
429 
430  for (int i=0; i<=1; i++) { // loop on the stars
431 
432  const Cmp& a2 = (et[i]->get_a_car())() ;
433  const Cmp& se = (et[i]->get_s_euler())() ;
434  const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
435  const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
436 
437  const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
438  const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
439  const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
440  const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
441 
442  Cmp a1 = sqrt(a2) ;
443  a1.std_base_scal() ;
444 
445  Cmp source = 2. * a2 * a1 * se ;
446  vir_pres += source.integrale() ;
447 
448  source = 1.5 * a1 * (ak2_auto + ak2_comp) / qpig ;
449  source.std_base_scal() ;
450  vir_extr += source.integrale() ;
451 
452  Tenseur sprod1 = flat_scalar_prod(dbe_auto, dbe_auto+dbe_comp) ;
453  Tenseur sprod2 = flat_scalar_prod(dnu_auto, dnu_auto+dnu_comp) ;
454  Tenseur sprod3 = flat_scalar_prod(dbe_auto, dnu_auto+dnu_comp) ;
455 
456  source = a1 * ( sprod1() - sprod2() - 2.*sprod3() )/qpig ;
457  vir_grav += source.integrale() ;
458 
459  } // End of the loop on the stars
460 
461 
462  *p_virial_gb = (vir_pres + vir_extr + vir_grav) / mass_adm() ;
463 
464  }
465  else { // Newtonian case
466  // --------------
467 
468  *p_virial_gb = virial() ;
469 
470 
471  } // End of the Newtonian case
472 
473  } // End of the case where a new computation was necessary
474 
475  return *p_virial_gb ;
476 
477 }
478 
479 
480  //------------------------------------------------//
481  // Virial error by Friedman, Uryu, and Shibata //
482  //------------------------------------------------//
483 
484 double Binaire::virial_fus() const {
485 
486  using namespace Unites ;
487 
488  if (p_virial_fus == 0x0) { // a new computation is requireed
489 
490  p_virial_fus = new double ;
491 
492  if (star1.is_relativistic()) { // Relativistic case
493  // -----------------
494 
495  assert( star2.is_relativistic() ) ;
496 
497  *p_virial_fus = 0 ;
498 
499  double vir_pres = 0. ;
500  double vir_extr = 0. ;
501  double vir_grav = 0. ;
502 
503  for (int i=0; i<=1; i++) { // loop on the stars
504 
505  const Cmp& lapse = (et[i]->get_nnn())() ;
506  const Cmp& a2 = (et[i]->get_a_car())() ;
507  const Cmp& se = (et[i]->get_s_euler())() ;
508  const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
509  const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
510 
511  const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
512  const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
513  const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
514  const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
515 
516  Cmp a1 = sqrt(a2) ;
517  a1.std_base_scal() ;
518 
519  Cmp source = 2. * lapse * a2 * a1 * se ;
520  vir_pres += source.integrale() ;
521 
522  source = 1.5 * lapse * a1 * (ak2_auto + ak2_comp) / qpig ;
523  vir_extr += source.integrale() ;
524 
525  Tenseur sprod = flat_scalar_prod( dbe_auto, dbe_auto+dbe_comp )
526  - flat_scalar_prod( dnu_auto, dnu_auto+dnu_comp ) ;
527 
528  source = lapse * a1 * sprod() / qpig ;
529  vir_grav += source.integrale() ;
530 
531  } // End of the loop on the stars
532 
533 
534  *p_virial_fus = (vir_pres + vir_extr + vir_grav) / mass_adm() ;
535 
536  }
537  else { // Newtonian case
538  // --------------
539 
540  *p_virial_fus = virial() ;
541 
542 
543  } // End of the Newtonian case
544 
545  } // End of the case where a new computation was necessary
546 
547  return *p_virial_fus ;
548 
549 }
550 
551 }
Coord xa
Absolute x coordinate.
Definition: map.h:742
double total_ener() const
Total energy (excluding the rest mass energy).
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:662
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition: etoile.h:736
Etoile_bin star1
First star of the system.
Definition: binaire.h:118
double mass_adm() const
Total ADM mass.
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene prototypes.
Definition: app_hor.h:67
const Tenseur & get_akcar_comp() const
Returns the part of the scalar generated by shift_auto and shift_comp , i.e.
Definition: etoile.h:1213
Standard units of space, time and mass.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
const Tenseur & get_d_beta_auto() const
Returns the gradient of beta_auto (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1144
double * p_virial
Virial theorem error.
Definition: binaire.h:154
Base class for coordinate mappings.
Definition: map.h:682
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
const Tenseur & get_logn_auto() const
Returns the logarithm of the part of the lapse N generated principaly by the star.
Definition: etoile.h:704
double virial_gb() const
Estimates the relative error on the virial theorem calculated by E.Gourgoulhon and S...
const Tenseur & get_d_logn_auto() const
Returns the gradient of logn_auto (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1124
double * p_total_ener
Total energy of the system.
Definition: binaire.h:151
const Tenseur & get_press() const
Returns the fluid pressure.
Definition: etoile.h:685
const Tenseur & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: etoile.h:1119
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): { et[0]} contains the address of { star1} and...
Definition: binaire.h:127
double * p_virial_fus
Virial theorem error by J.L.Friedman, K.Uryu, and M.Shibata.
Definition: binaire.h:160
const Tenseur & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:691
double integrale() const
Computes the integral over all space of *this .
Definition: cmp_integ.C:58
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition: etoile.h:730
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Tenseur & get_d_beta_comp() const
Returns the gradient of beta_comp (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1149
const Tbl & angu_mom() const
Total angular momentum.
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: etoile.h:688
Tbl * p_angu_mom
Total angular momentum of the system.
Definition: binaire.h:148
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
const Tenseur & get_d_logn_comp() const
Returns the gradient of logn_comp (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1134
Coord ya
Absolute y coordinate.
Definition: map.h:743
Bases of the spectral expansions.
Definition: base_val.h:325
const Tenseur & get_akcar_auto() const
Returns the part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:1207
double mass_kom() const
Total Komar mass.
double virial_fus() const
Estimates the relative error on the virial theorem calculated by J.L.Friedman, K.Uryu, and M.Shibata (PRD accepted, gr-qc/0108070)
double * p_mass_adm
Total ADM mass of the system.
Definition: binaire.h:142
Coord za
Absolute z coordinate.
Definition: map.h:744
double * p_mass_kom
Total Komar mass of the system.
Definition: binaire.h:145
double * p_virial_gb
Virial theorem error by E.Gourgoulhon and S.Bonazzola.
Definition: binaire.h:157
const Tenseur & get_ener() const
Returns the proper total energy density.
Definition: etoile.h:682
double virial() const
Estimates the relative error on the virial theorem (for a relativistic one, it returns $|1 - M_{ Koma...
Basic array class.
Definition: tbl.h:164
const Tenseur & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:697
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
virtual double mass_b() const
Baryon mass.
const Tenseur & get_nbar() const
Returns the proper baryon density.
Definition: etoile.h:679
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:670
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Etoile_bin star2
Second star of the system.
Definition: binaire.h:121