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
29char 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
102namespace Lorene {
103double 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
135double 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
201double 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
229double 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
285const 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
607double 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
627double 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
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 ya
Absolute y coordinate.
Definition map.h:731
Coord r
r coordinate centered on the grid
Definition map.h:718
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
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
double integrale() const
Computes the integral over all space of *this .
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.
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 Scalar & get_lnq_auto() const
Returns the part of the vector field generated principally by the star.
Definition star.h:828
const Vector & get_dcov_logn() const
Returns the covariant derivative of .
Definition star.h:837
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 Sym_tensor & get_hij() const
Return the total deviation of the inverse conformal metric from the inverse flat metric.
Definition star.h:867
const Scalar & get_kcar_auto() const
Returns the part of generated by beta_auto.
Definition star.h:897
const Vector & get_dcon_logn() const
Returns the contravariant derivative of .
Definition star.h:841
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
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 Sym_tensor & get_tkij_auto() const
Returns the part of the extrinsic curvature tensor generated by beta_auto.
Definition star.h:886
const Metric & get_gtilde() const
Return the conformal 3-metric .
Definition star.h:862
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_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition star.h:385
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
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition star.h:376
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
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.