23char bhole_glob_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_glob.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $" ;
82#include "utilitaires.h"
115 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
117 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
118 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
131 for (
int i=0 ; i<3 ; i++)
154 for (
int i=0 ; i<3 ; i++)
166 return moment_un+same_orient*moment_deux ;
177 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
179 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
180 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
184 double* bornes =
new double [nzones+1] ;
186 for (
int i=nzones-1 ; i>0 ; i--) {
187 bornes[i] = courant ;
191 bornes[nzones] = __infinity ;
215 Tenseur shift_tot (shift_un+shift_deux) ;
217 shift_tot.
annule(0, nzones-2) ;
223 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
224 for (
int i=0 ; i<3 ; i++) {
225 k_total.
set(i, i) = grad(i, i)-trace()/3. ;
226 for (
int j=i+1 ; j<3 ; j++)
227 k_total.
set(i, j) = (grad(i, j)+grad(j, i))/2. ;
230 for (
int lig=0 ; lig<3 ; lig++)
231 for (
int col=lig ; col<3 ; col++)
236 for (
int i=0 ; i<3 ; i++)
237 vecteur_un.
set(i) = k_total(0, i) ;
239 Cmp integrant_un (vecteur_un(0)) ;
243 for (
int i=0 ; i<3 ; i++)
244 vecteur_deux.
set(i) = k_total(1, i) ;
246 Cmp integrant_deux (vecteur_deux(0)) ;
269 double pente = 2./(x_un-x_deux) ;
270 double constante = - (x_un+x_deux)/(x_un-x_deux) ;
275 double air_un, theta_un, phi_un ;
276 double air_deux, theta_deux, phi_deux ;
278 double* coloc =
new double[nr] ;
279 double* coef =
new double[nr] ;
280 int* deg =
new int[3] ;
281 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
283 for (
int i=0 ; i<nr ; i++) {
284 ksi = -
cos (M_PI*i/(nr-1)) ;
285 xabs = (ksi-constante)/pente ;
291 hole2.
psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
295 cfrcheb(deg, deg, coloc, deg, coef) ;
298 double* som =
new double[nr] ;
300 for (
int i=2 ; i<nr ; i+=2)
301 som[i] = 1./(i+1)-1./(i-1) ;
302 for (
int i=1 ; i<nr ; i+=2)
306 for (
int i=0 ; i<nr ; i++)
307 res += som[i]*coef[i] ;
319Tbl Bhole_binaire::linear_momentum_systeme_inf()
const {
326 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
328 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
329 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
333 double* bornes =
new double [nzones+1] ;
335 for (
int i=nzones-1 ; i>0 ; i--) {
336 bornes[i] = courant ;
340 bornes[nzones] = __infinity ;
347 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
348 k_total.set_etat_qcq() ;
350 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
351 shift_un.set_etat_qcq() ;
353 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
354 shift_deux.set_etat_qcq() ;
364 Tenseur shift_tot (shift_un+shift_deux) ;
365 shift_tot.set_std_base() ;
366 shift_tot.annule(0, nzones-2) ;
369 Tenseur shift_old (shift_tot) ;
370 shift_tot.inc2_dzpuis() ;
371 shift_tot.dec2_dzpuis() ;
372 for (
int i=0 ; i< 3 ; i++)
376 Tenseur grad (shift_tot.gradient()) ;
377 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
378 for (
int i=0 ; i<3 ; i++) {
379 k_total.set(i, i) = grad(i, i)-trace()/3. ;
380 for (
int j=i+1 ; j<3 ; j++)
381 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
385 for (
int lig=0 ; lig<3 ; lig++)
386 for (
int col=lig ; col<3 ; col++)
387 k_total.set(lig, col).mult_r_zec() ;
389 for (
int comp=0 ; comp<3 ; comp++) {
390 Tenseur vecteur (mapping, 1, CON, mapping.get_bvect_cart()) ;
391 vecteur.set_etat_qcq() ;
392 for (
int i=0 ; i<3 ; i++)
393 vecteur.set(i) = k_total(i, comp) ;
394 vecteur.change_triad (mapping.get_bvect_spher()) ;
395 Cmp integrant (vecteur(0)) ;
397 res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
double moment_systeme_inf()
Calculates the angular momentum of the black hole using the formula at infinity : where .
double omega
Position of the axis of rotation.
Bhole hole1
Black hole one.
Bhole hole2
Black hole two.
double komar_systeme() const
Calculates the Komar mass of the system using : .
double adm_systeme() const
Calculates the ADM mass of the system using : .
double moment_systeme_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and den...
double distance_propre(const int nr=65) const
Calculation of the proper distance between the two spheres of inversion, along the x axis.
Tenseur psi_auto
Part of generated by the hole.
Tenseur shift_auto
Part of generated by the hole.
Tenseur n_auto
Part of N generated by the hole.
double rayon
Radius of the horizon in LORENE's units.
Map_af & mp
Affine mapping.
Tenseur_sym tkij_tot
Total .
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Valeur va
The numerical value of the Cmp
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void import_symy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
void import_asymy(const Cmp &ci)
Assignment to another Cmp defined on a different mapping.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord ya
Absolute y coordinate.
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X,...
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
double get_ori_x() const
Returns the x coordinate of the origin.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Coord xa
Absolute x coordinate.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_nzone() const
Returns the number of domains.
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void annule(int l)
Sets the Tenseur to zero in a given domain.
void dec2_dzpuis()
dzpuis -= 2 ;
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void inc2_dzpuis()
dzpuis += 2 ;
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
const Valeur & mult_st() const
Returns applied to *this.
const Valeur & mult_cp() const
Returns applied to *this.
const Valeur & mult_sp() const
Returns applied to *this.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).