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
30char 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
87namespace Lorene {
88double 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
136double 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
198const 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
289double 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
345double 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
408double 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
481double 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}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Binaire::p_angu_mom
Tbl * p_angu_mom
Total angular momentum of the system.
Definition binaire.h:145
Lorene::Binaire::p_virial_gb
double * p_virial_gb
Virial theorem error by E.Gourgoulhon and S.Bonazzola.
Definition binaire.h:154
Lorene::Binaire::mass_adm
double mass_adm() const
Total ADM mass.
Definition binaire_global.C:88
Lorene::Binaire::virial_gb
double virial_gb() const
Estimates the relative error on the virial theorem calculated by E.Gourgoulhon and S....
Definition binaire_global.C:408
Lorene::Binaire::et
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
Lorene::Binaire::angu_mom
const Tbl & angu_mom() const
Total angular momentum.
Definition binaire_global.C:198
Lorene::Binaire::p_mass_adm
double * p_mass_adm
Total ADM mass of the system.
Definition binaire.h:139
Lorene::Binaire::p_total_ener
double * p_total_ener
Total energy of the system.
Definition binaire.h:148
Lorene::Binaire::p_virial
double * p_virial
Virial theorem error.
Definition binaire.h:151
Lorene::Binaire::star2
Etoile_bin star2
Second star of the system.
Definition binaire.h:118
Lorene::Binaire::mass_kom
double mass_kom() const
Total Komar mass.
Definition binaire_global.C:136
Lorene::Binaire::p_mass_kom
double * p_mass_kom
Total Komar mass of the system.
Definition binaire.h:142
Lorene::Binaire::virial
double virial() const
Estimates the relative error on the virial theorem (for a relativistic one, it returns $|1 - M_{\rm K...
Definition binaire_global.C:345
Lorene::Binaire::p_virial_fus
double * p_virial_fus
Virial theorem error by J.L.Friedman, K.Uryu, and M.Shibata.
Definition binaire.h:157
Lorene::Binaire::star1
Etoile_bin star1
First star of the system.
Definition binaire.h:115
Lorene::Binaire::total_ener
double total_ener() const
Total energy (excluding the rest mass energy).
Definition binaire_global.C:289
Lorene::Binaire::virial_fus
double virial_fus() const
Estimates the relative error on the virial theorem calculated by J.L.Friedman, K.Uryu,...
Definition binaire_global.C:481
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Lorene::Cmp::integrale
double integrale() const
Computes the integral over all space of *this .
Definition cmp_integ.C:55
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::flat_scalar_prod
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...
Definition tenseur_operateur.C:653
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Unites
Standard units of space, time and mass.