LORENE
hole_bhns_global.C
1/*
2 * Methods of class Hole_bhns to compute global quantities
3 *
4 * (see file hole_bhns.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005,2007 Keisuke Taniguchi
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char hole_bhns_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_global.C,v 1.5 2014/10/13 08:53:00 j_novak Exp $" ;
29
30/*
31 * $Id: hole_bhns_global.C,v 1.5 2014/10/13 08:53:00 j_novak Exp $
32 * $Log: hole_bhns_global.C,v $
33 * Revision 1.5 2014/10/13 08:53:00 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2014/10/06 15:13:10 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.3 2008/07/02 21:10:15 k_taniguchi
40 * A bug removed.
41 *
42 * Revision 1.2 2008/05/15 19:07:26 k_taniguchi
43 * Introduction of the quasilocal spin angular momentum.
44 *
45 * Revision 1.1 2007/06/22 01:25:15 k_taniguchi
46 * *** empty log message ***
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_global.C,v 1.5 2014/10/13 08:53:00 j_novak Exp $
50 *
51 */
52
53// C++ headers
54//#include <>
55
56// C headers
57#include <cmath>
58
59// Lorene headers
60#include "hole_bhns.h"
61#include "unites.h"
62#include "utilitaires.h"
63
64 //-----------------------------------------//
65 // Irreducible mass of BH //
66 //-----------------------------------------//
67
68namespace Lorene {
70
71 // Fundamental constants and units
72 // -------------------------------
73 using namespace Unites ;
74
75 if (p_mass_irr_bhns == 0x0) { // a new computation is required
76
77 Scalar psi4(mp) ;
78 psi4 = pow(confo_tot, 4.) ;
79 psi4.std_spectral_base() ;
80 psi4.annule_domain(0) ;
81 psi4.raccord(1) ;
82
83 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
84
85 Map_af& mp_aff= dynamic_cast<Map_af&>(mp) ;
86
87 double a_ah = mp_aff.integrale_surface(psi4, radius_ah) ;
88 double mirr = sqrt(a_ah/16./M_PI) / ggrav ;
89
90 p_mass_irr_bhns = new double( mirr ) ;
91
92 }
93
94 return *p_mass_irr_bhns ;
95
96}
97
98 //----------------------------------------------------------//
99 // Quasilocal spin angular momentum of BH //
100 //----------------------------------------------------------//
101
102double Hole_bhns::spin_am_bhns(const Tbl& xi_i, const double& phi_i,
103 const double& theta_i, const int& nrk_phi,
104 const int& nrk_theta) const {
105
106 // Fundamental constants and units
107 // -------------------------------
108 using namespace Unites ;
109
110 if (p_spin_am_bhns == 0x0) { // a new computation is required
111
112 double mass = ggrav * mass_bh ;
113
114 Scalar rr(mp) ;
115 rr = mp.r ;
116 rr.std_spectral_base() ;
117
118 Scalar st(mp) ;
119 st = mp.sint ;
120 st.std_spectral_base() ;
121 Scalar ct(mp) ;
122 ct = mp.cost ;
123 ct.std_spectral_base() ;
124 Scalar sp(mp) ;
125 sp = mp.sinp ;
126 sp.std_spectral_base() ;
127 Scalar cp(mp) ;
128 cp = mp.cosp ;
129 cp.std_spectral_base() ;
130
131 Vector ll(mp, CON, mp.get_bvect_cart()) ;
132 ll.set_etat_qcq() ;
133 ll.set(1) = st % cp ;
134 ll.set(2) = st % sp ;
135 ll.set(3) = ct ;
136 ll.std_spectral_base() ;
137
138 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
139
140 if (kerrschild) {
141
142 cout << "Not yet prepared!!!" << endl ;
143 abort() ;
144
145 }
146 else { // Isotropic coordinates
147
148 // Sets C/M^2 for each case of the lapse boundary condition
149 // --------------------------------------------------------
150 double cc ;
151
152 if (bc_lapconf_nd) { // Neumann boundary condition
153 if (bc_lapconf_fs) { // First condition
154 // d(\alpha \psi)/dr = 0
155 // ---------------------
156 cc = 2. * (sqrt(13.) - 1.) / 3. ;
157 }
158 else { // Second condition
159 // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
160 // -----------------------------------------
161 cc = 4. / 3. ;
162 }
163 }
164 else { // Dirichlet boundary condition
165 if (bc_lapconf_fs) { // First condition
166 // (\alpha \psi) = 1/2
167 // -------------------
168 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
169 abort() ;
170 }
171 else { // Second condition
172 // (\alpha \psi) = 1/sqrt(2.) \psi_KS
173 // ----------------------------------
174 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
175 abort() ;
176 // cc = 2. * sqrt(2.) ;
177 }
178 }
179
180 Scalar r_are(mp) ;
182 r_are.std_spectral_base() ;
183
184 // Killing vector of the spherical components
185 Vector killing_spher(mp, COV, mp.get_bvect_spher()) ;
186 killing_spher.set_etat_qcq() ;
187 killing_spher = killing_vect(xi_i, phi_i, theta_i,
188 nrk_phi, nrk_theta) ;
189 killing_spher.std_spectral_base() ;
190
191 killing_spher.set(2) = confo_tot * confo_tot * radius_ah
192 * killing_spher(2) ;
193 killing_spher.set(3) = confo_tot * confo_tot * radius_ah
194 * killing_spher(3) ;
195 // killing_spher(3) is already divided by sin(theta)
196 killing_spher.std_spectral_base() ;
197
198 // Killing vector of the Cartesian components
199 Vector killing(mp, COV, mp.get_bvect_cart()) ;
200 killing.set_etat_qcq() ;
201 killing.set(1) = (killing_spher(2) * ct * cp - killing_spher(3) * sp)
202 / radius_ah ;
203 killing.set(2) = (killing_spher(2) * ct * sp + killing_spher(3) * cp)
204 / radius_ah ;
205 killing.set(3) = - killing_spher(2) * st / radius_ah ;
206 killing.std_spectral_base() ;
207
208 // Surface integral <- dzpuis should be 0
209 // --------------------------------------
210 // Source terms in the surface integral
211 Scalar source_1(mp) ;
212 source_1 = (ll(1) * (taij_tot_rs(1,1) * killing(1)
213 + taij_tot_rs(1,2) * killing(2)
214 + taij_tot_rs(1,3) * killing(3))
215 + ll(2) * (taij_tot_rs(2,1) * killing(1)
216 + taij_tot_rs(2,2) * killing(2)
217 + taij_tot_rs(2,3) * killing(3))
218 + ll(3) * (taij_tot_rs(3,1) * killing(1)
219 + taij_tot_rs(3,2) * killing(2)
220 + taij_tot_rs(3,3) * killing(3)))
221 / pow(confo_tot, 4.) ;
222 source_1.std_spectral_base() ;
223 source_1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
224
225 Scalar source_2(mp) ;
226 source_2 = -2. * pow(confo_tot, 3.) * mass * mass * cc
227 * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
228 * (ll(1)*killing(1) + ll(2)*killing(2) + ll(3)*killing(3))
229 / lapconf_tot / pow(r_are*rr, 3.) ;
230 source_2.std_spectral_base() ;
231
232 Scalar source_surf(mp) ;
233 source_surf = source_1 + source_2 ;
234 source_surf.std_spectral_base() ;
235 source_surf.annule_domain(0) ;
236 source_surf.raccord(1) ;
237
238 Map_af& mp_aff= dynamic_cast<Map_af&>(mp) ;
239
240 double spin = mp_aff.integrale_surface(source_surf, radius_ah) ;
241 double spin_angmom = 0.5 * spin / qpig ;
242
243 p_spin_am_bhns = new double( spin_angmom ) ;
244
245 }
246
247 }
248
249 return *p_spin_am_bhns ;
250
251}
252}
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
double spin_am_bhns(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Spin angular momentum of the black hole.
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Definition hole_bhns.h:190
virtual double mass_irr_bhns() const
Irreducible mass of the black hole.
Vector killing_vect(const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Compute the Killing vector of a black hole normalized so that its affine length is 2 M_PI.
bool bc_lapconf_fs
true for the first type BC for the lapconf function, false for the second type BH
Definition hole_bhns.h:78
bool bc_lapconf_nd
true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH
Definition hole_bhns.h:73
Scalar lapconf_tot
Total lapconf function.
Definition hole_bhns.h:101
Scalar confo_tot
Total conformal factor.
Definition hole_bhns.h:169
double * p_spin_am_bhns
Irreducible mass of BH.
Definition hole_bhns.h:248
Affine radial mapping.
Definition map.h:2027
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Coord cosp
Definition map.h:724
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 sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord sinp
Definition map.h:723
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord cost
Definition map.h:722
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
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
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Basic array class.
Definition tbl.h:161
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
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.