LORENE
binary_global.C
1 /*
2  * Methods of class Binary to compute global quantities
3  *
4  * (see file binary.h for documentation)
5  */
6 
7 /*
8  * Copyright (c) 2004 Francois Limousin
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 char binary_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary/binary_global.C,v 1.16 2014/10/13 08:52:45 j_novak Exp $" ;
30 
31 /*
32  * $Id: binary_global.C,v 1.16 2014/10/13 08:52:45 j_novak Exp $
33  * $Log: binary_global.C,v $
34  * Revision 1.16 2014/10/13 08:52:45 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.15 2006/08/01 14:26:50 f_limousin
38  * Small changes
39  *
40  * Revision 1.14 2006/04/11 14:25:15 f_limousin
41  * New version of the code : improvement of the computation of some
42  * critical sources, estimation of the dirac gauge, helical symmetry...
43  *
44  * Revision 1.13 2005/09/18 13:13:41 f_limousin
45  * Extension of vphi in the compactified domain for the computation
46  * of J_ADM by a volume integral.
47  *
48  * Revision 1.12 2005/09/15 14:41:04 e_gourgoulhon
49  * The total angular momentum is now computed via a volume integral.
50  *
51  * Revision 1.11 2005/09/13 19:38:31 f_limousin
52  * Reintroduction of the resolution of the equations in cartesian coordinates.
53  *
54  * Revision 1.10 2005/04/08 12:36:45 f_limousin
55  * Just to avoid warnings...
56  *
57  * Revision 1.9 2005/02/17 17:35:00 f_limousin
58  * Change the name of some quantities to be consistent with other classes
59  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
60  *
61  * Revision 1.8 2004/07/21 11:46:24 f_limousin
62  * Add function mass_adm_vol() to compute the ADM mass of the system
63  * with a volume integral instead of a surface one.
64  *
65  * Revision 1.7 2004/05/25 14:25:53 f_limousin
66  * Add the virial theorem for conformally flat configurations.
67  *
68  * Revision 1.6 2004/03/31 12:44:54 f_limousin
69  * Minor modifs.
70  *
71  * Revision 1.5 2004/03/25 10:29:01 j_novak
72  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
73  *
74  * Revision 1.4 2004/02/27 10:25:30 f_limousin
75  * Modif. to avoid an error in compilation.
76  *
77  * Revision 1.3 2004/02/27 10:03:04 f_limousin
78  * The computation of mass_adm() and mass_komar() is now OK !
79  *
80  * Revision 1.2 2004/01/20 15:21:36 f_limousin
81  * First version
82  *
83  *
84  * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_global.C,v 1.16 2014/10/13 08:52:45 j_novak Exp $
85  *
86  */
87 
88 
89 // Headers C
90 #include "math.h"
91 
92 // Headers Lorene
93 #include "nbr_spx.h"
94 #include "binary.h"
95 #include "unites.h"
96 #include "metric.h"
97 
98  //---------------------------------//
99  // ADM mass //
100  //---------------------------------//
101 
102 namespace Lorene {
103 double Binary::mass_adm() const {
104 
105  using namespace Unites ;
106  if (p_mass_adm == 0x0) { // a new computation is requireed
107 
108  p_mass_adm = new double ;
109 
110  *p_mass_adm = 0 ;
111 
112  const Map_af map0 (et[0]->get_mp()) ;
113  const Metric& flat = (et[0]->get_flat()) ;
114 
115  Vector dpsi(0.5*(et[0]->get_lnq() -
116  et[0]->get_logn()).derive_cov(flat)) ;
117 
118  Vector ww (0.125*(contract(et[0]->get_hij().derive_cov(flat), 1, 2)
119  - (et[0]->get_hij().trace(flat)).derive_con(flat))) ;
120 
121  dpsi.change_triad(map0.get_bvect_spher()) ;
122  ww.change_triad(map0.get_bvect_spher()) ;
123 
124  // ww = 0 in Dirac gauge (Eq 174 of BGGN)
125  Scalar integrand (dpsi(1) + 0*ww(1)) ;
126 
127  *p_mass_adm = map0.integrale_surface_infini (integrand) / (-qpig/2.) ;
128 
129  } // End of the case where a new computation was necessary
130 
131  return *p_mass_adm ;
132 
133 }
134 
135 double Binary::mass_adm_vol() const {
136 
137  using namespace Unites ;
138 
139  double massadm ;
140  massadm = 0. ;
141 
142  for (int i=0; i<=1; i++) { // loop on the stars
143 
144  // Declaration of all fields
145  const Scalar& psi4 = et[i]->get_psi4() ;
146  Scalar psi (pow(psi4, 0.25)) ;
147  psi.std_spectral_base() ;
148  const Scalar& ener_euler = et[i]->get_ener_euler() ;
149  const Scalar& kcar_auto = et[i]->get_kcar_auto() ;
150  const Scalar& kcar_comp = et[i]->get_kcar_comp() ;
151  const Metric& gtilde = et[i]->get_gtilde() ;
152  const Metric& flat = et[i]->get_flat() ;
153  const Sym_tensor& hij = et[i]->get_hij() ;
154  const Sym_tensor& hij_auto = et[i]->get_hij_auto() ;
155  const Vector& dcov_logn = et[i]->get_dcov_logn() ;
156  const Vector& dcov_phi = et[i]->get_dcov_phi() ;
157  const Vector& dcov_lnq = 2*dcov_phi + dcov_logn ;
158  const Scalar& lnq_auto = et[i]->get_lnq_auto() ;
159  const Scalar& logn_auto = et[i]->get_logn_auto() ;
160  const Scalar& phi_auto = 0.5 * (lnq_auto - logn_auto) ;
161 
162  const Tensor& dcov_hij_auto = hij_auto.derive_cov(flat) ;
163  const Tensor& dcov_gtilde = gtilde.cov().derive_cov(flat) ;
164  const Tensor& dcov_phi_auto = phi_auto.derive_cov(flat) ;
165  const Tensor& dcov_logn_auto = logn_auto.derive_cov(flat) ;
166  const Tensor& dcov_lnq_auto = lnq_auto.derive_cov(flat) ;
167  Tensor dcovdcov_lnq_auto = lnq_auto.derive_cov(flat).derive_cov(flat) ;
168  dcovdcov_lnq_auto.inc_dzpuis() ;
169  Tensor dcovdcov_logn_auto = logn_auto.derive_cov(flat).derive_cov(flat) ;
170  dcovdcov_logn_auto.inc_dzpuis() ;
171 
172  // Source in IWM approximation
173  Scalar source = - psi4 % (qpig*ener_euler + (kcar_auto + kcar_comp)/4.)
174  - 0*2*contract(contract(gtilde.con(), 0, dcov_phi, 0),
175  0, dcov_phi_auto, 0, true) ;
176 
177  // Source = 0 in IWM
178  source += 4*contract(hij, 0, 1, dcov_logn * dcov_phi_auto, 0, 1) +
179  2*contract(hij, 0, 1, dcov_phi * dcov_phi_auto, 0, 1) +
180  0.0625 * contract(gtilde.con(), 0, 1, contract(
181  dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) -
182  0.125 * contract(gtilde.con(), 0, 1, contract(dcov_hij_auto,
183  0, 1, dcov_gtilde, 0, 2), 0, 1) -
184  contract(hij,0,1,dcovdcov_lnq_auto + dcov_lnq_auto*dcov_lnq,0,1) +
185  contract(hij,0,1,dcovdcov_logn_auto + dcov_logn_auto*dcov_logn,0,1) ;
186 
187  source = source * psi ;
188 
189  source.std_spectral_base() ;
190 
191  massadm += - source.integrale()/qpig ;
192  }
193 
194  return massadm ;
195 }
196 
197  //---------------------------------//
198  // Komar mass //
199  //---------------------------------//
200 
201 double Binary::mass_kom() const {
202 
203  using namespace Unites ;
204 
205  if (p_mass_kom == 0x0) { // a new computation is requireed
206 
207  p_mass_kom = new double ;
208 
209  *p_mass_kom = 0 ;
210 
211  const Tensor& logn = et[0]->get_logn() ;
212  const Metric& flat = (et[0]->get_flat()) ;
213  const Sym_tensor& hij = (et[0]->get_hij()) ;
214  Map_af map0 (et[0]->get_mp()) ;
215 
216  Vector vect = logn.derive_con(flat) +
217  contract(hij, 1, logn.derive_cov(flat), 0) ;
218  vect.change_triad(map0.get_bvect_spher()) ;
219  Scalar integrant (vect(1)) ;
220 
221  *p_mass_kom = map0.integrale_surface_infini (integrant) / qpig ;
222 
223  } // End of the case where a new computation was necessary
224 
225  return *p_mass_kom ;
226 
227 }
228 
229 double Binary::mass_kom_vol() const {
230 
231  using namespace Unites ;
232 
233  double masskom ;
234  masskom = 0. ;
235 
236  for (int i=0; i<=1; i++) { // loop on the stars
237 
238  // Declaration of all fields
239  const Scalar& psi4 = et[i]->get_psi4() ;
240  const Scalar& ener_euler = et[i]->get_ener_euler() ;
241  const Scalar& s_euler = et[i]->get_s_euler() ;
242  const Scalar& kcar_auto = et[i]->get_kcar_auto() ;
243  const Scalar& kcar_comp = et[i]->get_kcar_comp() ;
244  const Metric& gtilde = et[i]->get_gtilde() ;
245  const Metric& flat = et[i]->get_flat() ;
246  const Sym_tensor& hij = et[i]->get_hij() ;
247  const Scalar& logn = et[i]->get_logn_auto() + et[i]->get_logn_comp() ;
248  const Scalar& logn_auto = et[i]->get_logn_auto() ;
249  Scalar nn = exp(logn) ;
250  nn.std_spectral_base() ;
251 
252  const Tensor& dcov_logn_auto = logn_auto.derive_cov(flat) ;
253  const Vector& dcov_logn = et[i]->get_dcov_logn() ;
254  const Vector& dcon_logn = et[i]->get_dcon_logn() ;
255  const Vector& dcov_phi = et[i]->get_dcov_phi() ;
256  Tensor dcovdcov_logn_auto = (logn_auto.derive_cov(flat))
257  .derive_cov(flat) ;
258  dcovdcov_logn_auto.inc_dzpuis() ;
259 
260  Scalar source = qpig * psi4 % (ener_euler + s_euler) ;
261  source += psi4 % (kcar_auto + kcar_comp) ;
262  source += - 0*contract(dcov_logn_auto, 0, dcon_logn, 0, true)
263  - 2. * contract(contract(gtilde.con(), 0, dcov_phi, 0), 0,
264  dcov_logn_auto, 0, true) ;
265  source += - contract(hij, 0, 1, dcovdcov_logn_auto +
266  dcov_logn_auto*dcov_logn, 0, 1) ;
267 
268  source = source / qpig * nn ;
269 
270  source.std_spectral_base() ;
271 
272  masskom += source.integrale() ;
273 
274  }
275 
276  return masskom ;
277 
278 }
279 
280 
281  //---------------------------------//
282  // Total angular momentum //
283  //---------------------------------//
284 
285 const Tbl& Binary::angu_mom() const {
286 
287  using namespace Unites ;
288 
289  /*
290  if (p_angu_mom == 0x0) { // a new computation is requireed
291 
292  p_angu_mom = new Tbl(3) ;
293 
294  p_angu_mom->annule_hard() ; // fills the double array with zeros
295 
296  const Sym_tensor& kij_auto = et[0]->get_tkij_auto() ;
297  const Sym_tensor& kij_comp = et[0]->get_tkij_comp() ;
298  const Tensor& psi4 = et[0]->get_psi4() ;
299  const Map_af map0 (kij_auto.get_mp()) ;
300 
301  Sym_tensor kij = (kij_auto + kij_comp) / psi4 ;
302  kij.change_triad(map0.get_bvect_cart()) ;
303 
304  // X component
305  // -----------
306 
307  Vector vect_x(et[0]->get_mp(), CON, map0.get_bvect_cart()) ;
308 
309  for (int i=1; i<=3; i++) {
310 
311  Scalar kij_1 = kij(3, i) ;
312  Scalar kij_2 = kij(2, i) ;
313 
314  kij_1.mult_rsint() ;
315  Valeur vtmp = kij_1.get_spectral_va().mult_sp() ;
316  kij_1.set_spectral_va() = vtmp ;
317 
318  kij_2.mult_r() ;
319  vtmp = kij_2.get_spectral_va().mult_ct() ;
320  kij_2.set_spectral_va() = vtmp ;
321 
322  vect_x.set(i) = kij_1 - kij_2 ;
323  }
324 
325  vect_x.change_triad(map0.get_bvect_spher()) ;
326  Scalar integrant_x (vect_x(1)) ;
327 
328  p_angu_mom->set(0) = map0.integrale_surface_infini (integrant_x)
329  / (8*M_PI) ;
330 
331  // Y component
332  // -----------
333 
334  Vector vect_y(et[0]->get_mp(), CON, map0.get_bvect_cart()) ;
335 
336  for (int i=1; i<=3; i++) {
337 
338  Scalar kij_1 = kij(1, i) ;
339  Scalar kij_2 = kij(3, i) ;
340 
341  kij_1.mult_r() ;
342  Valeur vtmp = kij_1.get_spectral_va().mult_ct() ;
343  kij_1.set_spectral_va() = vtmp ;
344 
345  kij_2.mult_rsint() ;
346  vtmp = kij_2.get_spectral_va().mult_cp() ;
347  kij_2.set_spectral_va() = vtmp ;
348 
349  vect_y.set(i) = kij_1 - kij_2 ;
350  }
351 
352  vect_y.change_triad(map0.get_bvect_spher()) ;
353  Scalar integrant_y (vect_y(1)) ;
354 
355  p_angu_mom->set(1) = map0.integrale_surface_infini (integrant_y)
356  / (8*M_PI) ;
357 
358  // Z component
359  // -----------
360 
361  Vector vect_z(et[0]->get_mp(), CON, map0.get_bvect_cart()) ;
362 
363  for (int i=1; i<=3; i++) {
364 
365  Scalar kij_1 = kij(2, i) ;
366  Scalar kij_2 = kij(1, i) ;
367 
368  kij_1.mult_rsint() ;
369  Valeur vtmp = kij_1.get_spectral_va().mult_cp() ;
370  kij_1.set_spectral_va() = vtmp ;
371 
372  kij_2.mult_rsint() ;
373  vtmp = kij_2.get_spectral_va().mult_sp() ;
374  kij_2.set_spectral_va() = vtmp ;
375 
376  vect_z.set(i) = kij_1 - kij_2 ;
377  }
378 
379  vect_z.change_triad(map0.get_bvect_spher()) ;
380  Scalar integrant_z (vect_z(1)) ;
381 
382  p_angu_mom->set(2) = map0.integrale_surface_infini (integrant_z)
383  ;// (8*M_PI) ;
384 
385 
386  } // End of the case where a new computation was necessary
387  */
388 
389 
390  /*
391  if (p_angu_mom == 0x0) { // a new computation is requireed
392  p_angu_mom = new Tbl(3) ;
393  p_angu_mom->annule_hard() ; // fills the double array with zeros
394  p_angu_mom->set(0) = 0. ;
395  p_angu_mom->set(1) = 0. ;
396 
397  // Alignement
398  double orientation_un = et[0]->get_mp().get_rot_phi() ;
399  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
400  double orientation_deux = et[1]->get_mp().get_rot_phi() ;
401  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
402  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
403 
404  // Construction of an auxiliar grid and mapping
405  int nzones = et[0]->get_mp().get_mg()->get_nzone() ;
406  double* bornes = new double [nzones+1] ;
407  double courant = (et[0]->get_mp().get_ori_x()-et[0]->get_mp().get_ori_x())+1 ;
408  for (int i=nzones-1 ; i>0 ; i--) {
409  bornes[i] = courant ;
410  courant /= 2. ;
411  }
412  bornes[0] = 0 ;
413  bornes[nzones] = __infinity ;
414 
415  Map_af mapping (*(et[0]->get_mp().get_mg()), bornes) ;
416 
417  delete [] bornes ;
418 
419  // Construction of k_total
420  Sym_tensor k_total (mapping, CON, mapping.get_bvect_cart()) ;
421 
422  Vector shift_un (mapping, CON, mapping.get_bvect_cart()) ;
423  Vector shift_deux (mapping, CON, mapping.get_bvect_cart()) ;
424 
425  Vector beta_un (et[0]->get_beta_auto()) ;
426  Vector beta_deux (et[1]->get_beta_auto()) ;
427  beta_un.change_triad(et[0]->get_mp().get_bvect_cart()) ;
428  beta_deux.change_triad(et[1]->get_mp().get_bvect_cart()) ;
429  beta_un.std_spectral_base() ;
430  beta_deux.std_spectral_base() ;
431 
432  shift_un.set(1).import(beta_un(1)) ;
433  shift_un.set(2).import(beta_un(2)) ;
434  shift_un.set(3).import(beta_un(3)) ;
435 
436  shift_deux.set(1).import(same_orient*beta_deux(1)) ;
437  shift_deux.set(2).import(same_orient*beta_deux(2)) ;
438  shift_deux.set(3).import(beta_deux(3)) ;
439 
440  Vector shift_tot (shift_un+shift_deux) ;
441  shift_tot.std_spectral_base() ;
442  shift_tot.annule(0, nzones-2) ;
443 
444 
445  // Substract the residuals
446  shift_tot.inc_dzpuis(2) ;
447  shift_tot.dec_dzpuis(2) ;
448 
449 
450  Sym_tensor temp_gamt (et[0]->get_gtilde().cov()) ;
451  temp_gamt.change_triad(mapping.get_bvect_cart()) ;
452  Metric gamt_cart (temp_gamt) ;
453 
454  k_total = shift_tot.ope_killing_conf(gamt_cart) / 2. ;
455 
456  for (int lig=1 ; lig<=3 ; lig++)
457  for (int col=lig ; col<=3 ; col++)
458  k_total.set(lig, col).mult_r_ced() ;
459 
460  Vector vecteur_un (mapping, CON, mapping.get_bvect_cart()) ;
461  for (int i=1 ; i<=3 ; i++)
462  vecteur_un.set(i) = k_total(1, i) ;
463  vecteur_un.change_triad (mapping.get_bvect_spher()) ;
464  Scalar integrant_un (vecteur_un(1)) ;
465 
466  Vector vecteur_deux (mapping, CON, mapping.get_bvect_cart()) ;
467  for (int i=1 ; i<=3 ; i++)
468  vecteur_deux.set(i) = k_total(2, i) ;
469  vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
470  Scalar integrant_deux (vecteur_deux(1)) ;
471 
472  // Multiplication by y and x :
473  integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
474  .mult_st() ;
475  integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
476  .mult_sp() ;
477 
478  integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
479  .mult_st() ;
480  integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
481  .mult_cp() ;
482 
483  p_angu_mom->set(2) = mapping.integrale_surface_infini (-integrant_un
484  +integrant_deux) / (2*qpig) ;
485 
486  }
487 
488  */
489 
490  if (p_angu_mom == 0x0) { // a new computation is requireed
491 
492  p_angu_mom = new Tbl(3) ;
493  p_angu_mom->annule_hard() ; // fills the double array with zeros
494 
495  // Reference Cartesian vector basis of the Absolute frame
496  Base_vect_cart bvect_ref(0.) ; // 0. = parallel to the Absolute frame
497 
498  for (int i=0; i<=1; i++) { // loop on the stars
499 
500  const Map& mp = et[i]->get_mp() ;
501  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
502 
503  // Function exp(-(r-r_0)^2/sigma^2)
504  // --------------------------------
505 
506  double r0 = mp.val_r(nzm1-1, 1, 0, 0) ;
507  double sigma = 1.*r0 ;
508 
509  Scalar rr (mp) ;
510  rr = mp.r ;
511 
512  Scalar ff (mp) ;
513  ff = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
514  for (int ii=0; ii<nzm1; ii++)
515  ff.set_domain(ii) = 1. ;
516  ff.set_outer_boundary(nzm1, 0) ;
517  ff.std_spectral_base() ;
518 
519  // Azimuthal vector d/dphi
520  Vector vphi(mp, CON, bvect_ref) ;
521  Scalar yya (mp) ;
522  yya = mp.ya ;
523  Scalar xxa (mp) ;
524  xxa = mp.xa ;
525  vphi.set(1) = - yya * ff ; // phi^X
526  vphi.set(2) = xxa * ff ;
527  vphi.set(3) = 0 ;
528 
529  vphi.set(1).set_outer_boundary(nzm1, 0) ;
530  vphi.set(2).set_outer_boundary(nzm1, 0) ;
531 
532  vphi.std_spectral_base() ;
533  vphi.change_triad(mp.get_bvect_cart()) ;
534 
535  // Matter part
536  // -----------
537  const Scalar& ee = et[i]->get_ener_euler() ; // E
538  const Scalar& pp = et[i]->get_press() ; // p
539  const Scalar& psi4 = et[i]->get_psi4() ; // Psi^4
540  Scalar rho = pow(psi4, double(2.5)) * (ee+pp) ;
541  rho.std_spectral_base() ;
542 
543  Vector jmom = rho * (et[i]->get_u_euler()) ;
544 
545  const Metric& gtilde = et[i]->get_gtilde() ;
546  const Metric_flat flat (mp.flat_met_cart()) ;
547 
548  Vector vphi_cov = vphi.up_down(gtilde) ;
549 
550  Scalar integrand = contract(jmom, 0, vphi_cov, 0) ;
551 
552  p_angu_mom->set(2) += integrand.integrale() ;
553 
554  // Extrinsic curvature part (0 if IWM)
555  // -----------------------------------
556 
557  const Sym_tensor& aij = et[i]->get_tkij_auto() ;
558  rho = pow(psi4, double(1.5)) ;
559  rho.std_spectral_base() ;
560 
561  // Construction of D_k \Phi^i
562  Itbl type (2) ;
563  type.set(0) = CON ;
564  type.set(1) = COV ;
565 
566  Tensor dcov_vphi (mp, 2, type, mp.get_bvect_cart()) ;
567  dcov_vphi.set(1,1) = 0. ;
568  dcov_vphi.set(2,1) = ff ;
569  dcov_vphi.set(3,1) = 0. ;
570  dcov_vphi.set(2,2) = 0. ;
571  dcov_vphi.set(3,2) = 0. ;
572  dcov_vphi.set(3,3) = 0. ;
573  dcov_vphi.set(1,2) = -ff ;
574  dcov_vphi.set(1,3) = 0. ;
575  dcov_vphi.set(2,3) = 0. ;
576  dcov_vphi.inc_dzpuis(2) ;
577 
578  Connection gamijk (gtilde, flat) ;
579  const Tensor& deltaijk = gamijk.get_delta() ;
580 
581  // Computation of \tilde D_i \tilde \Phi_j
582  Sym_tensor kill_phi (mp, COV, mp.get_bvect_cart()) ;
583  kill_phi = contract(gtilde.cov(), 1, dcov_vphi +
584  contract(deltaijk, 2, vphi, 0), 0) +
585  contract(dcov_vphi + contract(deltaijk, 2, vphi, 0), 0,
586  gtilde.cov(), 1) ;
587 
588  integrand = rho * contract(aij, 0, 1, kill_phi, 0, 1) ;
589 
590  p_angu_mom->set(2) += integrand.integrale() / (4*qpig) ;
591 
592 
593  } // End of the loop on the stars
594 
595  } // End of the case where a new computation was necessary
596 
597  return *p_angu_mom ;
598 
599 }
600 
601 
602 
603  //---------------------------------//
604  // Total energy //
605  //---------------------------------//
606 
607 double Binary::total_ener() const {
608  /*
609  if (p_total_ener == 0x0) { // a new computation is requireed
610 
611  p_total_ener = new double ;
612 
613  *p_total_ener = mass_adm() - star1.mass_b() - star2.mass_b() ;
614 
615  } // End of the case where a new computation was necessary
616 
617  */
618  return *p_total_ener ;
619 
620 }
621 
622 
623  //---------------------------------//
624  // Error on the virial theorem //
625  //---------------------------------//
626 
627 double Binary::virial() const {
628 
629  if (p_virial == 0x0) { // a new computation is requireed
630 
631  p_virial = new double ;
632 
633  *p_virial = 1. - mass_kom() / mass_adm() ;
634 
635  }
636 
637  return *p_virial ;
638 
639 }
640 }
Cartesian vectorial bases (triads).
Definition: base_vect.h:201
double mass_adm_vol() const
Total ADM mass (computed by a volume integral)
double mass_kom_vol() const
Total Komar mass (computed by a volume integral)
Star_bin * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
Definition: binary.h:89
double * p_total_ener
Total energy of the system.
Definition: binary.h:114
Tbl * p_angu_mom
Total angular momentum of the system.
Definition: binary.h:111
double total_ener() const
Total energy (excluding the rest mass energy).
const Tbl & angu_mom() const
Total angular momentum.
double * p_mass_adm
Total ADM mass of the system.
Definition: binary.h:105
double * p_virial
Virial theorem error.
Definition: binary.h:117
double mass_adm() const
Total ADM mass.
double virial() const
Estimates the relative error on the virial theorem.
double * p_mass_kom
Total Komar mass of the system.
Definition: binary.h:108
double mass_kom() const
Total Komar mass.
Class Connection.
Definition: connection.h:113
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition: connection.h:271
Basic integer array class.
Definition: itbl.h:122
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Affine radial mapping.
Definition: map.h:2027
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Base class for coordinate mappings.
Definition: map.h:670
Coord ya
Absolute y coordinate.
Definition: map.h:731
Coord r
r coordinate centered on the grid
Definition: map.h:718
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 Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:331
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
Coord xa
Absolute x coordinate.
Definition: map.h:730
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Flat metric for tensor calculation.
Definition: metric.h:261
Metric for tensor calculation.
Definition: metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:61
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition: scalar_manip.C:315
const Vector & get_dcov_phi() const
Returns the covariant derivative of (logarithm of the conformal factor).
Definition: star.h:846
const Metric & get_flat() const
Return the flat metric defined on the mapping (Spherical components with respect to the mapping of th...
Definition: star.h:859
const Scalar & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: star.h:811
const Vector & get_dcov_logn() const
Returns the covariant derivative of .
Definition: star.h:837
const Metric & get_gtilde() const
Return the conformal 3-metric .
Definition: star.h:862
const Vector & get_dcon_logn() const
Returns the contravariant derivative of .
Definition: star.h:841
const Scalar & get_logn_auto() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: star.h:806
const Sym_tensor & get_tkij_auto() const
Returns the part of the extrinsic curvature tensor generated by beta_auto.
Definition: star.h:886
const Scalar & get_kcar_auto() const
Returns the part of generated by beta_auto.
Definition: star.h:897
const Scalar & get_lnq_auto() const
Returns the part of the vector field generated principally by the star.
Definition: star.h:828
const Sym_tensor & get_hij() const
Return the total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:867
const Sym_tensor & get_hij_auto() const
Return the deviation of the inverse conformal metric from the inverse flat metric principally genera...
Definition: star.h:873
const Scalar & get_kcar_comp() const
Returns the part of generated by beta_comp.
Definition: star.h:902
const Scalar & get_psi4() const
Return the conformal factor .
Definition: star.h:854
const Scalar & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Definition: star.h:379
const Scalar & get_press() const
Returns the fluid pressure.
Definition: star.h:373
const Vector & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:385
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: star.h:376
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
const Scalar & get_logn() const
Returns the logarithm of the lapse N.
Definition: star.h:396
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
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.
Definition: tensor.h:288
Tensor field of valence 1.
Definition: vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:816
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.