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