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 char 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 
68 namespace Lorene {
69 double 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 
103 double 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 
289 double 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 
505 double 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 
528 double 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 
634 const Tbl& Black_hole::angu_mom_bh() const {
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 }
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:666
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:481
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:220
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:783
Lorene prototypes.
Definition: app_hor.h:64
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:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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:61
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 double mass_adm() const
ADM mass.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
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:637
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:721
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:723
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:348
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:2027
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:791
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Coord y
y coordinate centered on the grid
Definition: map.h:727
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Coord cosp
Definition: map.h:724
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:726
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:161
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:372
Coord z
z coordinate centered on the grid
Definition: map.h:728
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:718
Coord cost
Definition: map.h:722