binaire_global.C
/*
* Methods of class Binaire to compute global quantities
*
* (see file binaire.h for documentation)
*/
/*
* Copyright (c) 2000-2001 Eric Gourgoulhon
* Copyright (c) 2000-2001 Keisuke Taniguchi
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char binaire_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_global.C,v 1.7 2014/10/13 08:52:44 j_novak Exp $" ;
/*
* $Id: binaire_global.C,v 1.7 2014/10/13 08:52:44 j_novak Exp $
* $Log: binaire_global.C,v $
* Revision 1.7 2014/10/13 08:52:44 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2004/03/25 10:28:59 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.5 2003/09/08 09:32:40 e_gourgoulhon
* Corrected a problem of spectral basis initialisation in virial_gb() and
* virial_fus(): introduced the new variable a1.
*
* Revision 1.4 2001/12/20 14:18:40 k_taniguchi
* Addition of the Komar mass, the virial error by Gourgoulhon and Bonazzola, and the virial error by Friedman, Uryu, and Shibata.
*
* Revision 1.3 2001/12/16 16:21:38 e_gourgoulhon
* #include "unites.h" is now local to Binaire::mass_adm(), in order
* not to make Lorene's units global variables.
*
* Revision 1.2 2001/12/14 09:45:14 k_taniguchi
* Correction of missing 16 pi G factor in the ADM mass
*
* Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
* LORENE
*
* Revision 2.3 2000/03/08 12:26:33 eric
* Ajout de l'appel a std_base_scal() sur le Cmp source dans le cas
* relativiste (masse ADM).
*
* Revision 2.2 2000/02/23 11:26:00 keisuke
* Changement de "virial relation".
*
* Revision 2.1 2000/02/18 15:48:55 eric
* *** empty log message ***
*
* Revision 2.0 2000/02/18 14:53:09 eric
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_global.C,v 1.7 2014/10/13 08:52:44 j_novak Exp $
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "binaire.h"
#include "unites.h"
//---------------------------------//
// ADM mass //
//---------------------------------//
namespace Lorene {
double Binaire::mass_adm() const {
using namespace Unites ;
if (p_mass_adm == 0x0) { // a new computation is requireed
p_mass_adm = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_mass_adm = 0 ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp& a2 = (et[i]->get_a_car())() ;
const Cmp& ee = (et[i]->get_ener_euler())() ;
const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
Cmp source = pow(a2, 1.25) * ee
+ pow(a2, 0.25) * (ak2_auto + ak2_comp) / (4.*qpig) ;
source.std_base_scal() ;
*p_mass_adm += source.integrale() ;
}
}
else { // Newtonian case
// --------------
*p_mass_adm = star1.mass_b() + star2.mass_b() ;
}
} // End of the case where a new computation was necessary
return *p_mass_adm ;
}
//---------------------------------//
// Komar mass //
//---------------------------------//
double Binaire::mass_kom() const {
using namespace Unites ;
if (p_mass_kom == 0x0) { // a new computation is requireed
p_mass_kom = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_mass_kom = 0 ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp& lapse = (et[i]->get_nnn())() ;
const Cmp& a2 = (et[i]->get_a_car())() ;
const Cmp& ee = (et[i]->get_ener_euler())() ;
const Cmp& se = (et[i]->get_s_euler())() ;
const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
Tenseur dndb = flat_scalar_prod(dnu_auto,
dbe_auto + dbe_comp) ;
Tenseur dndn = flat_scalar_prod(dnu_auto,
dnu_auto + dnu_comp) ;
Cmp source = lapse * ( a2 * (ee + se)
+ (ak2_auto + ak2_comp)/qpig
- dndb()/qpig + dndn()/qpig ) ;
source.std_base_scal() ;
*p_mass_kom += source.integrale() ;
}
}
else { // Newtonian case
// --------------
*p_mass_kom = star1.mass_b() + star2.mass_b() ;
}
} // End of the case where a new computation was necessary
return *p_mass_kom ;
}
//---------------------------------//
// Total angular momentum //
//---------------------------------//
const Tbl& Binaire::angu_mom() const {
if (p_angu_mom == 0x0) { // a new computation is requireed
p_angu_mom = new Tbl(3) ;
p_angu_mom->annule_hard() ; // fills the double array with zeros
for (int i=0; i<=1; i++) { // loop on the stars
const Map& mp = et[i]->get_mp() ;
Cmp xx(mp) ;
Cmp yy(mp) ;
Cmp zz(mp) ;
xx = mp.xa ;
yy = mp.ya ;
zz = mp.za ;
const Cmp& vx = (et[i]->get_u_euler())(0) ;
const Cmp& vy = (et[i]->get_u_euler())(1) ;
const Cmp& vz = (et[i]->get_u_euler())(2) ;
Cmp rho(mp) ;
if ( et[i]->is_relativistic() ) {
const Cmp& a2 = (et[i]->get_a_car())() ;
const Cmp& ee = (et[i]->get_ener_euler())() ;
const Cmp& pp = (et[i]->get_press())() ;
rho = pow(a2, 2.5) * (ee + pp) ;
}
else {
rho = (et[i]->get_nbar())() ;
}
Base_val** base = (et[i]->get_mp()).get_mg()->std_base_vect_cart() ;
// X component
// -----------
Cmp source = rho * ( yy * vz - zz * vy ) ;
(source.va).set_base( *(base[2]) ) ; // same basis as V^z
//## p_angu_mom->set(0) += source.integrale() ;
p_angu_mom->set(0) += 0 ;
// y component
// -----------
source = rho * ( zz * vx - xx * vz ) ;
(source.va).set_base( *(base[2]) ) ; // same basis as V^z
//## p_angu_mom->set(1) += source.integrale() ;
p_angu_mom->set(1) += 0 ;
// Z component
// -----------
source = rho * ( xx * vy - yy * vx ) ;
source.std_base_scal() ; // same basis as V^x (standard scalar
// field)
p_angu_mom->set(2) += source.integrale() ;
delete base[0] ;
delete base[1] ;
delete base[2] ;
delete [] base ;
} // End of the loop on the stars
} // End of the case where a new computation was necessary
return *p_angu_mom ;
}
//---------------------------------//
// Total energy //
//---------------------------------//
double Binaire::total_ener() const {
if (p_total_ener == 0x0) { // a new computation is requireed
p_total_ener = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_total_ener = mass_adm() - star1.mass_b() - star2.mass_b() ;
}
else { // Newtonian case
// --------------
*p_total_ener = 0 ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp e_int = (et[i]->get_ener())()
- (et[i]->get_nbar())() ;
const Cmp& rho = (et[i]->get_nbar())() ;
// Fluid velocity with respect to the inertial frame
const Tenseur& vit = et[i]->get_u_euler() ;
Cmp vit2 = flat_scalar_prod(vit, vit)() ;
// Gravitational potential
const Cmp nu = (et[i]->get_logn_auto())()
+ (et[i]->get_logn_comp())() ;
Cmp source = e_int + .5 * rho * vit2 + .5 * rho * nu ;
*p_total_ener += source.integrale() ;
} // End of the loop on the stars
} // End of Newtonian case
} // End of the case where a new computation was necessary
return *p_total_ener ;
}
//---------------------------------//
// Error on the virial theorem //
//---------------------------------//
double Binaire::virial() const {
if (p_virial == 0x0) { // a new computation is requireed
p_virial = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_virial = 1. - mass_kom() / mass_adm() ;
}
else { // Newtonian case
// --------------
*p_virial = 0 ;
double vir_mat = 0 ;
double vir_grav = 0 ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp& pp = (et[i]->get_press())() ;
const Cmp& rho = (et[i]->get_nbar())() ;
// Fluid velocity with respect to the inertial frame
const Tenseur& vit = et[i]->get_u_euler() ;
Cmp vit2 = flat_scalar_prod(vit, vit)() ;
// Gravitational potential
const Cmp nu = (et[i]->get_logn_auto())()
+ (et[i]->get_logn_comp())() ;
Cmp source = 3*pp + rho * vit2 ;
vir_mat += source.integrale() ;
source = .5 * rho * nu ;
vir_grav += source.integrale() ;
} // End of the loop on the stars
*p_virial = ( vir_mat + vir_grav ) / fabs(vir_grav) ;
} // End of the Newtonian case
} // End of the case where a new computation was necessary
return *p_virial ;
}
//----------------------------------------------//
// Virial error by Gourgoulhon and Bonazzola //
//----------------------------------------------//
double Binaire::virial_gb() const {
using namespace Unites ;
if (p_virial_gb == 0x0) { // a new computation is requireed
p_virial_gb = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_virial_gb = 0 ;
double vir_pres = 0. ;
double vir_extr = 0. ;
double vir_grav = 0. ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp& a2 = (et[i]->get_a_car())() ;
const Cmp& se = (et[i]->get_s_euler())() ;
const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
Cmp a1 = sqrt(a2) ;
a1.std_base_scal() ;
Cmp source = 2. * a2 * a1 * se ;
vir_pres += source.integrale() ;
source = 1.5 * a1 * (ak2_auto + ak2_comp) / qpig ;
source.std_base_scal() ;
vir_extr += source.integrale() ;
Tenseur sprod1 = flat_scalar_prod(dbe_auto, dbe_auto+dbe_comp) ;
Tenseur sprod2 = flat_scalar_prod(dnu_auto, dnu_auto+dnu_comp) ;
Tenseur sprod3 = flat_scalar_prod(dbe_auto, dnu_auto+dnu_comp) ;
source = a1 * ( sprod1() - sprod2() - 2.*sprod3() )/qpig ;
vir_grav += source.integrale() ;
} // End of the loop on the stars
*p_virial_gb = (vir_pres + vir_extr + vir_grav) / mass_adm() ;
}
else { // Newtonian case
// --------------
*p_virial_gb = virial() ;
} // End of the Newtonian case
} // End of the case where a new computation was necessary
return *p_virial_gb ;
}
//------------------------------------------------//
// Virial error by Friedman, Uryu, and Shibata //
//------------------------------------------------//
double Binaire::virial_fus() const {
using namespace Unites ;
if (p_virial_fus == 0x0) { // a new computation is requireed
p_virial_fus = new double ;
if (star1.is_relativistic()) { // Relativistic case
// -----------------
assert( star2.is_relativistic() ) ;
*p_virial_fus = 0 ;
double vir_pres = 0. ;
double vir_extr = 0. ;
double vir_grav = 0. ;
for (int i=0; i<=1; i++) { // loop on the stars
const Cmp& lapse = (et[i]->get_nnn())() ;
const Cmp& a2 = (et[i]->get_a_car())() ;
const Cmp& se = (et[i]->get_s_euler())() ;
const Cmp& ak2_auto = (et[i]->get_akcar_auto())() ;
const Cmp& ak2_comp = (et[i]->get_akcar_comp())() ;
const Tenseur& dnu_auto = et[i]->get_d_logn_auto() ;
const Tenseur& dnu_comp = et[i]->get_d_logn_comp() ;
const Tenseur& dbe_auto = et[i]->get_d_beta_auto() ;
const Tenseur& dbe_comp = et[i]->get_d_beta_comp() ;
Cmp a1 = sqrt(a2) ;
a1.std_base_scal() ;
Cmp source = 2. * lapse * a2 * a1 * se ;
vir_pres += source.integrale() ;
source = 1.5 * lapse * a1 * (ak2_auto + ak2_comp) / qpig ;
vir_extr += source.integrale() ;
Tenseur sprod = flat_scalar_prod( dbe_auto, dbe_auto+dbe_comp )
- flat_scalar_prod( dnu_auto, dnu_auto+dnu_comp ) ;
source = lapse * a1 * sprod() / qpig ;
vir_grav += source.integrale() ;
} // End of the loop on the stars
*p_virial_fus = (vir_pres + vir_extr + vir_grav) / mass_adm() ;
}
else { // Newtonian case
// --------------
*p_virial_fus = virial() ;
} // End of the Newtonian case
} // End of the case where a new computation was necessary
return *p_virial_fus ;
}
}
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.