LORENE
binhor_global.C
1 /*
2  * Copyright (c) 2005 Francois Limousin
3  * Jose Luis Jaramillo
4  *
5  * This file is part of LORENE.
6  *
7  * LORENE is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2 of the License, or
10  * (at your option) any later version.
11  *
12  * LORENE is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with LORENE; if not, write to the Free Software
19  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20  *
21  */
22 
23 
24 
25 
26 /*
27  * $Id: binhor_global.C,v 1.11 2016/12/05 16:17:46 j_novak Exp $
28  * $Log: binhor_global.C,v $
29  * Revision 1.11 2016/12/05 16:17:46 j_novak
30  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31  *
32  * Revision 1.10 2014/10/13 08:52:42 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.9 2014/10/06 15:13:01 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.8 2007/04/13 15:28:55 f_limousin
39  * Lots of improvements, generalisation to an arbitrary state of
40  * rotation, implementation of the spatial metric given by Samaya.
41  *
42  * Revision 1.7 2006/05/24 16:56:37 f_limousin
43  * Many small modifs.
44  *
45  * Revision 1.6 2005/09/24 08:31:38 f_limousin
46  * Improve the computation of moment_adm() and moment_hor().
47  *
48  * Revision 1.5 2005/09/13 18:33:15 f_limousin
49  * New function vv_bound_cart_bin(double) for computing binaries with
50  * berlin condition for the shift vector.
51  * Suppress all the symy and asymy in the importations.
52  *
53  * Revision 1.4 2005/06/09 16:12:04 f_limousin
54  * Implementation of amg_mom_adm().
55  *
56  * Revision 1.3 2005/04/29 14:02:44 f_limousin
57  * Important changes : manage the dependances between quantities (for
58  * instance psi and psi4). New function write_global(ost).
59  *
60  * Revision 1.2 2005/03/04 17:09:57 jl_jaramillo
61  * Change to avoid warnings
62  *
63  * Revision 1.1 2005/03/03 13:48:56 f_limousin
64  * First version
65  *
66  *
67  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_global.C,v 1.11 2016/12/05 16:17:46 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 "tensor.h"
80 #include "isol_hor.h"
81 #include "proto.h"
82 #include "utilitaires.h"
83 //#include "graphique.h"
84 
85 namespace Lorene {
86 double Bin_hor::adm_mass() const {
87 
88  Vector dpsi_un (hole1.psi_auto.derive_con(hole1.ff)) ;
89  Vector dpsi_deux (hole2.psi_auto.derive_con(hole2.ff)) ;
90 
91  Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
92  Vector ww (0.125*(hdirac1 - (hole1.hh.trace(hole1.ff)).
93  derive_con(hole1.ff))) ;
94 
95  double inf = hole1.mp.val_r(hole1.mp.get_mg()->get_nzone()-1, 1., 0., 0.) ;
96 
97  double masse = dpsi_un.flux(inf, hole1.ff) +
98  dpsi_deux.flux(inf, hole2.ff) +
99  ww.flux(inf, hole1.ff) ;
100  masse /= -2*M_PI ;
101  return masse ;
102 }
103 
104 double Bin_hor::komar_mass() const {
105 
106  Vector dnn_un (hole1.n_auto.derive_con(hole1.ff)) ;
107  Vector dnn_deux (hole2.n_auto.derive_con(hole2.ff)) ;
108 
109  Vector ww (contract(hole1.hh, 1, hole1.nn.derive_cov(hole1.ff), 0)) ;
110 
111  double inf = hole1.mp.val_r(hole1.mp.get_mg()->get_nzone()-1, 1., 0., 0.) ;
112 
113  double mass = dnn_un.flux(inf, hole1.ff) +
114  dnn_deux.flux(inf, hole2.ff) +
115  ww.flux(inf, hole1.ff) ;
116 
117  mass /= 4*M_PI ;
118  return mass ;
119 }
120 
121 double Bin_hor::ang_mom_hor() const {
122 
123  if (omega == 0)
124  return 0 ;
125  else {
126  // Alignement
127  double orientation_un = hole1.mp.get_rot_phi() ;
128  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
129  double orientation_deux = hole2.mp.get_rot_phi() ;
130  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
131  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
132 
133  // Integral on the first horizon :
134  Scalar xa_un (hole1.mp) ;
135  xa_un = hole1.mp.xa ;
136  xa_un.std_spectral_base() ;
137 
138  Scalar ya_un (hole1.mp) ;
139  ya_un = hole1.mp.ya ;
140  ya_un.std_spectral_base() ;
141 
142  Sym_tensor tkij_un (hole1.aa) ;
143  tkij_un.change_triad(hole1.mp.get_bvect_cart()) ;
144 
145  Vector vecteur_un (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
146  for (int i=1 ; i<=3 ; i++)
147  vecteur_un.set(i) = -ya_un*tkij_un(1, i)+
148  xa_un * tkij_un(2, i) ;
149  vecteur_un.annule_domain(hole1.mp.get_mg()->get_nzone()-1) ;
150  vecteur_un.change_triad (hole1.mp.get_bvect_spher()) ;
151 
152  Scalar integrant_un (hole1.get_psi4()*hole1.psi*hole1.psi
153  *vecteur_un(1)) ;
154  double moment_un = hole1.mp.integrale_surface
155  (integrant_un, hole1.radius+1e-12)/8/M_PI ;
156 
157  //Integral on the second horizon :
158  Scalar xa_deux (hole2.mp) ;
159  xa_deux = hole2.mp.xa ;
160  xa_deux.std_spectral_base() ;
161 
162  Scalar ya_deux (hole2.mp) ;
163  ya_deux = hole2.mp.ya ;
164  ya_deux.std_spectral_base() ;
165 
166  Sym_tensor tkij_deux (hole2.aa) ;
167  tkij_deux.change_triad(hole2.mp.get_bvect_cart()) ;
168 
169  Vector vecteur_deux (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
170  for (int i=1 ; i<=3 ; i++)
171  vecteur_deux.set(i) = -ya_deux*tkij_deux(1, i)+
172  xa_deux * tkij_deux(2, i) ;
173  vecteur_deux.annule_domain(hole2.mp.get_mg()->get_nzone()-1) ;
174  vecteur_deux.change_triad (hole2.mp.get_bvect_spher()) ;
175 
176  Scalar integrant_deux (hole2.get_psi4()*hole2.psi*hole2.psi
177  *vecteur_deux(1)) ;
178  double moment_deux = hole2.mp.integrale_surface
179  (integrant_deux, hole2.radius+1e-12)/8/M_PI ;
180 
181  return moment_un+same_orient*moment_deux ;
182  }
183 }
184 
185 
186 
187 double Bin_hor::ang_mom_adm() const {
188 /*
189  Scalar integrand_un (hole1.k_dd()(1,3) - hole1.gam_dd()(1,3) * hole1.trK) ;
190 
191  integrand_un.mult_rsint() ; // in order to pass from the triad components
192 
193  double mom = hole1.mp.integrale_surface_infini(integrand_un) ;
194  mom /= 8*M_PI ;
195  return mom ;
196 */
197 
198  if (omega == 0)
199  return 0 ;
200  else {
201  // Alignement
202  double orientation_un = hole1.mp.get_rot_phi() ;
203  assert ((orientation_un==0) || (orientation_un==M_PI)) ;
204  double orientation_deux = hole2.mp.get_rot_phi() ;
205  assert ((orientation_deux==0) || (orientation_deux==M_PI)) ;
206  int same_orient = (orientation_un == orientation_deux) ? 1 : -1 ;
207 
208  // Construction of an auxiliar grid and mapping
209  int nzones = hole1.mp.get_mg()->get_nzone() ;
210  double* bornes = new double [nzones+1] ;
211  double courant = (hole1.mp.get_ori_x()-hole2.mp.get_ori_x())+1 ;
212  for (int i=nzones-1 ; i>0 ; i--) {
213  bornes[i] = courant ;
214  courant /= 2. ;
215  }
216  bornes[0] = 0 ;
217  bornes[nzones] = __infinity ;
218 
219  Map_af mapping (*hole1.mp.get_mg(), bornes) ;
220 
221  delete [] bornes ;
222 
223  // Construction of k_total
224  Sym_tensor k_total (mapping, CON, mapping.get_bvect_cart()) ;
225 
226  Vector shift_un (mapping, CON, mapping.get_bvect_cart()) ;
227  Vector shift_deux (mapping, CON, mapping.get_bvect_cart()) ;
228 
229  Vector beta_un (hole1.beta_auto) ;
230  Vector beta_deux (hole2.beta_auto) ;
231  beta_un.change_triad(hole1.mp.get_bvect_cart()) ;
232  beta_deux.change_triad(hole2.mp.get_bvect_cart()) ;
233 
234  shift_un.set(1).import(beta_un(1)) ;
235  shift_un.set(2).import(beta_un(2)) ;
236  shift_un.set(3).import(beta_un(3)) ;
237 
238  shift_deux.set(1).import(same_orient*beta_deux(1)) ;
239  shift_deux.set(2).import(same_orient*beta_deux(2)) ;
240  shift_deux.set(3).import(beta_deux(3)) ;
241 
242  Vector shift_tot (shift_un+shift_deux) ;
243  shift_tot.std_spectral_base() ;
244  shift_tot.annule(0, nzones-2) ;
245 
246  // Substract the residuals
247  shift_tot.inc_dzpuis(2) ;
248  shift_tot.dec_dzpuis(2) ;
249 
250  const Metric_flat& flat0 (mapping.flat_met_cart()) ;
251 
252  k_total = shift_tot.ope_killing_conf(flat0) / 2. ;
253  //- flat0.con() * hole1.trK;
254 
255  for (int lig=1 ; lig<=3 ; lig++)
256  for (int col=lig ; col<=3 ; col++)
257  k_total.set(lig, col).mult_r_ced() ;
258 
259  Vector vecteur_un (mapping, CON, mapping.get_bvect_cart()) ;
260  for (int i=1 ; i<=3 ; i++)
261  vecteur_un.set(i) = k_total(1, i) ;
262  vecteur_un.change_triad (mapping.get_bvect_spher()) ;
263  Scalar integrant_un (vecteur_un(1)) ;
264 
265  Vector vecteur_deux (mapping, CON, mapping.get_bvect_cart()) ;
266  for (int i=1 ; i<=3 ; i++)
267  vecteur_deux.set(i) = k_total(2, i) ;
268  vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
269  Scalar integrant_deux (vecteur_deux(1)) ;
270 
271  // Multiplication by y and x :
272  integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
273  .mult_st() ;
274  integrant_un.set_spectral_va() = integrant_un.get_spectral_va()
275  .mult_sp() ;
276 
277  integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
278  .mult_st() ;
279  integrant_deux.set_spectral_va() = integrant_deux.get_spectral_va()
280  .mult_cp() ;
281 
282  double moment = mapping.integrale_surface_infini (-integrant_un
283  +integrant_deux) ;
284 
285  moment /= 8*M_PI ;
286 
287  return moment ;
288  }
289 }
290 
291 
292 /*
293 double Bin_hor::proper_distance(const int nr) const {
294 
295 
296  // On determine les rayons coordonnes des points limites de l'integrale :
297  double x_un = hole1.mp.get_ori_x() - hole1.rayon ;
298  double x_deux = hole2.mp.get_ori_x() + hole2.rayon ;
299 
300  // Les coefficients du changement de variable :
301  double pente = 2./(x_un-x_deux) ;
302  double constante = - (x_un+x_deux)/(x_un-x_deux) ;
303 
304 
305  double ksi ; // variable d'integration.
306  double xabs ; // x reel.
307  double air_un, theta_un, phi_un ; // coordonnee spheriques 1
308  double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
309 
310  double* coloc = new double[nr] ;
311  double* coef = new double[nr] ;
312  int* deg = new int[3] ;
313  deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
314 
315  for (int i=0 ; i<nr ; i++) {
316  ksi = -cos (M_PI*i/(nr-1)) ;
317  xabs = (ksi-constante)/pente ;
318 
319  hole1.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
320  hole2.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
321 
322  coloc[i] = pow (hole1.psi_auto().val_point (air_un, theta_un, phi_un) +
323  hole2.psi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
324  }
325 
326  // On prend les coefficients de la fonction
327  cfrcheb(deg, deg, coloc, deg, coef) ;
328 
329  // On integre
330  double* som = new double[nr] ;
331  som[0] = 2 ;
332  for (int i=2 ; i<nr ; i+=2)
333  som[i] = 1./(i+1)-1./(i-1) ;
334  for (int i=1 ; i<nr ; i+=2)
335  som[i] = 0 ;
336 
337  double res = 0 ;
338  for (int i=0 ; i<nr ; i++)
339  res += som[i]*coef[i] ;
340 
341  res /= pente ;
342 
343  delete [] deg ;
344  delete [] coef ;
345  delete [] coloc ;
346  delete [] som ;
347 
348  return res ;
349 
350 }
351 */
352 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
double flux(double radius, const Metric &met) const
Computes the flux of the vector accross a sphere r = const.
Definition: vector.C:813
double omega
Angular velocity.
Definition: isol_hor.h:1348
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Valeur & mult_sp() const
Returns applied to *this.
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
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Definition: map_af_radius.C:99
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
const Scalar & get_psi4() const
Conformal factor .
Definition: single_hor.C:291
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
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:817
Sym_tensor aa
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:971
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
Scalar psi
Conformal factor .
Definition: isol_hor.h:930
Sym_tensor hh
Deviation metric.
Definition: isol_hor.h:983
void mult_r_ced()
Multiplication by r in the compactified external domain (CED), the dzpuis flag is not changed...
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Valeur & mult_st() const
Returns applied to *this.
Tensor handling.
Definition: tensor.h:294
Scalar n_auto
Lapse function .
Definition: isol_hor.h:915
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
double komar_mass() const
Calculates the Komar mass of the system using : .
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
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:825
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
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar psi_auto
Conformal factor .
Definition: isol_hor.h:924
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
double ang_mom_adm() const
Calculates the angular momentum of the black hole.
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:906
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition: tensor.C:680
double adm_mass() const
Calculates the ADM mass of the system.
Definition: binhor_global.C:86
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Scalar nn
Lapse function .
Definition: isol_hor.h:921
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.