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