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 
23 
24 
25 /*
26  * $Id: bhole_glob.C,v 1.6 2016/12/05 16:17:45 j_novak Exp $
27  * $Log: bhole_glob.C,v $
28  * Revision 1.6 2016/12/05 16:17:45 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.5 2014/10/13 08:52:40 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.4 2014/10/06 15:12:58 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.3 2005/08/29 15:10:14 p_grandclement
38  * Addition of things needed :
39  * 1) For BBH with different masses
40  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
41  * WORKING YET !!!
42  *
43  * Revision 1.2 2002/10/16 14:36:33 j_novak
44  * Reorganization of #include instructions of standard C++, in order to
45  * use experimental version 3 of gcc.
46  *
47  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
48  * LORENE
49  *
50  * Revision 2.5 2001/06/28 12:11:19 eric
51  * Correction d'un memory leak dans Bhole_binaire::moment_systeme_inf()
52  * (ajout de delete [] bornes).
53  *
54  * Revision 2.4 2001/05/07 12:24:19 phil
55  * correction de calcul de J
56  *
57  * Revision 2.3 2001/05/07 09:12:09 phil
58  * *** empty log message ***
59  *
60  * Revision 2.2 2001/03/22 10:49:47 phil
61  * *** empty log message ***
62  *
63  * Revision 2.1 2001/03/22 10:44:20 phil
64  * *** empty log message ***
65  *
66  * Revision 2.0 2001/03/01 08:18:13 phil
67  * *** empty log message ***
68  *
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_glob.C,v 1.6 2016/12/05 16:17:45 j_novak Exp $
71  *
72  */
73 
74 
75 
76 //standard
77 #include <cstdlib>
78 #include <cmath>
79 
80 // Lorene
81 #include "nbr_spx.h"
82 #include "tenseur.h"
83 #include "bhole.h"
84 #include "proto.h"
85 #include "utilitaires.h"
86 #include "graphique.h"
87 
88 namespace Lorene {
90  Cmp der_un (hole1.psi_auto().dsdr()) ;
91  Cmp der_deux (hole2.psi_auto().dsdr()) ;
92 
93  double masse = hole1.mp.integrale_surface_infini(der_un) +
95 
96  masse /= -2*M_PI ;
97  return masse ;
98 }
99 
101  Cmp der_un (hole1.n_auto().dsdr()) ;
102  Cmp der_deux (hole2.n_auto().dsdr()) ;
103 
104  double masse = hole1.mp.integrale_surface_infini(der_un) +
105  hole2.mp.integrale_surface_infini(der_deux) ;
106 
107  masse /= 4*M_PI ;
108  return masse ;
109 }
110 
112 
113  if (omega == 0)
114  return 0 ;
115  else {
116  // Alignes ou non ?
117  double orientation_un = hole1.mp.get_rot_phi() ;
118  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
119  double orientation_deux = hole2.mp.get_rot_phi() ;
120  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
121  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
122 
123  // Integrale sur le premiere horizon :
124  Cmp xa_un (hole1.mp) ;
125  xa_un = hole1.mp.xa ;
126  xa_un.std_base_scal() ;
127 
128  Cmp ya_un (hole1.mp) ;
129  ya_un = hole1.mp.ya ;
130  ya_un.std_base_scal() ;
131 
132  Tenseur vecteur_un (hole1.mp, 1, CON, hole1.mp.get_bvect_cart()) ;
133  vecteur_un.set_etat_qcq() ;
134  for (int i=0 ; i<3 ; i++)
135  vecteur_un.set(i) = (-ya_un*hole1.tkij_tot(0, i)+
136  xa_un * hole1.tkij_tot(1, i)) ;
137  vecteur_un.set_std_base() ;
138  vecteur_un.annule(hole1.mp.get_mg()->get_nzone()-1) ;
139  vecteur_un.change_triad (hole1.mp.get_bvect_spher()) ;
140 
141  Cmp integrant_un (pow(hole1.psi_tot(), 6)*vecteur_un(0)) ;
142  integrant_un.std_base_scal() ;
143  double moment_un = hole1.mp.integrale_surface
144  (integrant_un, hole1.rayon)/8/M_PI ;
145 
146  //Integrale sur le second horizon :
147  Cmp xa_deux (hole2.mp) ;
148  xa_deux = hole2.mp.xa ;
149  xa_deux.std_base_scal() ;
150 
151  Cmp ya_deux (hole2.mp) ;
152  ya_deux = hole2.mp.ya ;
153  ya_deux.std_base_scal() ;
154 
155  Tenseur vecteur_deux (hole2.mp, 1, CON, hole2.mp.get_bvect_cart()) ;
156  vecteur_deux.set_etat_qcq() ;
157  for (int i=0 ; i<3 ; i++)
158  vecteur_deux.set(i) = -ya_deux*hole2.tkij_tot(0, i)+
159  xa_deux * hole2.tkij_tot(1, i) ;
160  vecteur_deux.set_std_base() ;
161  vecteur_deux.annule(hole2.mp.get_mg()->get_nzone()-1) ;
162  vecteur_deux.change_triad (hole2.mp.get_bvect_spher()) ;
163 
164  Cmp integrant_deux (pow(hole2.psi_tot(), 6)*vecteur_deux(0)) ;
165  integrant_deux.std_base_scal() ;
166  double moment_deux = hole2.mp.integrale_surface
167  (integrant_deux, hole2.rayon)/8/M_PI ;
168 
169  return moment_un+same_orient*moment_deux ;
170  }
171 }
172 
174 
175  if (omega == 0)
176  return 0 ;
177  else {
178  // Alignes ou non ?
179  double orientation_un = hole1.mp.get_rot_phi() ;
180  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
181  double orientation_deux = hole2.mp.get_rot_phi() ;
182  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
183  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
184 
185  // On construit une grille et un mapping auxiliaire :
186  int nzones = hole1.mp.get_mg()->get_nzone() ;
187  double* bornes = new double [nzones+1] ;
188  double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
189  for (int i=nzones-1 ; i>0 ; i--) {
190  bornes[i] = courant ;
191  courant /= 2. ;
192  }
193  bornes[0] = 0 ;
194  bornes[nzones] = __infinity ;
195 
196  Map_af mapping (*hole1.mp.get_mg(), bornes) ;
197 
198  delete [] bornes ;
199 
200  // On construit k_total dessus :
201  Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
202  k_total.set_etat_qcq() ;
203 
204  Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
205  shift_un.set_etat_qcq() ;
206 
207  Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
208  shift_deux.set_etat_qcq() ;
209 
210  shift_un.set(0).import_asymy(hole1.shift_auto(0)) ;
211  shift_un.set(1).import_symy(hole1.shift_auto(1)) ;
212  shift_un.set(2).import_asymy(hole1.shift_auto(2)) ;
213 
214  shift_deux.set(0).import_asymy(same_orient*hole2.shift_auto(0)) ;
215  shift_deux.set(1).import_symy(same_orient*hole2.shift_auto(1)) ;
216  shift_deux.set(2).import_asymy(hole2.shift_auto(2)) ;
217 
218  Tenseur shift_tot (shift_un+shift_deux) ;
219  shift_tot.set_std_base() ;
220  shift_tot.annule(0, nzones-2) ;
221  // On enleve les residus
222  shift_tot.inc2_dzpuis() ;
223  shift_tot.dec2_dzpuis() ;
224 
225  Tenseur grad (shift_tot.gradient()) ;
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. ;
231  }
232 
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() ;
236 
237  Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
238  vecteur_un.set_etat_qcq() ;
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)) ;
243 
244  Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
245  vecteur_deux.set_etat_qcq() ;
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)) ;
250 
251  // On multiplie par y et x :
252  integrant_un.va = integrant_un.va.mult_st() ;
253  integrant_un.va = integrant_un.va.mult_sp() ;
254 
255  integrant_deux.va = integrant_deux.va.mult_st() ;
256  integrant_deux.va = integrant_deux.va.mult_cp() ;
257 
258  double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
259 
260  moment /= 8*M_PI ;
261  return moment ;
262  }
263 }
264 
265 double Bhole_binaire::distance_propre(const int nr) const {
266 
267  // On determine les rayons coordonnes des points limites de l'integrale :
268  double x_un = hole1.mp.get_ori_x() - hole1.rayon ;
269  double x_deux = hole2.mp.get_ori_x() + hole2.rayon ;
270 
271  // Les coefficients du changement de variable :
272  double pente = 2./(x_un-x_deux) ;
273  double constante = - (x_un+x_deux)/(x_un-x_deux) ;
274 
275 
276  double ksi ; // variable d'integration.
277  double xabs ; // x reel.
278  double air_un, theta_un, phi_un ; // coordonnee spheriques 1
279  double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
280 
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 ;
285 
286  for (int i=0 ; i<nr ; i++) {
287  ksi = -cos (M_PI*i/(nr-1)) ;
288  xabs = (ksi-constante)/pente ;
289 
290  hole1.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
291  hole2.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
292 
293  coloc[i] = pow (hole1.psi_auto().val_point (air_un, theta_un, phi_un) +
294  hole2.psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
295  }
296 
297  // On prend les coefficients de la fonction
298  cfrcheb(deg, deg, coloc, deg, coef) ;
299 
300  // On integre
301  double* som = new double[nr] ;
302  som[0] = 2 ;
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)
306  som[i] = 0 ;
307 
308  double res = 0 ;
309  for (int i=0 ; i<nr ; i++)
310  res += som[i]*coef[i] ;
311 
312  res /= pente ;
313 
314  delete [] deg ;
315  delete [] coef ;
316  delete [] coloc ;
317  delete [] som ;
318 
319  return res ;
320 }
321 
322 Tbl Bhole_binaire::linear_momentum_systeme_inf() const {
323 
324  Tbl res (3) ;
325  res.set_etat_qcq() ;
326 
327  // Alignes ou non ?
328  double orientation_un = hole1.mp.get_rot_phi() ;
329  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
330  double orientation_deux = hole2.mp.get_rot_phi() ;
331  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
332  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
333 
334  // On construit une grille et un mapping auxiliaire :
335  int nzones = hole1.mp.get_mg()->get_nzone() ;
336  double* bornes = new double [nzones+1] ;
337  double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
338  for (int i=nzones-1 ; i>0 ; i--) {
339  bornes[i] = courant ;
340  courant /= 2. ;
341  }
342  bornes[0] = 0 ;
343  bornes[nzones] = __infinity ;
344 
345  Map_af mapping (*hole1.mp.get_mg(), bornes) ;
346 
347  delete [] bornes ;
348 
349  // On construit k_total dessus :
350  Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
351  k_total.set_etat_qcq() ;
352 
353  Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
354  shift_un.set_etat_qcq() ;
355 
356  Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
357  shift_deux.set_etat_qcq() ;
358 
359  shift_un.set(0).import_asymy(hole1.shift_auto(0)) ;
360  shift_un.set(1).import_symy(hole1.shift_auto(1)) ;
361  shift_un.set(2).import_asymy(hole1.shift_auto(2)) ;
362 
363  shift_deux.set(0).import_asymy(same_orient*hole2.shift_auto(0)) ;
364  shift_deux.set(1).import_symy(same_orient*hole2.shift_auto(1)) ;
365  shift_deux.set(2).import_asymy(hole2.shift_auto(2)) ;
366 
367  Tenseur shift_tot (shift_un+shift_deux) ;
368  shift_tot.set_std_base() ;
369  shift_tot.annule(0, nzones-2) ;
370 
371  // On enleve les residus
372  Tenseur shift_old (shift_tot) ;
373  shift_tot.inc2_dzpuis() ;
374  shift_tot.dec2_dzpuis() ;
375  for (int i=0 ; i< 3 ; i++)
376  cout << max(diffrelmax(shift_tot(i), shift_old(i))) << " " ;
377  cout << endl ;
378 
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. ;
385  }
386 
387 
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() ;
391 
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)) ;
399 
400  res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
401  }
402  return res ;
403 }
404 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:916
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
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 ;
Definition: tenseur.C:1146
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.
Definition: tenseur.C:1186
Class intended to describe tensors with a symmetry on the two last indices *** DEPRECATED : use class...
Definition: tenseur.h:1256
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
double moment_systeme_inf()
Calculates the angular momentum of the black hole using the formula at infinity : where ...
Definition: bhole_glob.C:173
Tenseur psi_auto
Part of generated by the hole.
Definition: bhole.h:290
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
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:265
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
Tenseur shift_auto
Part of generated by the hole.
Definition: bhole.h:297
void inc2_dzpuis()
dzpuis += 2 ;
Definition: tenseur.C:1159
double omega
Position of the axis of rotation.
Definition: bhole.h:769
Map_af & mp
Affine mapping.
Definition: bhole.h:273
Bhole hole1
Black hole one.
Definition: bhole.h:762
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Valeur & mult_st() const
Returns applied to *this.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
double rayon
Radius of the horizon in LORENE&#39;s units.
Definition: bhole.h:274
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
Tenseur_sym tkij_tot
Total .
Definition: bhole.h:308
Coord ya
Absolute y coordinate.
Definition: map.h:749
Affine radial mapping.
Definition: map.h:2048
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.
Definition: map.h:809
Tenseur psi_tot
Total .
Definition: bhole.h:292
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:111
double adm_systeme() const
Calculates the ADM mass of the system using : .
Definition: bhole_glob.C:89
double komar_systeme() const
Calculates the Komar mass of the system using : .
Definition: bhole_glob.C:100
Tenseur n_auto
Part of N generated by the hole.
Definition: bhole.h:286
Basic array class.
Definition: tbl.h:164
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
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:305
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558
Bhole hole2
Black hole two.
Definition: bhole.h:763