binhor_global.C
/*
* Copyright (c) 2005 Francois Limousin
* Jose Luis Jaramillo
*
* 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 binhor_glob_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_global.C,v 1.10 2014/10/13 08:52:42 j_novak Exp $" ;
/*
* $Id: binhor_global.C,v 1.10 2014/10/13 08:52:42 j_novak Exp $
* $Log: binhor_global.C,v $
* Revision 1.10 2014/10/13 08:52:42 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.9 2014/10/06 15:13:01 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.8 2007/04/13 15:28:55 f_limousin
* Lots of improvements, generalisation to an arbitrary state of
* rotation, implementation of the spatial metric given by Samaya.
*
* Revision 1.7 2006/05/24 16:56:37 f_limousin
* Many small modifs.
*
* Revision 1.6 2005/09/24 08:31:38 f_limousin
* Improve the computation of moment_adm() and moment_hor().
*
* Revision 1.5 2005/09/13 18:33:15 f_limousin
* New function vv_bound_cart_bin(double) for computing binaries with
* berlin condition for the shift vector.
* Suppress all the symy and asymy in the importations.
*
* Revision 1.4 2005/06/09 16:12:04 f_limousin
* Implementation of amg_mom_adm().
*
* Revision 1.3 2005/04/29 14:02:44 f_limousin
* Important changes : manage the dependances between quantities (for
* instance psi and psi4). New function write_global(ost).
*
* Revision 1.2 2005/03/04 17:09:57 jl_jaramillo
* Change to avoid warnings
*
* Revision 1.1 2005/03/03 13:48:56 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_global.C,v 1.10 2014/10/13 08:52:42 j_novak Exp $
*
*/
//standard
#include <cstdlib>
#include <cmath>
// Lorene
#include "nbr_spx.h"
#include "tensor.h"
#include "isol_hor.h"
#include "proto.h"
#include "utilitaires.h"
//#include "graphique.h"
namespace Lorene {
double Bin_hor::adm_mass() const {
Vector dpsi_un (hole1.psi_auto.derive_con(hole1.ff)) ;
Vector dpsi_deux (hole2.psi_auto.derive_con(hole2.ff)) ;
Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
Vector ww (0.125*(hdirac1 - (hole1.hh.trace(hole1.ff)).
derive_con(hole1.ff))) ;
double inf = hole1.mp.val_r(hole1.mp.get_mg()->get_nzone()-1, 1., 0., 0.) ;
double masse = dpsi_un.flux(inf, hole1.ff) +
dpsi_deux.flux(inf, hole2.ff) +
ww.flux(inf, hole1.ff) ;
masse /= -2*M_PI ;
return masse ;
}
double Bin_hor::komar_mass() const {
Vector dnn_un (hole1.n_auto.derive_con(hole1.ff)) ;
Vector dnn_deux (hole2.n_auto.derive_con(hole2.ff)) ;
Vector ww (contract(hole1.hh, 1, hole1.nn.derive_cov(hole1.ff), 0)) ;
double inf = hole1.mp.val_r(hole1.mp.get_mg()->get_nzone()-1, 1., 0., 0.) ;
double mass = dnn_un.flux(inf, hole1.ff) +
dnn_deux.flux(inf, hole2.ff) +
ww.flux(inf, hole1.ff) ;
mass /= 4*M_PI ;
return mass ;
}
double Bin_hor::ang_mom_hor() const {
if (omega == 0)
return 0 ;
else {
// Alignement
double orientation_un = hole1.mp.get_rot_phi() ;
assert ((orientation_un==0) || (orientation_un==M_PI)) ;
double orientation_deux = hole2.mp.get_rot_phi() ;
assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
// Integral on the first horizon :
Scalar xa_un (hole1.mp) ;
xa_un = hole1.mp.xa ;
xa_un.std_spectral_base() ;
Scalar ya_un (hole1.mp) ;
ya_un = hole1.mp.ya ;
ya_un.std_spectral_base() ;
Sym_tensor tkij_un (hole1.aa) ;
tkij_un.change_triad(hole1.mp.get_bvect_cart()) ;
Vector vecteur_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++)
vecteur_un.set(i) = -ya_un*tkij_un(1, i)+
xa_un * tkij_un(2, i) ;
vecteur_un.annule_domain(hole1.mp.get_mg()->get_nzone()-1) ;
vecteur_un.change_triad (hole1.mp.get_bvect_spher()) ;
Scalar integrant_un (hole1.get_psi4()*hole1.psi*hole1.psi
*vecteur_un(1)) ;
double moment_un = hole1.mp.integrale_surface
(integrant_un, hole1.radius+1e-12)/8/M_PI ;
//Integral on the second horizon :
Scalar xa_deux (hole2.mp) ;
xa_deux = hole2.mp.xa ;
xa_deux.std_spectral_base() ;
Scalar ya_deux (hole2.mp) ;
ya_deux = hole2.mp.ya ;
ya_deux.std_spectral_base() ;
Sym_tensor tkij_deux (hole2.aa) ;
tkij_deux.change_triad(hole2.mp.get_bvect_cart()) ;
Vector vecteur_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++)
vecteur_deux.set(i) = -ya_deux*tkij_deux(1, i)+
xa_deux * tkij_deux(2, i) ;
vecteur_deux.annule_domain(hole2.mp.get_mg()->get_nzone()-1) ;
vecteur_deux.change_triad (hole2.mp.get_bvect_spher()) ;
Scalar integrant_deux (hole2.get_psi4()*hole2.psi*hole2.psi
*vecteur_deux(1)) ;
double moment_deux = hole2.mp.integrale_surface
(integrant_deux, hole2.radius+1e-12)/8/M_PI ;
return moment_un+same_orient*moment_deux ;
}
}
double Bin_hor::ang_mom_adm() const {
/*
Scalar integrand_un (hole1.k_dd()(1,3) - hole1.gam_dd()(1,3) * hole1.trK) ;
integrand_un.mult_rsint() ; // in order to pass from the triad components
double mom = hole1.mp.integrale_surface_infini(integrand_un) ;
mom /= 8*M_PI ;
return mom ;
*/
if (omega == 0)
return 0 ;
else {
// Alignement
double orientation_un = hole1.mp.get_rot_phi() ;
assert ((orientation_un==0) || (orientation_un==M_PI)) ;
double orientation_deux = hole2.mp.get_rot_phi() ;
assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
// Construction of an auxiliar grid and mapping
int nzones = hole1.mp.get_mg()->get_nzone() ;
double* bornes = new double [nzones+1] ;
double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
for (int i=nzones-1 ; i>0 ; i--) {
bornes[i] = courant ;
courant /= 2. ;
}
bornes[0] = 0 ;
bornes[nzones] = __infinity ;
Map_af mapping (*hole1.mp.get_mg(), bornes) ;
delete [] bornes ;
// Construction of k_total
Sym_tensor k_total (mapping, CON, mapping.get_bvect_cart()) ;
Vector shift_un (mapping, CON, mapping.get_bvect_cart()) ;
Vector shift_deux (mapping, CON, mapping.get_bvect_cart()) ;
Vector beta_un (hole1.beta_auto) ;
Vector beta_deux (hole2.beta_auto) ;
beta_un.change_triad(hole1.mp.get_bvect_cart()) ;
beta_deux.change_triad(hole2.mp.get_bvect_cart()) ;
shift_un.set(1).import(beta_un(1)) ;
shift_un.set(2).import(beta_un(2)) ;
shift_un.set(3).import(beta_un(3)) ;
shift_deux.set(1).import(same_orient*beta_deux(1)) ;
shift_deux.set(2).import(same_orient*beta_deux(2)) ;
shift_deux.set(3).import(beta_deux(3)) ;
Vector shift_tot (shift_un+shift_deux) ;
shift_tot.std_spectral_base() ;
shift_tot.annule(0, nzones-2) ;
// Substract the residuals
shift_tot.inc_dzpuis(2) ;
shift_tot.dec_dzpuis(2) ;
const Metric_flat& flat0 (mapping.flat_met_cart()) ;
k_total = shift_tot.ope_killing_conf(flat0) / 2. ;
//- flat0.con() * hole1.trK;
for (int lig=1 ; lig<=3 ; lig++)
for (int col=lig ; col<=3 ; col++)
k_total.set(lig, col).mult_r_ced() ;
Vector vecteur_un (mapping, CON, mapping.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++)
vecteur_un.set(i) = k_total(1, i) ;
vecteur_un.change_triad (mapping.get_bvect_spher()) ;
Scalar integrant_un (vecteur_un(1)) ;
Vector vecteur_deux (mapping, CON, mapping.get_bvect_cart()) ;
for (int i=1 ; i<=3 ; i++)
vecteur_deux.set(i) = k_total(2, i) ;
vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
Scalar integrant_deux (vecteur_deux(1)) ;
// Multiplication by y and x :
integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
.mult_st() ;
integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
.mult_sp() ;
integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
.mult_st() ;
integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
.mult_cp() ;
double moment = mapping.integrale_surface_infini (-integrant_un
+integrant_deux) ;
moment /= 8*M_PI ;
return moment ;
}
}
/*
double Bin_hor::proper_distance(const int nr) const {
// On determine les rayons coordonnes des points limites de l'integrale :
double x_un = hole1.mp.get_ori_x() - hole1.rayon ;
double x_deux = hole2.mp.get_ori_x() + hole2.rayon ;
// Les coefficients du changement de variable :
double pente = 2./(x_un-x_deux) ;
double constante = - (x_un+x_deux)/(x_un-x_deux) ;
double ksi ; // variable d'integration.
double xabs ; // x reel.
double air_un, theta_un, phi_un ; // coordonnee spheriques 1
double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
double* coloc = new double[nr] ;
double* coef = new double[nr] ;
int* deg = new int[3] ;
deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
for (int i=0 ; i<nr ; i++) {
ksi = -cos (M_PI*i/(nr-1)) ;
xabs = (ksi-constante)/pente ;
hole1.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
hole2.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
coloc[i] = pow (hole1.psi_auto().val_point (air_un, theta_un, phi_un) +
hole2.psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
}
// On prend les coefficients de la fonction
cfrcheb(deg, deg, coloc, deg, coef) ;
// On integre
double* som = new double[nr] ;
som[0] = 2 ;
for (int i=2 ; i<nr ; i+=2)
som[i] = 1./(i+1)-1./(i-1) ;
for (int i=1 ; i<nr ; i+=2)
som[i] = 0 ;
double res = 0 ;
for (int i=0 ; i<nr ; i++)
res += som[i]*coef[i] ;
res /= pente ;
delete [] deg ;
delete [] coef ;
delete [] coloc ;
delete [] som ;
return res ;
}
*/
}
Lorene::Bin_hor::omega
double omega
Angular velocity.
Definition isol_hor.h:1348
Lorene::Bin_hor::hole1
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Lorene::Bin_hor::hole2
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
Lorene::Bin_hor::ang_mom_hor
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
Definition binhor_global.C:118
Lorene::Bin_hor::adm_mass
double adm_mass() const
Calculates the ADM mass of the system.
Definition binhor_global.C:83
Lorene::Bin_hor::ang_mom_adm
double ang_mom_adm() const
Calculates the angular momentum of the black hole.
Definition binhor_global.C:184
Lorene::Bin_hor::komar_mass
double komar_mass() const
Calculates the Komar mass of the system using : .
Definition binhor_global.C:101
Lorene::Map_af
Affine radial mapping.
Definition map.h:2027
Lorene::Map_af::integrale_surface_infini
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Definition map_af_integ_surf.C:196
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::std_spectral_base
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
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::mult_r_ced
void mult_r_ced()
Multiplication by r in the compactified external domain (CED), the dzpuis flag is not changed.
Definition scalar_r_manip.C:269
Lorene::Scalar::import
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition scalar_import.C:68
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur::mult_st
const Valeur & mult_st() const
Returns applied to *this.
Definition valeur_mult_st.C:97
Lorene::Valeur::mult_cp
const Valeur & mult_cp() const
Returns applied to *this.
Definition valeur_mult_cp.C:77
Lorene::Valeur::mult_sp
const Valeur & mult_sp() const
Returns applied to *this.
Definition valeur_mult_sp.C:77
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::flux
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition vector.C:807
Lorene::Vector::ope_killing_conf
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:75
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:77
Lorene::Tensor::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
Lorene::Tensor::inc_dzpuis
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
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64