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