LORENE
bhole_glob.C
1/*
2 * Copyright (c) 2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
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 $" ;
24
25/*
26 * $Id: bhole_glob.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: bhole_glob.C,v $
28 * Revision 1.5 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.4 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.3 2005/08/29 15:10:14 p_grandclement
35 * Addition of things needed :
36 * 1) For BBH with different masses
37 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
38 * WORKING YET !!!
39 *
40 * Revision 1.2 2002/10/16 14:36:33 j_novak
41 * Reorganization of #include instructions of standard C++, in order to
42 * use experimental version 3 of gcc.
43 *
44 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
45 * LORENE
46 *
47 * Revision 2.5 2001/06/28 12:11:19 eric
48 * Correction d'un memory leak dans Bhole_binaire::moment_systeme_inf()
49 * (ajout de delete [] bornes).
50 *
51 * Revision 2.4 2001/05/07 12:24:19 phil
52 * correction de calcul de J
53 *
54 * Revision 2.3 2001/05/07 09:12:09 phil
55 * *** empty log message ***
56 *
57 * Revision 2.2 2001/03/22 10:49:47 phil
58 * *** empty log message ***
59 *
60 * Revision 2.1 2001/03/22 10:44:20 phil
61 * *** empty log message ***
62 *
63 * Revision 2.0 2001/03/01 08:18:13 phil
64 * *** empty log message ***
65 *
66 *
67 * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_glob.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
68 *
69 */
70
71
72
73//standard
74#include <cstdlib>
75#include <cmath>
76
77// Lorene
78#include "nbr_spx.h"
79#include "tenseur.h"
80#include "bhole.h"
81#include "proto.h"
82#include "utilitaires.h"
83#include "graphique.h"
84
85namespace Lorene {
87 Cmp der_un (hole1.psi_auto().dsdr()) ;
88 Cmp der_deux (hole2.psi_auto().dsdr()) ;
89
90 double masse = hole1.mp.integrale_surface_infini(der_un) +
92
93 masse /= -2*M_PI ;
94 return masse ;
95}
96
98 Cmp der_un (hole1.n_auto().dsdr()) ;
99 Cmp der_deux (hole2.n_auto().dsdr()) ;
100
101 double masse = hole1.mp.integrale_surface_infini(der_un) +
103
104 masse /= 4*M_PI ;
105 return masse ;
106}
107
109
110 if (omega == 0)
111 return 0 ;
112 else {
113 // Alignes ou non ?
114 double orientation_un = hole1.mp.get_rot_phi() ;
115 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
116 double orientation_deux = hole2.mp.get_rot_phi() ;
117 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
118 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
119
120 // Integrale sur le premiere horizon :
121 Cmp xa_un (hole1.mp) ;
122 xa_un = hole1.mp.xa ;
123 xa_un.std_base_scal() ;
124
125 Cmp ya_un (hole1.mp) ;
126 ya_un = hole1.mp.ya ;
127 ya_un.std_base_scal() ;
128
129 Tenseur vecteur_un (hole1.mp, 1, CON, hole1.mp.get_bvect_cart()) ;
130 vecteur_un.set_etat_qcq() ;
131 for (int i=0 ; i<3 ; i++)
132 vecteur_un.set(i) = (-ya_un*hole1.tkij_tot(0, i)+
133 xa_un * hole1.tkij_tot(1, i)) ;
134 vecteur_un.set_std_base() ;
135 vecteur_un.annule(hole1.mp.get_mg()->get_nzone()-1) ;
136 vecteur_un.change_triad (hole1.mp.get_bvect_spher()) ;
137
138 Cmp integrant_un (pow(hole1.psi_tot(), 6)*vecteur_un(0)) ;
139 integrant_un.std_base_scal() ;
140 double moment_un = hole1.mp.integrale_surface
141 (integrant_un, hole1.rayon)/8/M_PI ;
142
143 //Integrale sur le second horizon :
144 Cmp xa_deux (hole2.mp) ;
145 xa_deux = hole2.mp.xa ;
146 xa_deux.std_base_scal() ;
147
148 Cmp ya_deux (hole2.mp) ;
149 ya_deux = hole2.mp.ya ;
150 ya_deux.std_base_scal() ;
151
152 Tenseur vecteur_deux (hole2.mp, 1, CON, hole2.mp.get_bvect_cart()) ;
153 vecteur_deux.set_etat_qcq() ;
154 for (int i=0 ; i<3 ; i++)
155 vecteur_deux.set(i) = -ya_deux*hole2.tkij_tot(0, i)+
156 xa_deux * hole2.tkij_tot(1, i) ;
157 vecteur_deux.set_std_base() ;
158 vecteur_deux.annule(hole2.mp.get_mg()->get_nzone()-1) ;
159 vecteur_deux.change_triad (hole2.mp.get_bvect_spher()) ;
160
161 Cmp integrant_deux (pow(hole2.psi_tot(), 6)*vecteur_deux(0)) ;
162 integrant_deux.std_base_scal() ;
163 double moment_deux = hole2.mp.integrale_surface
164 (integrant_deux, hole2.rayon)/8/M_PI ;
165
166 return moment_un+same_orient*moment_deux ;
167 }
168}
169
171
172 if (omega == 0)
173 return 0 ;
174 else {
175 // Alignes ou non ?
176 double orientation_un = hole1.mp.get_rot_phi() ;
177 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
178 double orientation_deux = hole2.mp.get_rot_phi() ;
179 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
180 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
181
182 // On construit une grille et un mapping auxiliaire :
183 int nzones = hole1.mp.get_mg()->get_nzone() ;
184 double* bornes = new double [nzones+1] ;
185 double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
186 for (int i=nzones-1 ; i>0 ; i--) {
187 bornes[i] = courant ;
188 courant /= 2. ;
189 }
190 bornes[0] = 0 ;
191 bornes[nzones] = __infinity ;
192
193 Map_af mapping (*hole1.mp.get_mg(), bornes) ;
194
195 delete [] bornes ;
196
197 // On construit k_total dessus :
198 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
199 k_total.set_etat_qcq() ;
200
201 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
202 shift_un.set_etat_qcq() ;
203
204 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
205 shift_deux.set_etat_qcq() ;
206
207 shift_un.set(0).import_asymy(hole1.shift_auto(0)) ;
208 shift_un.set(1).import_symy(hole1.shift_auto(1)) ;
209 shift_un.set(2).import_asymy(hole1.shift_auto(2)) ;
210
211 shift_deux.set(0).import_asymy(same_orient*hole2.shift_auto(0)) ;
212 shift_deux.set(1).import_symy(same_orient*hole2.shift_auto(1)) ;
213 shift_deux.set(2).import_asymy(hole2.shift_auto(2)) ;
214
215 Tenseur shift_tot (shift_un+shift_deux) ;
216 shift_tot.set_std_base() ;
217 shift_tot.annule(0, nzones-2) ;
218 // On enleve les residus
219 shift_tot.inc2_dzpuis() ;
220 shift_tot.dec2_dzpuis() ;
221
222 Tenseur grad (shift_tot.gradient()) ;
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. ;
228 }
229
230 for (int lig=0 ; lig<3 ; lig++)
231 for (int col=lig ; col<3 ; col++)
232 k_total.set(lig, col).mult_r_zec() ;
233
234 Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
235 vecteur_un.set_etat_qcq() ;
236 for (int i=0 ; i<3 ; i++)
237 vecteur_un.set(i) = k_total(0, i) ;
238 vecteur_un.change_triad (mapping.get_bvect_spher()) ;
239 Cmp integrant_un (vecteur_un(0)) ;
240
241 Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
242 vecteur_deux.set_etat_qcq() ;
243 for (int i=0 ; i<3 ; i++)
244 vecteur_deux.set(i) = k_total(1, i) ;
245 vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
246 Cmp integrant_deux (vecteur_deux(0)) ;
247
248 // On multiplie par y et x :
249 integrant_un.va = integrant_un.va.mult_st() ;
250 integrant_un.va = integrant_un.va.mult_sp() ;
251
252 integrant_deux.va = integrant_deux.va.mult_st() ;
253 integrant_deux.va = integrant_deux.va.mult_cp() ;
254
255 double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
256
257 moment /= 8*M_PI ;
258 return moment ;
259 }
260}
261
262double Bhole_binaire::distance_propre(const int nr) const {
263
264 // On determine les rayons coordonnes des points limites de l'integrale :
265 double x_un = hole1.mp.get_ori_x() - hole1.rayon ;
266 double x_deux = hole2.mp.get_ori_x() + hole2.rayon ;
267
268 // Les coefficients du changement de variable :
269 double pente = 2./(x_un-x_deux) ;
270 double constante = - (x_un+x_deux)/(x_un-x_deux) ;
271
272
273 double ksi ; // variable d'integration.
274 double xabs ; // x reel.
275 double air_un, theta_un, phi_un ; // coordonnee spheriques 1
276 double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
277
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 ;
282
283 for (int i=0 ; i<nr ; i++) {
284 ksi = -cos (M_PI*i/(nr-1)) ;
285 xabs = (ksi-constante)/pente ;
286
287 hole1.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
288 hole2.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
289
290 coloc[i] = pow (hole1.psi_auto().val_point (air_un, theta_un, phi_un) +
291 hole2.psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
292 }
293
294 // On prend les coefficients de la fonction
295 cfrcheb(deg, deg, coloc, deg, coef) ;
296
297 // On integre
298 double* som = new double[nr] ;
299 som[0] = 2 ;
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)
303 som[i] = 0 ;
304
305 double res = 0 ;
306 for (int i=0 ; i<nr ; i++)
307 res += som[i]*coef[i] ;
308
309 res /= pente ;
310
311 delete [] deg ;
312 delete [] coef ;
313 delete [] coloc ;
314 delete [] som ;
315
316 return res ;
317}
318
319Tbl Bhole_binaire::linear_momentum_systeme_inf() const {
320
321 Tbl res (3) ;
322 res.set_etat_qcq() ;
323
324 // Alignes ou non ?
325 double orientation_un = hole1.mp.get_rot_phi() ;
326 assert ((orientation_un==0) || (orientation_un==M_PI)) ;
327 double orientation_deux = hole2.mp.get_rot_phi() ;
328 assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
329 int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
330
331 // On construit une grille et un mapping auxiliaire :
332 int nzones = hole1.mp.get_mg()->get_nzone() ;
333 double* bornes = new double [nzones+1] ;
334 double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
335 for (int i=nzones-1 ; i>0 ; i--) {
336 bornes[i] = courant ;
337 courant /= 2. ;
338 }
339 bornes[0] = 0 ;
340 bornes[nzones] = __infinity ;
341
342 Map_af mapping (*hole1.mp.get_mg(), bornes) ;
343
344 delete [] bornes ;
345
346 // On construit k_total dessus :
347 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
348 k_total.set_etat_qcq() ;
349
350 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
351 shift_un.set_etat_qcq() ;
352
353 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
354 shift_deux.set_etat_qcq() ;
355
356 shift_un.set(0).import_asymy(hole1.shift_auto(0)) ;
357 shift_un.set(1).import_symy(hole1.shift_auto(1)) ;
358 shift_un.set(2).import_asymy(hole1.shift_auto(2)) ;
359
360 shift_deux.set(0).import_asymy(same_orient*hole2.shift_auto(0)) ;
361 shift_deux.set(1).import_symy(same_orient*hole2.shift_auto(1)) ;
362 shift_deux.set(2).import_asymy(hole2.shift_auto(2)) ;
363
364 Tenseur shift_tot (shift_un+shift_deux) ;
365 shift_tot.set_std_base() ;
366 shift_tot.annule(0, nzones-2) ;
367
368 // On enleve les residus
369 Tenseur shift_old (shift_tot) ;
370 shift_tot.inc2_dzpuis() ;
371 shift_tot.dec2_dzpuis() ;
372 for (int i=0 ; i< 3 ; i++)
373 cout << max(diffrelmax(shift_tot(i), shift_old(i))) << " " ;
374 cout << endl ;
375
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. ;
382 }
383
384
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() ;
388
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)) ;
396
397 res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
398 }
399 return res ;
400}
401}
double moment_systeme_inf()
Calculates the angular momentum of the black hole using the formula at infinity : where .
Definition bhole_glob.C:170
double omega
Position of the axis of rotation.
Definition bhole.h:769
Bhole hole1
Black hole one.
Definition bhole.h:762
Bhole hole2
Black hole two.
Definition bhole.h:763
double komar_systeme() const
Calculates the Komar mass of the system using : .
Definition bhole_glob.C:97
double adm_systeme() const
Calculates the ADM mass of the system using : .
Definition bhole_glob.C:86
double moment_systeme_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and den...
Definition bhole_glob.C:108
double distance_propre(const int nr=65) const
Calculation of the proper distance between the two spheres of inversion, along the x axis.
Definition bhole_glob.C:262
Tenseur psi_auto
Part of generated by the hole.
Definition bhole.h:290
Tenseur shift_auto
Part of generated by the hole.
Definition bhole.h:297
Tenseur psi_tot
Total .
Definition bhole.h:292
Tenseur n_auto
Part of N generated by the hole.
Definition bhole.h:286
double rayon
Radius of the horizon in LORENE's units.
Definition bhole.h:274
Map_af & mp
Affine mapping.
Definition bhole.h:273
Tenseur_sym tkij_tot
Total .
Definition bhole.h:308
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
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.
Affine radial mapping.
Definition map.h:2027
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,...
Definition map.h:791
Coord ya
Absolute y coordinate.
Definition map.h:731
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,...
Definition map.C:302
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
double get_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Basic array class.
Definition tbl.h:161
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition tenseur.h:1253
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:900
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
void inc2_dzpuis()
dzpuis += 2 ;
Definition tenseur.C:1143
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:668
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.
Definition cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene prototypes.
Definition app_hor.h:64