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