85 #include "utilitaires.h" 86 #include "graphique.h" 118 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
120 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
121 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
134 for (
int i=0 ; i<3 ; i++)
137 vecteur_un.set_std_base() ;
157 for (
int i=0 ; i<3 ; i++)
169 return moment_un+same_orient*moment_deux ;
180 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
182 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
183 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
187 double* bornes =
new double [nzones+1] ;
189 for (
int i=nzones-1 ; i>0 ; i--) {
190 bornes[i] = courant ;
194 bornes[nzones] = __infinity ;
201 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
204 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
207 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
218 Tenseur shift_tot (shift_un+shift_deux) ;
220 shift_tot.
annule(0, nzones-2) ;
226 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
227 for (
int i=0 ; i<3 ; i++) {
228 k_total.set(i, i) = grad(i, i)-trace()/3. ;
229 for (
int j=i+1 ; j<3 ; j++)
230 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
233 for (
int lig=0 ; lig<3 ; lig++)
234 for (
int col=lig ; col<3 ; col++)
235 k_total.set(lig, col).mult_r_zec() ;
237 Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
239 for (
int i=0 ; i<3 ; i++)
240 vecteur_un.set(i) = k_total(0, i) ;
241 vecteur_un.change_triad (mapping.get_bvect_spher()) ;
242 Cmp integrant_un (vecteur_un(0)) ;
244 Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
246 for (
int i=0 ; i<3 ; i++)
247 vecteur_deux.set(i) = k_total(1, i) ;
248 vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
249 Cmp integrant_deux (vecteur_deux(0)) ;
258 double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
272 double pente = 2./(x_un-x_deux) ;
273 double constante = - (x_un+x_deux)/(x_un-x_deux) ;
278 double air_un, theta_un, phi_un ;
279 double air_deux, theta_deux, phi_deux ;
281 double* coloc =
new double[nr] ;
282 double* coef =
new double[nr] ;
283 int* deg =
new int[3] ;
284 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
286 for (
int i=0 ; i<nr ; i++) {
287 ksi = -
cos (M_PI*i/(nr-1)) ;
288 xabs = (ksi-constante)/pente ;
294 hole2.
psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
298 cfrcheb(deg, deg, coloc, deg, coef) ;
301 double* som =
new double[nr] ;
303 for (
int i=2 ; i<nr ; i+=2)
304 som[i] = 1./(i+1)-1./(i-1) ;
305 for (
int i=1 ; i<nr ; i+=2)
309 for (
int i=0 ; i<nr ; i++)
310 res += som[i]*coef[i] ;
322 Tbl Bhole_binaire::linear_momentum_systeme_inf()
const {
329 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
331 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
332 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
336 double* bornes =
new double [nzones+1] ;
338 for (
int i=nzones-1 ; i>0 ; i--) {
339 bornes[i] = courant ;
343 bornes[nzones] = __infinity ;
350 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
351 k_total.set_etat_qcq() ;
353 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
354 shift_un.set_etat_qcq() ;
356 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
357 shift_deux.set_etat_qcq() ;
367 Tenseur shift_tot (shift_un+shift_deux) ;
368 shift_tot.set_std_base() ;
369 shift_tot.annule(0, nzones-2) ;
372 Tenseur shift_old (shift_tot) ;
373 shift_tot.inc2_dzpuis() ;
374 shift_tot.dec2_dzpuis() ;
375 for (
int i=0 ; i< 3 ; i++)
379 Tenseur grad (shift_tot.gradient()) ;
380 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
381 for (
int i=0 ; i<3 ; i++) {
382 k_total.set(i, i) = grad(i, i)-trace()/3. ;
383 for (
int j=i+1 ; j<3 ; j++)
384 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
388 for (
int lig=0 ; lig<3 ; lig++)
389 for (
int col=lig ; col<3 ; col++)
390 k_total.set(lig, col).mult_r_zec() ;
392 for (
int comp=0 ; comp<3 ; comp++) {
393 Tenseur vecteur (mapping, 1, CON, mapping.get_bvect_cart()) ;
394 vecteur.set_etat_qcq() ;
395 for (
int i=0 ; i<3 ; i++)
396 vecteur.set(i) = k_total(i, comp) ;
397 vecteur.change_triad (mapping.get_bvect_spher()) ;
398 Cmp integrant (vecteur(0)) ;
400 res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
Coord xa
Absolute x coordinate.
void annule(int l)
Sets the Tenseur to zero in a given domain.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
void dec2_dzpuis()
dzpuis -= 2 ;
const Valeur & mult_sp() const
Returns applied to *this.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double get_ori_x() const
Returns the x coordinate of the origin.
double moment_systeme_inf()
Calculates the angular momentum of the black hole using the formula at infinity : where ...
Tenseur psi_auto
Part of generated by the hole.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
double distance_propre(const int nr=65) const
Calculation of the proper distance between the two spheres of inversion, along the x axis...
Cmp cos(const Cmp &)
Cosine.
Tenseur shift_auto
Part of generated by the hole.
void inc2_dzpuis()
dzpuis += 2 ;
double omega
Position of the axis of rotation.
Map_af & mp
Affine mapping.
Bhole hole1
Black hole one.
int get_nzone() const
Returns the number of domains.
const Valeur & mult_st() const
Returns applied to *this.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
double rayon
Radius of the horizon in LORENE's units.
Cmp pow(const Cmp &, int)
Power .
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Tenseur_sym tkij_tot
Total .
Coord ya
Absolute y coordinate.
const Valeur & mult_cp() const
Returns applied to *this.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
double moment_systeme_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and den...
double adm_systeme() const
Calculates the ADM mass of the system using : .
double komar_systeme() const
Calculates the Komar mass of the system using : .
Tenseur n_auto
Part of N generated by the hole.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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...
Valeur va
The numerical value of the Cmp.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Bhole hole2
Black hole two.