LORENE
blackhole_global.C
1/*
2 * Methods of class Black_hole to compute global quantities
3 *
4 * (see file blackhole.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 blackhole_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_global.C,v 1.5 2014/10/13 08:52:46 j_novak Exp $" ;
29
30/*
31 * $Id: blackhole_global.C,v 1.5 2014/10/13 08:52:46 j_novak Exp $
32 * $Log: blackhole_global.C,v $
33 * Revision 1.5 2014/10/13 08:52:46 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2014/10/06 15:13:02 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.3 2008/07/02 20:45:58 k_taniguchi
40 * Addition of routines to compute angular momentum.
41 *
42 * Revision 1.2 2008/05/15 19:28:03 k_taniguchi
43 * Change of some parameters.
44 *
45 * Revision 1.1 2007/06/22 01:19:51 k_taniguchi
46 * *** empty log message ***
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_global.C,v 1.5 2014/10/13 08:52:46 j_novak Exp $
50 *
51 */
52
53// C++ headers
54//#include <>
55
56// C headers
57#include <cmath>
58
59// Lorene headers
60#include "blackhole.h"
61#include "unites.h"
62#include "utilitaires.h"
63
64 //-----------------------------------------//
65 // Irreducible mass of BH //
66 //-----------------------------------------//
67
68namespace Lorene {
69double Black_hole::mass_irr() const {
70
71 // Fundamental constants and units
72 // -------------------------------
73 using namespace Unites ;
74
75 if (p_mass_irr == 0x0) { // a new computation is required
76
77 Scalar psi4(mp) ;
78 psi4 = pow(confo, 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(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 = new double( mirr ) ;
91
92 }
93
94 return *p_mass_irr ;
95
96}
97
98
99 //---------------------------//
100 // ADM mass //
101 //---------------------------//
102
103double Black_hole::mass_adm() const {
104
105 // Fundamental constants and units
106 // -------------------------------
107 using namespace Unites ;
108
109 if (p_mass_adm == 0x0) { // a new computation is required
110
111 double madm ;
112 double integ_s ;
113 double integ_v ;
114
115 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
116 Map_af mp_aff(mp) ;
117
118 Scalar source_surf(mp) ;
119 source_surf.set_etat_qcq() ;
120 Scalar source_volm(mp) ;
121 source_volm.set_etat_qcq() ;
122
123 Scalar rr(mp) ;
124 rr = mp.r ;
125 rr.std_spectral_base() ;
126 Scalar st(mp) ;
127 st = mp.sint ;
128 st.std_spectral_base() ;
129 Scalar ct(mp) ;
130 ct = mp.cost ;
131 ct.std_spectral_base() ;
132 Scalar sp(mp) ;
133 sp = mp.sinp ;
134 sp.std_spectral_base() ;
135 Scalar cp(mp) ;
136 cp = mp.cosp ;
137 cp.std_spectral_base() ;
138
139 Vector ll(mp, CON, mp.get_bvect_cart()) ;
140 ll.set_etat_qcq() ;
141 ll.set(1) = st * cp ;
142 ll.set(2) = st * sp ;
143 ll.set(3) = ct ;
144 ll.std_spectral_base() ;
145
146 Scalar lldconf = confo.dsdr() ;
147 lldconf.std_spectral_base() ;
148
149 if (kerrschild) {
150
151 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
152 abort() ;
153 /*
154 Scalar divshf(mp) ;
155 divshf = shift_rs(1).deriv(1) + shift_rs(2).deriv(2)
156 + shift_rs(3).deriv(3) ;
157 divshf.std_spectral_base() ;
158 divshf.dec_dzpuis(2) ;
159
160 Scalar llshift(mp) ;
161 llshift = ll(1)%shift_rs(1) + ll(2)%shift_rs(2)
162 + ll(3)%shift_rs(3) ;
163 llshift.std_spectral_base() ;
164
165 Scalar lldllsh = llshift.dsdr() ;
166 lldllsh.std_spectral_base() ;
167 lldllsh.dec_dzpuis(2) ;
168
169 double mass = ggrav * mass_bh ;
170
171 // Surface integral
172 // ----------------
173 source_surf = confo/rr - 2.*pow(confo,3.)*lapse_bh*mass/lapse/rr/rr
174 - pow(confo,3.)*(divshf - 3.*lldllsh
175 + 2.*mass*lapse_bh%lapse_bh%llshift/rr/rr
176 + 4.*mass*pow(lapse_bh,3.)*lapse_rs
177 *(1.+3.*mass/rr)/rr/rr)
178 /6./lapse/lapse_bh ;
179
180 source_surf.std_spectral_base() ;
181 source_surf.annule_domain(0) ;
182 source_surf.raccord(1) ;
183
184 integ_s = mp_aff.integrale_surface(source_surf, radius_ah) ;
185
186 // Volume integral
187 // ---------------
188 Scalar lldlldco = lldconf.dsdr() ;
189 lldlldco.std_spectral_base() ;
190
191 Scalar tmp1 = 2.*mass*mass*pow(lapse_bh,4.)*confo/pow(rr,4.) ;
192 tmp1.std_spectral_base() ;
193 tmp1.inc_dzpuis(4) ;
194
195 Scalar tmp2 = 2.*mass*mass*pow(lapse_bh,6.)
196 *(1.+3.*mass/rr)*(1.+3.*mass/rr)*pow(confo,5.)/3./pow(rr,4.) ;
197 tmp2.std_spectral_base() ;
198 tmp2.inc_dzpuis(4) ;
199
200 Scalar tmp3 = 4.*mass*lapse_bh*lapse_bh*lldlldco/rr ;
201 tmp3.std_spectral_base() ;
202 tmp3.inc_dzpuis(1) ;
203
204 Scalar tmp4 = 2.*mass*pow(lapse_bh,4.)*lldconf
205 *(3.+8.*mass/rr)/rr/rr ;
206 tmp4.std_spectral_base() ;
207 tmp4.inc_dzpuis(2) ;
208
209 source_volm = 0.25 * taij_quad / pow(confo,7.) - tmp1 - tmp2
210 - tmp3 - tmp4 ;
211
212 source_volm.annule_domain(0) ;
213 source_volm.std_spectral_base() ;
214
215 integ_v = source_volm.integrale() ;
216
217 // ADM mass
218 // --------
219 madm = mass_bh + integ_s / qpig + integ_v / qpig ;
220
221 // Another ADM mass
222 // ----------------
223 double mmm = mass_bh
224 - 2.*(mp_aff.integrale_surface_infini(lldconf))/qpig ;
225
226 cout << "Another ADM mass : " << mmm / msol << endl ;
227 */
228 }
229 else { // Isotropic coordinates with the maximal slicing
230
231 Scalar divshf(mp) ;
232 divshf = shift(1).deriv(1) + shift(2).deriv(2)
233 + shift(3).deriv(3) ;
234 divshf.std_spectral_base() ;
235 divshf.dec_dzpuis(2) ;
236
237 Scalar llshift(mp) ;
238 llshift = ll(1)%shift(1) + ll(2)%shift(2) + ll(3)%shift(3) ;
239 llshift.std_spectral_base() ;
240
241 Scalar lldllsh = llshift.dsdr() ;
242 lldllsh.std_spectral_base() ;
243 lldllsh.dec_dzpuis(2) ;
244
245 // Surface integral
246 // ----------------
247 source_surf = confo/rr
248 - pow(confo,4.) * (divshf - 3.*lldllsh) / lapconf / 6. ;
249
250 source_surf.std_spectral_base() ;
251 source_surf.annule_domain(0) ;
252 source_surf.raccord(1) ;
253
254 integ_s = mp_aff.integrale_surface(source_surf, radius_ah) ;
255
256 // Volume integral
257 // ---------------
258 source_volm = 0.25 * taij_quad / pow(confo,7.) ;
259
260 source_volm.std_spectral_base() ;
261 source_volm.annule_domain(0) ;
262
263 integ_v = source_volm.integrale() ;
264
265 // ADM mass
266 // --------
267 madm = integ_s / qpig + integ_v / qpig ;
268
269 // Another ADM mass
270 // ----------------
271 double mmm = - 2.*(mp_aff.integrale_surface_infini(lldconf))/qpig ;
272
273 cout << "Another ADM mass : " << mmm / msol << endl ;
274
275 }
276
277 p_mass_adm = new double( madm ) ;
278
279 }
280
281 return *p_mass_adm ;
282
283}
284
285 //-----------------------------//
286 // Komar mass //
287 //-----------------------------//
288
289double Black_hole::mass_kom() const {
290
291 // Fundamental constants and units
292 // -------------------------------
293 using namespace Unites ;
294
295 if (p_mass_kom == 0x0) { // a new computation is required
296
297 double mkom ;
298 double integ_s ;
299 double integ_v ;
300
301 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
302 Map_af mp_aff(mp) ;
303
304 Scalar source_surf(mp) ;
305 source_surf.set_etat_qcq() ;
306 Scalar source_volm(mp) ;
307 source_volm.set_etat_qcq() ;
308
309 Scalar rr(mp) ;
310 rr = mp.r ;
311 rr.std_spectral_base() ;
312 Scalar st(mp) ;
313 st = mp.sint ;
314 st.std_spectral_base() ;
315 Scalar ct(mp) ;
316 ct = mp.cost ;
317 ct.std_spectral_base() ;
318 Scalar sp(mp) ;
319 sp = mp.sinp ;
320 sp.std_spectral_base() ;
321 Scalar cp(mp) ;
322 cp = mp.cosp ;
323 cp.std_spectral_base() ;
324
325 Vector ll(mp, CON, mp.get_bvect_cart()) ;
326 ll.set_etat_qcq() ;
327 ll.set(1) = st * cp ;
328 ll.set(2) = st * sp ;
329 ll.set(3) = ct ;
330 ll.std_spectral_base() ;
331
332 Vector dlcnf(mp, CON, mp.get_bvect_cart()) ;
333 dlcnf.set_etat_qcq() ;
334 for (int i=1; i<=3; i++)
335 dlcnf.set(i) = confo.deriv(i) / confo ;
336
337 dlcnf.std_spectral_base() ;
338
339 if (kerrschild) {
340
341 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
342 abort() ;
343 /*
344 Scalar llshift(mp) ;
345 llshift = ll(1)%shift_rs(1) + ll(2)%shift_rs(2)
346 + ll(3)%shift_rs(3) ;
347 llshift.std_spectral_base() ;
348
349 Vector dlap(mp, CON, mp.get_bvect_cart()) ;
350 dlap.set_etat_qcq() ;
351
352 for (int i=1; i<=3; i++)
353 dlap.set(i) = lapse_rs.deriv(i) ;
354
355 dlap.std_spectral_base() ;
356
357 double mass = ggrav * mass_bh ;
358
359 // Surface integral
360 // ----------------
361 Scalar lldlap_bh(mp) ;
362 lldlap_bh = pow(lapse_bh,3.) * mass / rr / rr ;
363 lldlap_bh.std_spectral_base() ;
364 lldlap_bh.annule_domain(0) ;
365 lldlap_bh.inc_dzpuis(2) ;
366
367 Scalar lldlap_rs = lapse_rs.dsdr() ;
368 lldlap_rs.std_spectral_base() ;
369
370 source_surf = lldlap_rs + lldlap_bh ;
371 source_surf.std_spectral_base() ;
372 source_surf.dec_dzpuis(2) ;
373 source_surf.annule_domain(0) ;
374 source_surf.raccord(1) ;
375
376 integ_s = mp_aff.integrale_surface(source_surf, radius_ah) ;
377
378 // Volume integral
379 // ---------------
380 Scalar lldlldlap = lldlap_rs.dsdr() ;
381 lldlldlap.std_spectral_base() ;
382
383 Scalar lldlcnf = lconfo.dsdr() ;
384 lldlcnf.std_spectral_base() ;
385
386 Scalar tmp1(mp) ;
387 tmp1 = dlap(1)%dlcnf(1) + dlap(2)%dlcnf(2) + dlap(3)%dlcnf(3)
388 -2.*mass*lapse_bh%lapse_bh%lldlap_rs%lldlcnf/rr ;
389 tmp1.std_spectral_base() ;
390
391 Scalar tmp2(mp) ;
392 tmp2 = 4.*mass*mass*pow(lapse_bh,6.)*(1.+3.*mass/rr)*(1.+3.*mass/rr)
393 *lapse_rs*pow(confo,4.)/3./pow(rr,4.) ;
394 tmp2.std_spectral_base() ;
395 tmp2.inc_dzpuis(4) ;
396
397 Scalar tmp3(mp) ;
398 tmp3 = -2.*mass*pow(lapse_bh,5.)*llshift
399 *(2.+10.*mass/rr+9.*mass*mass/rr/rr)*pow(confo,4.)/pow(rr,3.) ;
400 tmp3.std_spectral_base() ;
401 tmp3.inc_dzpuis(4) ;
402
403 Scalar tmp4(mp) ;
404 tmp4 = 2.*mass*lapse_bh%lapse_bh%lldlldlap/rr ;
405 tmp4.std_spectral_base() ;
406 tmp4.inc_dzpuis(1) ;
407
408 Scalar tmp5(mp) ;
409 tmp5 = mass*pow(lapse_bh,4.)*lldlap_rs*(3.+8.*mass/rr)/rr/rr ;
410 tmp5.std_spectral_base() ;
411 tmp5.inc_dzpuis(2) ;
412
413 Scalar tmp6(mp) ;
414 tmp6 = -2.*pow(lapse_bh,5.)*mass*lldlcnf/rr/rr ;
415 tmp6.std_spectral_base() ;
416 tmp6.inc_dzpuis(2) ;
417
418 Scalar tmp_bh(mp) ;
419 tmp_bh = -pow(lapse_bh,7.)*mass*mass
420 *( 4.*(5.+24.*mass/rr+18.*mass*mass/rr/rr)*pow(confo,4.)/3.
421 + 1. - 6.*mass/rr) / pow(rr, 4.) ;
422 tmp_bh.std_spectral_base() ;
423 tmp_bh.inc_dzpuis(4) ;
424
425 source_volm = lapse % taij_quad / pow(confo,8.) - 2.*tmp1
426 + tmp2 + tmp3 + tmp4 + tmp5 + tmp6 + tmp_bh ;
427
428 source_volm.annule_domain(0) ;
429 source_volm.std_spectral_base() ;
430
431 integ_v = source_volm.integrale() ;
432
433 // Komar mass
434 // ----------
435 mkom = integ_s / qpig + integ_v / qpig ;
436
437 // Another Komar mass
438 // ------------------
439 double mmm = (mp_aff.integrale_surface_infini(lldlap_rs+lldlap_bh))
440 / qpig ;
441
442 cout << "Another Komar mass : " << mmm / msol << endl ;
443 */
444 }
445 else { // Isotropic coordinates with the maximal slicing
446
447 // Surface integral
448 // ----------------
449 Scalar lldlap = lapconf.dsdr() / confo
450 - lapconf * confo.dsdr() / confo / confo ;
451 lldlap.std_spectral_base() ;
452
453 source_surf = lldlap ;
454
455 source_surf.std_spectral_base() ;
456 source_surf.dec_dzpuis(2) ;
457 source_surf.annule_domain(0) ;
458 source_surf.raccord(1) ;
459
460 integ_s = mp_aff.integrale_surface(source_surf, radius_ah) ;
461
462 // Volume integral
463 // ---------------
464 Vector dlap(mp, CON, mp.get_bvect_cart()) ;
465 dlap.set_etat_qcq() ;
466
467 for (int i=1; i<=3; i++)
468 dlap.set(i) = lapconf.deriv(i) / confo
469 - lapconf * confo.deriv(i) / confo / confo ;
470
471 dlap.std_spectral_base() ;
472
473 source_volm = lapconf % taij_quad / pow(confo,9.)
474 - 2.*(dlap(1)%dlcnf(1) + dlap(2)%dlcnf(2) + dlap(3)%dlcnf(3)) ;
475
476 source_volm.std_spectral_base() ;
477 source_volm.annule_domain(0) ;
478
479 integ_v = source_volm.integrale() ;
480
481 // Komar mass
482 // ----------
483 mkom = integ_s / qpig + integ_v / qpig ;
484
485 // Another Komar mass
486 double mmm = (mp_aff.integrale_surface_infini(lldlap)) / qpig ;
487
488 cout << "Another Komar mass : " << mmm / msol << endl ;
489
490 }
491
492 p_mass_kom = new double( mkom ) ;
493
494 }
495
496 return *p_mass_kom ;
497
498}
499
500
501 //------------------------------------//
502 // Apparent horizon //
503 //------------------------------------//
504
505double Black_hole::rad_ah() const {
506
507 if (p_rad_ah == 0x0) { // a new computation is required
508
509 Scalar rr(mp) ;
510 rr = mp.r ;
511 rr.std_spectral_base() ;
512
513 double rad = rr.val_grid_point(1,0,0,0) ;
514
515 p_rad_ah = new double( rad ) ;
516
517 }
518
519 return *p_rad_ah ;
520
521}
522
523
524 //-------------------------------------------//
525 // Quasi-local spin angular momentum //
526 //-------------------------------------------//
527
528double Black_hole::spin_am_bh(bool bclapconf_nd, bool bclapconf_fs,
529 const Tbl& xi_i, const double& phi_i,
530 const double& theta_i, const int& nrk_phi,
531 const int& nrk_theta) const {
532
533 // Fundamental constants and units
534 // -------------------------------
535 using namespace Unites ;
536
537 if (p_spin_am_bh == 0x0) { // a new computation is required
538
539 Scalar st(mp) ;
540 st = mp.sint ;
541 st.std_spectral_base() ;
542 Scalar ct(mp) ;
543 ct = mp.cost ;
544 ct.std_spectral_base() ;
545 Scalar sp(mp) ;
546 sp = mp.sinp ;
547 sp.std_spectral_base() ;
548 Scalar cp(mp) ;
549 cp = mp.cosp ;
550 cp.std_spectral_base() ;
551
552 Vector ll(mp, CON, mp.get_bvect_cart()) ;
553 ll.set_etat_qcq() ;
554 ll.set(1) = st % cp ;
555 ll.set(2) = st % sp ;
556 ll.set(3) = ct ;
557 ll.std_spectral_base() ;
558
559 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
560
561 if (kerrschild) {
562
563 cout << "Not yet prepared!!!" << endl ;
564 abort() ;
565
566 }
567 else { // Isotropic coordinates
568
569 // Killing vector of the spherical components
570 Vector killing_spher(mp, COV, mp.get_bvect_spher()) ;
571 killing_spher.set_etat_qcq() ;
572 killing_spher = killing_vect_bh(xi_i, phi_i, theta_i,
573 nrk_phi, nrk_theta) ;
574 killing_spher.std_spectral_base() ;
575
576 killing_spher.set(2) = confo * confo * radius_ah * killing_spher(2) ;
577 killing_spher.set(3) = confo * confo * radius_ah * killing_spher(3) ;
578 // killing_spher(3) is already divided by sin(theta)
579 killing_spher.std_spectral_base() ;
580
581 // Killing vector of the Cartesian components
582 Vector killing(mp, COV, mp.get_bvect_cart()) ;
583 killing.set_etat_qcq() ;
584
585 killing.set(1) = (killing_spher(2) * ct * cp - killing_spher(3) * sp)
586 / radius_ah ;
587 killing.set(2) = (killing_spher(2) * ct * sp + killing_spher(3) * cp)
588 / radius_ah ;
589 killing.set(3) = - killing_spher(2) * st / radius_ah ;
590 killing.std_spectral_base() ;
591
592 // Surface integral <- dzpuis should be 0
593 // --------------------------------------
594 // Source terms in the surface integral
595 Scalar source_1(mp) ;
596 source_1 = (ll(1) * (taij(1,1) * killing(1)
597 + taij(1,2) * killing(2)
598 + taij(1,3) * killing(3))
599 + ll(2) * (taij(2,1) * killing(1)
600 + taij(2,2) * killing(2)
601 + taij(2,3) * killing(3))
602 + ll(3) * (taij(3,1) * killing(1)
603 + taij(3,2) * killing(2)
604 + taij(3,3) * killing(3)))
605 / pow(confo, 4.) ;
606 source_1.std_spectral_base() ;
607 source_1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
608
609 Scalar source_surf(mp) ;
610 source_surf = source_1 ;
611 source_surf.std_spectral_base() ;
612 source_surf.annule_domain(0) ;
613 source_surf.raccord(1) ;
614
615 Map_af mp_aff(mp) ;
616
617 double spin = mp_aff.integrale_surface(source_surf, radius_ah) ;
618 double spin_angmom = 0.5 * spin / qpig ;
619
620 p_spin_am_bh = new double( spin_angmom ) ;
621
622 }
623
624 }
625
626 return *p_spin_am_bh ;
627
628}
629
630 //------------------------------------------//
631 // Total angular momentum //
632 //------------------------------------------//
633
635
636 // Fundamental constants and units
637 // -------------------------------
638 using namespace Unites ;
639
640 if (p_angu_mom_bh == 0x0) { // a new computation is required
641
642 p_angu_mom_bh = new Tbl(3) ;
643 p_angu_mom_bh->annule_hard() ; // fills the double array with zeros
644
645 double integ_bh_s_x ;
646 double integ_bh_s_y ;
647 double integ_bh_s_z ;
648
649 double radius_ah = mp.val_r(1,-1.,M_PI/2.,0.) ;
650 Map_af mp_aff(mp) ;
651
652 Scalar xx(mp) ;
653 xx = mp.x ;
654 xx.std_spectral_base() ;
655 Scalar yy(mp) ;
656 yy = mp.y ;
657 yy.std_spectral_base() ;
658 Scalar zz(mp) ;
659 zz = mp.z ;
660 zz.std_spectral_base() ;
661
662 Scalar st(mp) ;
663 st = mp.sint ;
664 st.std_spectral_base() ;
665 Scalar ct(mp) ;
666 ct = mp.cost ;
667 ct.std_spectral_base() ;
668 Scalar sp(mp) ;
669 sp = mp.sinp ;
670 sp.std_spectral_base() ;
671 Scalar cp(mp) ;
672 cp = mp.cosp ;
673 cp.std_spectral_base() ;
674
675 Vector ll(mp, CON, mp.get_bvect_cart()) ;
676 ll.set_etat_qcq() ;
677 ll.set(1) = st % cp ;
678 ll.set(2) = st % sp ;
679 ll.set(3) = ct ;
680 ll.std_spectral_base() ;
681
682 Vector jbh_x(mp, CON, mp.get_bvect_cart()) ;
683 jbh_x.set_etat_qcq() ;
684 for (int i=1; i<=3; i++)
685 jbh_x.set(i) = yy * taij(3,i) - zz * taij(2,i) ;
686
687 jbh_x.std_spectral_base() ;
688
689 Vector jbh_y(mp, CON, mp.get_bvect_cart()) ;
690 jbh_y.set_etat_qcq() ;
691 for (int i=1; i<=3; i++)
692 jbh_y.set(i) = zz * taij(1,i) - xx * taij(3,i) ;
693
694 jbh_y.std_spectral_base() ;
695
696 Vector jbh_z(mp, CON, mp.get_bvect_cart()) ;
697 jbh_z.set_etat_qcq() ;
698 for (int i=1; i<=3; i++)
699 jbh_z.set(i) = xx * taij(2,i) - yy * taij(1,i) ;
700
701 jbh_z.std_spectral_base() ;
702
703 /*
704 if (kerrschild) {
705
706 cout << "Not yet prepared!!!" << endl ;
707 abort() ;
708
709 }
710 else { // Isotropic coordinates
711
712 // Sets C/M^2 for each case of the lapse boundary condition
713 // --------------------------------------------------------
714 double cc ;
715
716 if (bclapconf_nd) { // Neumann boundary condition
717 if (bclapconf_fs) { // First condition
718 // d(\alpha \psi)/dr = 0
719 // ---------------------
720 cc = 2. * (sqrt(13.) - 1.) / 3. ;
721 }
722 else { // Second condition
723 // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
724 // -----------------------------------------
725 cc = 4. / 3. ;
726 }
727 }
728 else { // Dirichlet boundary condition
729 if (bclapconf_fs) { // First condition
730 // (\alpha \psi) = 1/2
731 // -------------------
732 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
733 abort() ;
734 }
735 else { // Second condition
736 // (\alpha \psi) = 1/sqrt(2.) \psi_KS
737 // ----------------------------------
738 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
739 abort() ;
740 }
741 }
742
743 Scalar r_are(mp) ;
744 r_are = r_coord(bclapconf_nd, bclapconf_fs) ;
745 r_are.std_spectral_base() ;
746
747 }
748 */
749
750 //-------------//
751 // X component //
752 //-------------//
753
754 //-------------------------------------
755 // Integration over the BH map
756 //-------------------------------------
757
758 // Surface integral <- dzpuis should be 0
759 // --------------------------------------
760 Scalar sou_bh_sx(mp) ;
761 sou_bh_sx = jbh_x(1)%ll(1) + jbh_x(2)%ll(2) + jbh_x(3)%ll(3) ;
762 sou_bh_sx.std_spectral_base() ;
763 sou_bh_sx.dec_dzpuis(2) ; // dzpuis : 2 -> 0
764 sou_bh_sx.annule_domain(0) ;
765 sou_bh_sx.raccord(1) ;
766
767 integ_bh_s_x = mp_aff.integrale_surface(sou_bh_sx, radius_ah) ;
768
769 p_angu_mom_bh->set(0) = 0.5 * integ_bh_s_x / qpig ;
770
771 //-------------//
772 // Y component //
773 //-------------//
774
775 //-------------------------------------
776 // Integration over the BH map
777 //-------------------------------------
778
779 // Surface integral <- dzpuis should be 0
780 // --------------------------------------
781 Scalar sou_bh_sy(mp) ;
782 sou_bh_sy = jbh_y(1)%ll(1) + jbh_y(2)%ll(2) + jbh_y(3)%ll(3) ;
783 sou_bh_sy.std_spectral_base() ;
784 sou_bh_sy.dec_dzpuis(2) ; // dzpuis : 2 -> 0
785 sou_bh_sy.annule_domain(0) ;
786 sou_bh_sy.raccord(1) ;
787
788 integ_bh_s_y = mp_aff.integrale_surface(sou_bh_sy, radius_ah) ;
789
790 p_angu_mom_bh->set(1) = 0.5 * integ_bh_s_y / qpig ;
791
792 //-------------//
793 // Z component //
794 //-------------//
795
796 //-------------------------------------
797 // Integration over the BH map
798 //-------------------------------------
799
800 // Surface integral <- dzpuis should be 0
801 // --------------------------------------
802 Scalar sou_bh_sz(mp) ;
803 sou_bh_sz = jbh_z(1)%ll(1) + jbh_z(2)%ll(2) + jbh_z(3)%ll(3) ;
804 sou_bh_sz.std_spectral_base() ;
805 sou_bh_sz.dec_dzpuis(2) ; // dzpuis : 2 -> 0
806 sou_bh_sz.annule_domain(0) ;
807 sou_bh_sz.raccord(1) ;
808
809 integ_bh_s_z = mp_aff.integrale_surface(sou_bh_sz, radius_ah) ;
810
811 p_angu_mom_bh->set(2) = 0.5 * integ_bh_s_z / qpig ;
812
813 }
814
815 return *p_angu_mom_bh ;
816
817}
818}
virtual double mass_kom() const
Komar mass.
double spin_am_bh(bool bclapconf_nd, bool bclapconf_fs, const Tbl &xi_i, const double &phi_i, const double &theta_i, const int &nrk_phi, const int &nrk_theta) const
Spin angular momentum.
Scalar taij_quad
Part of the scalar generated by .
Definition blackhole.h:135
Sym_tensor taij
Trace of the extrinsic curvature.
Definition blackhole.h:127
virtual double mass_irr() const
Irreducible mass of the black hole.
double * p_spin_am_bh
Radius of the apparent horizon.
Definition blackhole.h:159
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
virtual double rad_ah() const
Radius of the apparent horizon.
Tbl * p_angu_mom_bh
Spin angular momentum.
Definition blackhole.h:162
Vector killing_vect_bh(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.
virtual double mass_adm() const
ADM mass.
double * p_rad_ah
Komar mass.
Definition blackhole.h:157
double * p_mass_kom
ADM mass.
Definition blackhole.h:155
const Tbl & angu_mom_bh() const
Total angular momentum.
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
double * p_mass_adm
Irreducible mass of the black hole.
Definition blackhole.h:153
double * p_mass_irr
Conformal metric .
Definition blackhole.h:151
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 .
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 y
y coordinate centered on the grid
Definition map.h:727
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 z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
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.
const Scalar & deriv(int i) const
Returns of *this , where .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
double integrale() const
Computes the integral over all space 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
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
const Scalar & dsdr() const
Returns of *this .
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
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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.