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