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 
71 namespace Lorene {
72 double 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 
106 double 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 
292 double 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 
508 double 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 
531 double 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 
637 const 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 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
virtual double mass_irr() const
Irreducible mass of the black hole.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
const Tbl & angu_mom_bh() const
Total angular momentum.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Scalar taij_quad
Part of the scalar generated by .
Definition: blackhole.h:135
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
double * p_mass_adm
Irreducible mass of the black hole.
Definition: blackhole.h:153
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
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
Standard units of space, time and mass.
double * p_mass_kom
ADM mass.
Definition: blackhole.h:155
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
double integrale() const
Computes the integral over all space of *this .
Definition: scalar_integ.C:64
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 double mass_adm() const
ADM mass.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Tensor field of valence 1.
Definition: vector.h:188
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
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Coord sint
Definition: map.h:739
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.
double * p_mass_irr
Conformal metric .
Definition: blackhole.h:151
double * p_rad_ah
Komar mass.
Definition: blackhole.h:157
Coord sinp
Definition: map.h:741
virtual double rad_ah() const
Radius of the apparent horizon.
const Scalar & deriv(int i) const
Returns of *this , where .
Definition: scalar_deriv.C:359
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
double * p_spin_am_bh
Radius of the apparent horizon.
Definition: blackhole.h:159
virtual double mass_kom() const
Komar mass.
Tbl * p_angu_mom_bh
Spin angular momentum.
Definition: blackhole.h:162
Affine radial mapping.
Definition: map.h:2048
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
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Coord y
y coordinate centered on the grid
Definition: map.h:745
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Coord cosp
Definition: map.h:742
Sym_tensor taij
Trace of the extrinsic curvature.
Definition: blackhole.h:127
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Coord x
x coordinate centered on the grid
Definition: map.h:744
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...
Basic array class.
Definition: tbl.h:164
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition: blackhole.h:97
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.
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Coord z
z coordinate centered on the grid
Definition: map.h:746
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...
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord cost
Definition: map.h:740