LORENE
bin_bhns_global.C
1 /*
2  * Methods of class Bin_bhns to compute global quantities
3  *
4  * (see file bin_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005-2007 Keisuke Taniguchi
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char bin_bhns_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_global.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
29 
30 /*
31  * $Id: bin_bhns_global.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
32  * $Log: bin_bhns_global.C,v $
33  * Revision 1.4 2014/10/13 08:52:41 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.3 2014/10/06 15:13:00 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.2 2008/05/15 18:59:27 k_taniguchi
40  * Introduction of new global quantities.
41  *
42  * Revision 1.1 2007/06/22 01:09:31 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_global.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "bin_bhns.h"
58 #include "blackhole.h"
59 #include "unites.h"
60 #include "utilitaires.h"
61 #include "nbr_spx.h"
62 
63  //----------------------------//
64  // ADM mass //
65  //----------------------------//
66 
67 namespace Lorene {
69 
70  // Fundamental constants and units
71  // -------------------------------
72  using namespace Unites ;
73 
74  if (p_mass_adm_bhns_surf == 0x0) { // a new computation is required
75 
76  double madm ;
77 
78  const Map& mp_bh = hole.get_mp() ;
79  const Map& mp_ns = star.get_mp() ;
80 
81  Map_af mp_aff(mp_bh) ;
82  Map_af mp_ns_aff(mp_ns) ;
83 
84  Scalar rr(mp_bh) ;
85  rr = mp_bh.r ;
86  rr.std_spectral_base() ;
87 
88  double mass = ggrav * hole.get_mass_bh() ;
89 
90  if (hole.is_kerrschild()) {
91 
92  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
93  abort() ;
94 
95  }
96  else { // Isotropic coordinates with the maximal slicing
97 
98  //-------------------------------------
99  // Integration over the BH map
100  //-------------------------------------
101 
102  // Sets C/M^2 for each case of the lapse boundary condition
103  // --------------------------------------------------------
104  double cc ;
105 
106  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
107  if (hole.has_bc_lapconf_fs()) { // First condition
108  // d(\alpha \psi)/dr = 0
109  // ---------------------
110  cc = 2. * (sqrt(13.) - 1.) / 3. ;
111  }
112  else { // Second condition
113  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
114  // -----------------------------------------
115  cc = 4. / 3. ;
116  }
117  }
118  else { // Dirichlet boundary condition
119  if (hole.has_bc_lapconf_fs()) { // First condition
120  // (\alpha \psi) = 1/2
121  // -------------------
122  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
123  abort() ;
124  }
125  else { // Second condition
126  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
127  // ----------------------------------
128  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
129  abort() ;
130  // cc = 2. * sqrt(2.) ;
131  }
132  }
133 
134  Scalar r_are(mp_bh) ;
137  r_are.std_spectral_base() ;
138 
139  // ADM mass by surface integral at infinity : dzpuis should be 2
140  // ----------------------------------------
141  const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
142 
143  Scalar lldconf_iso = confo_bh_auto_rs.dsdr() ; // dzpuis = 2
144  lldconf_iso.std_spectral_base() ;
145 
146  Scalar anoth(mp_bh) ;
147  anoth = 0.5 * sqrt(r_are)
148  * (sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
149  - 1.) / rr ;
150  anoth.std_spectral_base() ;
151  anoth.annule_domain(0) ;
152  anoth.raccord(1) ;
153  anoth.inc_dzpuis(2) ;
154 
155  const Scalar& confo_ns_auto = star.get_confo_auto() ;
156 
157  Scalar lldconf_ns = confo_ns_auto.dsdr() ; // dzpuis = 2
158  lldconf_ns.std_spectral_base() ;
159 
160  madm =
161  - 2.*(mp_aff.integrale_surface_infini(lldconf_iso+anoth))/qpig
162  - 2.*(mp_ns_aff.integrale_surface_infini(lldconf_ns))/qpig ;
163 
164  cout << "ADM mass (surface) : " << madm / msol << " [Mo]"
165  << endl ;
166 
167  }
168 
169  p_mass_adm_bhns_surf = new double( madm ) ;
170 
171  }
172 
173  return *p_mass_adm_bhns_surf ;
174 
175 }
176 
177 
178  //----------------------------//
179  // ADM mass //
180  //----------------------------//
181 
182 double Bin_bhns::mass_adm_bhns_vol() const {
183 
184  // Fundamental constants and units
185  // -------------------------------
186  using namespace Unites ;
187 
188  if (p_mass_adm_bhns_vol == 0x0) { // a new computation is required
189 
190  double madm ;
191  double integ_bh_s ;
192  double integ_bh_v ;
193  double integ_ns_v ;
194 
195  const Map& mp_bh = hole.get_mp() ;
196  const Map& mp_ns = star.get_mp() ;
197 
198  double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
199  Map_af mp_aff(mp_bh) ;
200 
201  Map_af mp_ns_aff(mp_ns) ;
202 
203  Scalar source_bh_surf(mp_bh) ;
204  source_bh_surf.set_etat_qcq() ;
205 
206  Scalar source_bh_volm(mp_bh) ;
207  source_bh_volm.set_etat_qcq() ;
208 
209  Scalar source_ns_volm(mp_ns) ;
210  source_ns_volm.set_etat_qcq() ;
211 
212  Scalar rr(mp_bh) ;
213  rr = mp_bh.r ;
214  rr.std_spectral_base() ;
215  Scalar st(mp_bh) ;
216  st = mp_bh.sint ;
217  st.std_spectral_base() ;
218  Scalar ct(mp_bh) ;
219  ct = mp_bh.cost ;
220  ct.std_spectral_base() ;
221  Scalar sp(mp_bh) ;
222  sp = mp_bh.sinp ;
223  sp.std_spectral_base() ;
224  Scalar cp(mp_bh) ;
225  cp = mp_bh.cosp ;
226  cp.std_spectral_base() ;
227 
228  Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
229  ll.set_etat_qcq() ;
230  ll.set(1) = st % cp ;
231  ll.set(2) = st % sp ;
232  ll.set(3) = ct ;
233  ll.std_spectral_base() ;
234 
235  const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
236  const Vector& shift_comp = hole.get_shift_comp() ;
237  const Tensor& dshift_comp = hole.get_d_shift_comp() ;
238 
239  Scalar divshift(mp_bh) ; // dzpuis = 2
240  divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
241  + shift_auto_rs(3).deriv(3) + dshift_comp(1,1)
242  + dshift_comp(2,2) + dshift_comp(3,3) ;
243  divshift.std_spectral_base() ;
244 
245  Scalar llshift(mp_bh) ; // dzpuis = 0
246  llshift = ll(1) % (shift_auto_rs(1) + shift_comp(1))
247  + ll(2) % (shift_auto_rs(2) + shift_comp(2))
248  + ll(3) % (shift_auto_rs(3) + shift_comp(3)) ;
249  llshift.std_spectral_base() ;
250 
251  Scalar llshift_auto(mp_bh) ; // dzpuis = 0
252  llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
253  + ll(3)%shift_auto_rs(3) ;
254  llshift_auto.std_spectral_base() ;
255 
256  Scalar lldllsh = llshift_auto.dsdr()
257  + ll(1) * (ll(1) % dshift_comp(1,1) + ll(2) % dshift_comp(1,2)
258  + ll(3) % dshift_comp(1,3))
259  + ll(2) * (ll(1) % dshift_comp(2,1) + ll(2) % dshift_comp(2,2)
260  + ll(3) % dshift_comp(2,3))
261  + ll(3) * (ll(1) % dshift_comp(3,1) + ll(2) % dshift_comp(3,2)
262  + ll(3) % dshift_comp(3,3)) ; // dzpuis = 2
263  lldllsh.std_spectral_base() ;
264 
265  const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
266  const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
267  const Scalar& lapconf_bh_comp = hole.get_lapconf_comp() ;
268  const Scalar& confo_bh = hole.get_confo_tot() ;
269  const Scalar& confo_bh_auto = hole.get_confo_auto() ;
270  const Scalar& confo_bh_comp = hole.get_confo_comp() ;
271  const Vector& dconfo_bh_comp = hole.get_d_confo_comp() ;
272  const Scalar& taij_quad_tot_rs = hole.get_taij_quad_tot_rs() ;
273  const Scalar& taij_quad_tot_rot = hole.get_taij_quad_tot_rot() ;
274 
275  const Scalar& taij_quad_auto_bh = hole.get_taij_quad_auto() ;
276  const Scalar& taij_quad_comp_bh = hole.get_taij_quad_comp() ;
277 
278  const Scalar& confo_ns = star.get_confo_tot() ;
279  const Scalar& confo_ns_auto = star.get_confo_auto() ;
280  const Scalar& ener_euler = star.get_ener_euler() ;
281 
282  const Scalar& taij_quad_auto_ns = star.get_taij_quad_auto() ;
283 
284  Scalar lldconf = confo_bh_auto.dsdr() + ll(1)%dconfo_bh_comp(1)
285  + ll(2)%dconfo_bh_comp(2) + ll(3)%dconfo_bh_comp(3) ; // dzpuis = 2
286  lldconf.std_spectral_base() ;
287 
288  double mass = ggrav * hole.get_mass_bh() ;
289 
290  if (hole.is_kerrschild()) {
291 
292  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
293  abort() ;
294 
295  /*
296  Scalar lap_bh(mp_bh) ;
297  lap_bh = 1./sqrt(1.+2.*mass/rr) ;
298  lap_bh.std_spectral_base() ;
299 
300  Scalar lap_bh2(mp_bh) ;
301  lap_bh2 = 1./(1.+2.*mass/rr) ;
302  lap_bh2.std_spectral_base() ;
303 
304  Scalar omelld(mp_bh) ;
305  omelld = omega * (ll(2) * (mp_bh.get_ori_x() - x_rot)
306  - ll(1) * (mp_bh.get_ori_y() - y_rot)) ;
307  omelld.std_spectral_base() ;
308 
309  Scalar lldlldconf(mp_bh) ; // dzpuis = 3
310  lldlldconf = lldconf.dsdr() ;
311  lldlldconf.std_spectral_base() ;
312 
313  //-------------------------------------
314  // Integration over the BH map
315  //-------------------------------------
316 
317  // Surface integral <- dzpuis should be 0
318  // --------------------------------------
319  Scalar divshift_zero(divshift) ;
320  divshift_zero.dec_dzpuis(2) ;
321 
322  Scalar lldllsh_zero(lldllsh) ;
323  lldllsh_zero.dec_dzpuis(2) ;
324 
325  source_bh_surf = confo_bh
326  * (1. - 2.*mass*lap_bh*confo_bh*confo_bh/lapse_bh/rr) / rr
327  - pow(confo_bh, 3.)
328  * ( divshift_zero - 3.*lldllsh_zero
329  + 2. * lap_bh2 * mass * (llshift + omelld) / rr / rr
330  + 4.*mass*lap_bh2*lap_bh*(1.+3.*mass/rr)
331  *(lapse_bh_auto_rs + lapse_bh_comp)/rr/rr
332  ) / 6. / lap_bh / lapse_bh ;
333 
334  source_bh_surf.std_spectral_base() ;
335  source_bh_surf.annule_domain(0) ;
336  source_bh_surf.raccord(1) ;
337 
338  integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
339 
340  // Volume integral <- dzpuis should be 4
341  // -------------------------------------
342  Scalar sou_bh1(mp_bh) ;
343  sou_bh1 = 2.*lap_bh2*mass*lldlldconf/rr ;
344  sou_bh1.std_spectral_base() ;
345  sou_bh1.inc_dzpuis(1) ;
346 
347  Scalar sou_bh2(mp_bh) ;
348  sou_bh2 = lap_bh2*lap_bh2*mass*(3.+8.*mass/rr)*lldconf/rr/rr ;
349  sou_bh2.std_spectral_base() ;
350  sou_bh2.inc_dzpuis(2) ;
351 
352  Scalar sou_bh3(mp_bh) ;
353  sou_bh3 = pow(lap_bh2,3.)*mass*mass*confo_bh
354  * ( (1.-lap_bh2/lapse_bh/lapse_bh)
355  *(4.+12.*mass/rr+9.*mass*mass/rr/rr)*pow(confo_bh,4.)
356  +3.*(1.+2.*mass/rr)*(1.-pow(confo_bh,4.)) )
357  /3./pow(rr,4.) ;
358  sou_bh3.std_spectral_base() ;
359  sou_bh3.inc_dzpuis(4) ;
360 
361  source_bh_volm = 0.25 * (taij_quad_tot_rs + taij_quad_tot_rot)
362  / pow(confo_bh,7.)
363  - 2. * (sou_bh1 + sou_bh2 + sou_bh3) ;
364 
365  source_bh_volm.std_spectral_base() ;
366  source_bh_volm.annule_domain(0) ;
367 
368  integ_bh_v = source_bh_volm.integrale() ;
369 
370  //-------------------------------------
371  // Integration over the NS map
372  //-------------------------------------
373 
374  // Volume integral <- dzpuis should be 4
375  // -------------------------------------
376  source_ns_volm = pow(confo_ns, 5.) * ener_euler ;
377  source_ns_volm.std_spectral_base() ;
378  source_ns_volm.inc_dzpuis(4) ;
379 
380  integ_ns_v = source_ns_volm.integrale() ;
381 
382  cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
383  << " integ_bh_v : "
384  << integ_bh_v/ qpig / msol
385  << " integ_ns_v : " << integ_ns_v/ msol << endl ;
386 
387  //------------------
388  // ADM mass
389  //------------------
390  madm = hole.get_mass_bh()
391  + (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
392 
393  cout << "ADM mass : " << madm / msol << " [Mo]"
394  << endl ;
395 
396  // ADM mass by surface integral at infinity : dzpuis should be 2
397  // ----------------------------------------
398  double mmm = hole.get_mass_bh()
399  - 2.*(mp_aff.integrale_surface_infini(lldconf))/qpig ;
400 
401  cout << "Another ADM mass : " << mmm / msol << " [Mo]"
402  << endl ;
403  */
404  }
405  else { // Isotropic coordinates with the maximal slicing
406 
407  //-------------------------------------
408  // Integration over the BH map
409  //-------------------------------------
410 
411  // Sets C/M^2 for each case of the lapse boundary condition
412  // --------------------------------------------------------
413  double cc ;
414 
415  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
416  if (hole.has_bc_lapconf_fs()) { // First condition
417  // d(\alpha \psi)/dr = 0
418  // ---------------------
419  cc = 2. * (sqrt(13.) - 1.) / 3. ;
420  }
421  else { // Second condition
422  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
423  // -----------------------------------------
424  cc = 4. / 3. ;
425  }
426  }
427  else { // Dirichlet boundary condition
428  if (hole.has_bc_lapconf_fs()) { // First condition
429  // (\alpha \psi) = 1/2
430  // -------------------
431  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
432  abort() ;
433  }
434  else { // Second condition
435  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
436  // ----------------------------------
437  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
438  abort() ;
439  // cc = 2. * sqrt(2.) ;
440  }
441  }
442 
443  Scalar r_are(mp_bh) ;
446  r_are.std_spectral_base() ;
447 
448  // Surface integral <- dzpuis should be 0
449  // --------------------------------------
450  Scalar divshift_zero(divshift) ;
451  divshift_zero.dec_dzpuis(2) ;
452 
453  Scalar lldllsh_zero(lldllsh) ;
454  lldllsh_zero.dec_dzpuis(2) ;
455 
456  source_bh_surf = confo_bh / rr
457  - pow(confo_bh, 4.) * (divshift_zero - 3.*lldllsh_zero)
458  /6./lapconf_bh
459  - pow(confo_bh, 4.) * mass * mass * cc
460  * sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
461  / lapconf_bh / pow(r_are*rr,3.) ;
462 
463  source_bh_surf.std_spectral_base() ;
464  source_bh_surf.annule_domain(0) ;
465  source_bh_surf.raccord(1) ;
466 
467  integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
468 
469  // Volume integral <- dzpuis should be 4
470  // -------------------------------------
471  Scalar sou_bh1(mp_bh) ;
472  sou_bh1 = 1.5 * pow(confo_bh,7.) * pow(mass*mass*cc,2.)
473  * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
474  / lapconf_bh / lapconf_bh / pow(r_are*rr,6.) ;
475  sou_bh1.std_spectral_base() ;
476  sou_bh1.inc_dzpuis(4) ;
477 
478  source_bh_volm = 0.25 * taij_quad_auto_bh / pow(confo_bh,7.)
479  + 0.25 * taij_quad_comp_bh
480  * (pow(confo_bh/(confo_bh_comp+0.5),7.)
481  *pow((lapconf_bh_comp+0.5)/lapconf_bh,2.)
482  - 1.) / pow(confo_bh_comp+0.5,7.)
483  + sou_bh1 ;
484 
485  source_bh_volm.std_spectral_base() ;
486  source_bh_volm.annule_domain(0) ;
487 
488  integ_bh_v = source_bh_volm.integrale() ;
489 
490  //-------------------------------------
491  // Integration over the NS map
492  //-------------------------------------
493 
494  // Volume integral <- dzpuis should be 4
495  // -------------------------------------
496  Scalar sou_ns1(mp_ns) ;
497  sou_ns1 = pow(confo_ns, 5.) * ener_euler ;
498  sou_ns1.std_spectral_base() ;
499  sou_ns1.inc_dzpuis(4) ;
500 
501  source_ns_volm = 0.25 * taij_quad_auto_ns
502  / pow(confo_ns_auto+0.5, 7.) / qpig + sou_ns1 ;
503 
504  source_ns_volm.std_spectral_base() ;
505 
506  integ_ns_v = source_ns_volm.integrale() ;
507 
508  cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
509  << " integ_bh_v : "
510  << integ_bh_v/ qpig / msol
511  << " integ_ns_v : " << integ_ns_v/ msol << endl ;
512 
513  //------------------
514  // ADM mass
515  //------------------
516  madm = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
517 
518  cout << "ADM mass (volume) : " << madm / msol << " [Mo]"
519  << endl ;
520 
521  }
522 
523  p_mass_adm_bhns_vol = new double( madm ) ;
524 
525  }
526 
527  return *p_mass_adm_bhns_vol ;
528 
529 }
530 
531 
532 
533  //------------------------------//
534  // Komar mass //
535  //------------------------------//
536 
538 
539  // Fundamental constants and units
540  // -------------------------------
541  using namespace Unites ;
542 
543  if (p_mass_kom_bhns_surf == 0x0) { // a new computation is required
544 
545  double mkom ;
546 
547  const Map& mp_bh = hole.get_mp() ;
548  const Map& mp_ns = star.get_mp() ;
549 
550  Map_af mp_aff(mp_bh) ;
551  Map_af mp_ns_aff(mp_ns) ;
552 
553  Scalar rr(mp_bh) ;
554  rr = mp_bh.r ;
555  rr.std_spectral_base() ;
556 
557  double mass = ggrav * hole.get_mass_bh() ;
558 
559  if (hole.is_kerrschild()) {
560 
561  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
562  abort() ;
563 
564  }
565  else { // Isotropic coordinates with the maximal slicing
566 
567  //-------------------------------------
568  // Integration over the BH map
569  //-------------------------------------
570 
571  // Sets C/M^2 for each case of the lapse boundary condition
572  // --------------------------------------------------------
573  double cc ;
574 
575  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
576  if (hole.has_bc_lapconf_fs()) { // First condition
577  // d(\alpha \psi)/dr = 0
578  // ---------------------
579  cc = 2. * (sqrt(13.) - 1.) / 3. ;
580  }
581  else { // Second condition
582  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
583  // -----------------------------------------
584  cc = 4. / 3. ;
585  }
586  }
587  else { // Dirichlet boundary condition
588  if (hole.has_bc_lapconf_fs()) { // First condition
589  // (\alpha \psi) = 1/2
590  // -------------------
591  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
592  abort() ;
593  }
594  else { // Second condition
595  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
596  // ----------------------------------
597  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
598  abort() ;
599  // cc = 2. * sqrt(2.) ;
600  }
601  }
602 
603  Scalar r_are(mp_bh) ;
606  r_are.std_spectral_base() ;
607 
608  // Komar mass by surface integral at infinity : dzpuis should be 2
609  // ------------------------------------------
610  const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
611  const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
612 
613  Scalar lldlap_bh(mp_bh) ;
614  lldlap_bh = lapconf_bh_auto_rs.dsdr()
615  - confo_bh_auto_rs.dsdr() ; // dzpuis = 2
616  lldlap_bh.std_spectral_base() ;
617 
618  Scalar anoth(mp_bh) ;
619  anoth = sqrt(r_are) * (1. - 1.5 *cc*cc*pow(mass/r_are/rr,4.)
620  - sqrt(1. - 2.*mass/r_are/rr
621  + cc*cc*pow(mass/r_are/rr,4.))) / rr ;
622  anoth.std_spectral_base() ;
623  anoth.annule_domain(0) ;
624  anoth.raccord(1) ;
625  anoth.inc_dzpuis(2) ;
626 
627  const Scalar& lapconf_ns_auto = star.get_lapconf_auto() ;
628  const Scalar& confo_ns_auto = star.get_confo_auto() ;
629 
630  Scalar lldlap_ns(mp_ns) ;
631  lldlap_ns = lapconf_ns_auto.dsdr() - confo_ns_auto.dsdr() ;
632  lldlap_ns.std_spectral_base() ; // dzpuis = 2
633 
634  mkom =
635  (mp_aff.integrale_surface_infini(lldlap_bh+anoth))/qpig
636  + (mp_ns_aff.integrale_surface_infini(lldlap_ns))/qpig ;
637 
638  cout << "Komar mass (surface) : " << mkom / msol << " [Mo]"
639  << endl ;
640 
641  }
642 
643  p_mass_kom_bhns_surf = new double( mkom ) ;
644 
645  }
646 
647  return *p_mass_kom_bhns_surf ;
648 
649 }
650 
651 
652 
653  //------------------------------//
654  // Komar mass //
655  //------------------------------//
656 
657 double Bin_bhns::mass_kom_bhns_vol() const {
658 
659  // Fundamental constants and units
660  // -------------------------------
661  using namespace Unites ;
662 
663  if (p_mass_kom_bhns_vol == 0x0) { // a new computation is required
664 
665  double mkom ;
666  double integ_bh_s ;
667  double integ_bh_v ;
668  double integ_ns_v ;
669 
670  const Map& mp_bh = hole.get_mp() ;
671  const Map& mp_ns = star.get_mp() ;
672 
673  double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
674  Map_af mp_aff(mp_bh) ;
675 
676  Map_af mp_ns_aff(mp_ns) ;
677 
678  Scalar source_bh_surf(mp_bh) ;
679  source_bh_surf.set_etat_qcq() ;
680 
681  Scalar source_bh_volm(mp_bh) ;
682  source_bh_volm.set_etat_qcq() ;
683 
684  Scalar source_ns_volm(mp_ns) ;
685  source_ns_volm.set_etat_qcq() ;
686 
687  Scalar rr(mp_bh) ;
688  rr = mp_bh.r ;
689  rr.std_spectral_base() ;
690  Scalar st(mp_bh) ;
691  st = mp_bh.sint ;
692  st.std_spectral_base() ;
693  Scalar ct(mp_bh) ;
694  ct = mp_bh.cost ;
695  ct.std_spectral_base() ;
696  Scalar sp(mp_bh) ;
697  sp = mp_bh.sinp ;
698  sp.std_spectral_base() ;
699  Scalar cp(mp_bh) ;
700  cp = mp_bh.cosp ;
701  cp.std_spectral_base() ;
702 
703  Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
704  ll.set_etat_qcq() ;
705  ll.set(1) = st % cp ;
706  ll.set(2) = st % sp ;
707  ll.set(3) = ct ;
708  ll.std_spectral_base() ;
709 
710  const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
711  const Scalar& lapconf_bh_auto_rs = hole.get_lapconf_auto_rs() ;
712  const Scalar& lapconf_bh_comp = hole.get_lapconf_comp() ;
713  const Vector& dlapconf_bh_comp = hole.get_d_lapconf_comp() ;
714  const Scalar& confo_bh = hole.get_confo_tot() ;
715  const Scalar& confo_bh_auto_rs = hole.get_confo_auto_rs() ;
716  const Scalar& confo_bh_comp = hole.get_confo_comp() ;
717  const Vector& dconfo_bh_comp = hole.get_d_confo_comp() ;
718  const Scalar& taij_quad_tot_rs = hole.get_taij_quad_tot_rs() ;
719  const Scalar& taij_quad_tot_rot = hole.get_taij_quad_tot_rot() ;
720 
721  const Scalar& taij_quad_auto_bh = hole.get_taij_quad_auto() ;
722  const Scalar& taij_quad_comp_bh = hole.get_taij_quad_comp() ;
723 
724  const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
725  const Vector& shift_comp = hole.get_shift_comp() ;
726  const Tensor& dshift_comp = hole.get_d_shift_comp() ;
727 
728  Scalar divshift(mp_bh) ; // dzpuis = 2
729  divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
730  + shift_auto_rs(3).deriv(3) + dshift_comp(1,1)
731  + dshift_comp(2,2) + dshift_comp(3,3) ;
732  divshift.std_spectral_base() ;
733 
734  Scalar llshift_auto(mp_bh) ; // dzpuis = 0
735  llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
736  + ll(3)%shift_auto_rs(3) ;
737  llshift_auto.std_spectral_base() ;
738 
739  Scalar lldllsh = llshift_auto.dsdr()
740  + ll(1) * (ll(1) % dshift_comp(1,1) + ll(2) % dshift_comp(1,2)
741  + ll(3) % dshift_comp(1,3))
742  + ll(2) * (ll(1) % dshift_comp(2,1) + ll(2) % dshift_comp(2,2)
743  + ll(3) % dshift_comp(2,3))
744  + ll(3) * (ll(1) % dshift_comp(3,1) + ll(2) % dshift_comp(3,2)
745  + ll(3) % dshift_comp(3,3)) ; // dzpuis = 2
746  lldllsh.std_spectral_base() ;
747 
748  const Scalar& lapconf_ns = star.get_lapconf_tot() ;
749  const Scalar& ener_euler = star.get_ener_euler() ;
750  const Scalar& s_euler = star.get_s_euler() ;
751 
752  const Scalar& confo_ns = star.get_confo_tot() ;
753  const Scalar& lapconf_ns_auto = star.get_lapconf_auto() ;
754  const Scalar& confo_ns_auto = star.get_confo_auto() ;
755  const Vector& dconfo_ns_comp = star.get_d_confo_comp() ;
756  const Scalar& taij_quad_auto_ns = star.get_taij_quad_auto() ;
757 
758  const Vector& dlapconf_bh_auto_rs = hole.get_d_lapconf_auto_rs() ;
759  /*
760  Vector dlc_bh_auto_rs(mp_bh, COV, mp_bh.get_bvect_cart()) ;
761  dlc_bh_auto_rs.set_etat_qcq() ;
762  for (int i=1; i<=3; i++) {
763  dlc_bh_auto_rs.set(i) = lapconf_bh_auto_rs.deriv(i) ;
764  }
765  dlc_bh_auto_rs.std_spectral_base() ;
766  */
767 
768  const Vector& dconfo_bh_auto_rs = hole.get_d_confo_auto_rs() ;
769  /*
770  Vector dc_bh_auto_rs(mp_bh, COV, mp_bh.get_bvect_cart()) ;
771  dc_bh_auto_rs.set_etat_qcq() ;
772  for (int i=1; i<=3; i++) {
773  dc_bh_auto_rs.set(i) = confo_bh_auto_rs.deriv(i) ;
774  }
775  dc_bh_auto_rs.std_spectral_base() ;
776  */
777 
778  const Vector& dlapconf_ns_auto = star.get_d_lapconf_auto() ;
779  /*
780  Vector dlc_ns_auto(mp_ns, COV, mp_ns.get_bvect_cart()) ;
781  dlc_ns_auto.set_etat_qcq() ;
782  for (int i=1; i<=3; i++) {
783  dlc_ns_auto.set(i) = lapconf_ns_auto.deriv(i) ;
784  }
785  dlc_ns_auto.std_spectral_base() ;
786  */
787 
788  const Vector& dconfo_ns_auto = star.get_d_confo_auto() ;
789  /*
790  Vector dc_ns_auto(mp_ns, COV, mp_ns.get_bvect_cart()) ;
791  dc_ns_auto.set_etat_qcq() ;
792  for (int i=1; i<=3; i++) {
793  dc_ns_auto.set(i) = confo_ns_auto.deriv(i) ;
794  }
795  dc_ns_auto.std_spectral_base() ;
796  */
797  double mass = ggrav * hole.get_mass_bh() ;
798 
799  if (hole.is_kerrschild()) {
800 
801  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
802  abort() ;
803  /*
804  const Vector& shift_auto_rs = hole.get_shift_auto_rs() ;
805  const Vector& shift_comp = hole.get_shift_comp() ;
806 
807  Scalar llshift(mp_bh) ; // dzpuis = 0
808  llshift = ll(1) % (shift_auto_rs(1) + shift_comp(1))
809  + ll(2) % (shift_auto_rs(2) + shift_comp(2))
810  + ll(3) % (shift_auto_rs(3) + shift_comp(3)) ;
811  llshift.std_spectral_base() ;
812 
813  Scalar lldlldlap(mp_bh) ; // dzpuis = 3
814  lldlldlap = lldlap.dsdr() ;
815  lldlldlap.std_spectral_base() ;
816 
817  Scalar lap_bh2(mp_bh) ;
818  lap_bh2 = 1./(1.+2.*mass/rr) ;
819  lap_bh2.std_spectral_base() ;
820 
821  Scalar omelld(mp_bh) ;
822  omelld = omega * (ll(2) * (mp_bh.get_ori_x() - x_rot)
823  - ll(1) * (mp_bh.get_ori_y() - y_rot)) ;
824  omelld.std_spectral_base() ;
825 
826  //-------------------------------------
827  // Integration over the BH map
828  //-------------------------------------
829 
830  // Surface integral <- dzpuis should be 0
831  // --------------------------------------
832  source_bh_surf = lldlap ;
833 
834  source_bh_surf.std_spectral_base() ;
835  source_bh_surf.annule_domain(0) ;
836  source_bh_surf.raccord(1) ;
837  source_bh_surf.dec_dzpuis(2) ;
838 
839  integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
840 
841  // Volume integral <- dzpuis should be 4
842  // -------------------------------------
843  Scalar sou_bh1(mp_bh) ;
844  sou_bh1 = -2. * pow(lap_bh2,2.5) * mass * lldlconf / rr / rr ;
845  sou_bh1.std_spectral_base() ;
846  sou_bh1.inc_dzpuis(2) ;
847 
848  Scalar sou_bh2(mp_bh) ;
849  sou_bh2 = 4.*mass*mass*pow(lap_bh2,3.)*(1.+3.*mass/rr)
850  *(1.+3.*mass/rr)*(lapse_bh_auto_rs+lapse_bh_comp)
851  *pow(confo_bh,4.)/3./pow(rr,4.) ;
852  sou_bh2.std_spectral_base() ;
853  sou_bh2.inc_dzpuis(4) ;
854 
855  Scalar sou_bh3(mp_bh) ;
856  sou_bh3 = - 2.*mass*pow(lap_bh2,2.5)
857  *(2.+10.*mass/rr+9.*mass*mass/rr/rr)
858  *pow(confo_bh,4.)*(llshift+omelld)/pow(rr,3.) ;
859  sou_bh3.std_spectral_base() ;
860  sou_bh3.inc_dzpuis(4) ;
861 
862  Scalar sou_bh4(mp_bh) ;
863  sou_bh4 = 2. * lap_bh2 * mass * lldlldlap / rr ;
864  sou_bh4.std_spectral_base() ;
865  sou_bh4.inc_dzpuis(1) ;
866 
867  Scalar sou_bh5(mp_bh) ;
868  sou_bh5 = lap_bh2*lap_bh2*mass*(3.+8.*mass/rr)*lldlap/rr/rr ;
869  sou_bh5.std_spectral_base() ;
870  sou_bh5.inc_dzpuis(2) ;
871 
872  Scalar sou_bh6(mp_bh) ;
873  sou_bh6 = 4.*pow(lap_bh2,3.5)*mass*mass
874  *(2.*(sqrt(lap_bh2)/lapse_bh - 1.)*pow(confo_bh,4.)
875  *(4.+12.*mass/rr+9.*mass*mass/rr/rr) + 3.*(pow(confo_bh,4.)-1.))
876  /3./pow(rr,4.) ;
877  sou_bh6.std_spectral_base() ;
878  sou_bh6.inc_dzpuis(4) ;
879 
880  source_bh_volm = lapse_bh * (taij_quad_tot_rs + taij_quad_tot_rot)
881  / pow(confo_bh,8.)
882  - 2. * dlapdlcf + 4. * lap_bh2 * mass * lldlap * lldlconf / rr
883  + sou_bh1 + sou_bh2 + sou_bh3 + sou_bh4 + sou_bh5 + sou_bh6 ;
884 
885  source_bh_volm.std_spectral_base() ;
886  source_bh_volm.annule_domain(0) ;
887 
888  integ_bh_v = source_bh_volm.integrale() ;
889 
890  //-------------------------------------
891  // Integration over the NS map
892  //-------------------------------------
893 
894  // Volume integral <- dzpuis should be 4
895  // -------------------------------------
896  source_ns_volm = lapse_ns * psi4_ns * (ener_euler + s_euler) ;
897  source_ns_volm.std_spectral_base() ;
898  source_ns_volm.inc_dzpuis(4) ;
899 
900  integ_ns_v = source_ns_volm.integrale() ;
901 
902  cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
903  << " integ_bh_v : "
904  << integ_bh_v/ qpig / msol
905  << " integ_ns_v : " << integ_ns_v/ msol << endl ;
906 
907  //--------------------
908  // Komar mass
909  //--------------------
910  mkom = hole.get_mass_bh()
911  + (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
912 
913  cout << "Komar mass : " << mkom / msol << " [Mo]"
914  << endl ;
915 
916  // Komar mass by surface integral at infinity : dzpuis should be 2
917  // ------------------------------------------
918  double mmm = hole.get_mass_bh()
919  + (mp_aff.integrale_surface_infini(lldlap))/qpig ;
920 
921  cout << "Another Komar mass : " << mmm / msol << " [Mo]" << endl ;
922  */
923  }
924  else { // Isotropic coordinates with the maximal slicing
925 
926  //-------------------------------------
927  // Integration over the BH map
928  //-------------------------------------
929 
930  // Sets C/M^2 for each case of the lapse boundary condition
931  // --------------------------------------------------------
932  double cc ;
933 
934  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
935  if (hole.has_bc_lapconf_fs()) { // First condition
936  // d(\alpha \psi)/dr = 0
937  // ---------------------
938  cc = 2. * (sqrt(13.) - 1.) / 3. ;
939  }
940  else { // Second condition
941  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
942  // -----------------------------------------
943  cc = 4. / 3. ;
944  }
945  }
946  else { // Dirichlet boundary condition
947  if (hole.has_bc_lapconf_fs()) { // First condition
948  // (\alpha \psi) = 1/2
949  // -------------------
950  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
951  abort() ;
952  }
953  else { // Second condition
954  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
955  // ----------------------------------
956  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
957  abort() ;
958  // cc = 2. * sqrt(2.) ;
959  }
960  }
961 
962  Scalar r_are(mp_bh) ;
965  r_are.std_spectral_base() ;
966 
967  Scalar lldlapconf_is(mp_bh) ;
968  lldlapconf_is = ll(1)%dlapconf_bh_auto_rs(1)
969  + ll(2)%dlapconf_bh_auto_rs(2) + ll(3)%dlapconf_bh_auto_rs(3)
970  + ll(1)%dlapconf_bh_comp(1) + ll(2)%dlapconf_bh_comp(2)
971  + ll(3)%dlapconf_bh_comp(3) ;
972  // dzpuis = 2
973  lldlapconf_is.std_spectral_base() ;
974 
975  // Surface integral <- dzpuis should be 0
976  // --------------------------------------
977  Scalar divshift_zero(divshift) ;
978  divshift_zero.dec_dzpuis(2) ;
979 
980  Scalar lldllsh_zero(lldllsh) ;
981  lldllsh_zero.dec_dzpuis(2) ;
982 
983  Scalar sou_bh_s_psi(mp_bh) ;
984  sou_bh_s_psi = 0.5 * confo_bh / rr
985  - pow(confo_bh, 4.) * (divshift_zero - 3.*lldllsh_zero)
986  / 12. / lapconf_bh
987  - 0.5 * pow(confo_bh, 4.) * mass * mass * cc
988  * sqrt(1. -2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
989  / lapconf_bh / pow(r_are*rr,3.) ;
990 
991  sou_bh_s_psi.std_spectral_base() ;
992  sou_bh_s_psi.annule_domain(0) ;
993  sou_bh_s_psi.raccord(1) ;
994 
995  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
996  if (hole.has_bc_lapconf_fs()) { // First condition
997 
998  source_bh_surf = sou_bh_s_psi ;
999 
1000  source_bh_surf.std_spectral_base() ;
1001  source_bh_surf.annule_domain(0) ;
1002  source_bh_surf.raccord(1) ;
1003 
1004  }
1005  else {
1006 
1007  Scalar sou_bh_s1(mp_bh) ;
1008  sou_bh_s1 = 0.5 * lapconf_bh / rr ;
1009  sou_bh_s1.std_spectral_base() ;
1010  sou_bh_s1.annule_domain(0) ;
1011  sou_bh_s1.raccord(1) ;
1012 
1013  source_bh_surf = sou_bh_s1 + sou_bh_s_psi ;
1014 
1015  source_bh_surf.std_spectral_base() ;
1016  source_bh_surf.annule_domain(0) ;
1017  source_bh_surf.raccord(1) ;
1018 
1019  }
1020  }
1021  else {
1022 
1023  Scalar sou_bh_s1(mp_bh) ;
1024  sou_bh_s1 = lldlapconf_is ;
1025  sou_bh_s1.std_spectral_base() ;
1026  sou_bh_s1.dec_dzpuis(2) ;
1027 
1028  Scalar sou_bh_s2(mp_bh) ;
1029  sou_bh_s2 = 0.5 * sqrt(r_are)
1030  * (1. - 3.*cc*cc*pow(mass/r_are/rr,4.)
1031  -sqrt(1. - 2.*mass/r_are/rr
1032  + cc*cc*pow(mass/r_are/rr,4.))) / rr ;
1033 
1034  sou_bh_s2.std_spectral_base() ;
1035 
1036  source_bh_surf = sou_bh_s1 + sou_bh_s2 + sou_bh_s_psi ;
1037 
1038  source_bh_surf.std_spectral_base() ;
1039  source_bh_surf.annule_domain(0) ;
1040  source_bh_surf.raccord(1) ;
1041 
1042  }
1043 
1044  integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
1045 
1046  // Volume integral <- dzpuis should be 4
1047  // -------------------------------------
1048  Scalar sou_bh1(mp_bh) ;
1049  sou_bh1 = 0.75 * pow(mass*mass*cc,2.)
1050  * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
1051  * (7.*pow(confo_bh,6.)/lapconf_bh
1052  + pow(confo_bh,7.)/lapconf_bh/lapconf_bh)
1053  / pow(r_are*rr,6.) ;
1054 
1055  sou_bh1.std_spectral_base() ;
1056  sou_bh1.inc_dzpuis(4) ;
1057 
1058  Scalar sou_bh2(mp_bh) ;
1059  sou_bh2 = 0.125 * (7.*lapconf_bh/pow(confo_bh,8.)
1060  + 1./pow(confo_bh,7.)) * taij_quad_auto_bh ;
1061 
1062  sou_bh2.std_spectral_base() ;
1063 
1064  Scalar sou_bh3(mp_bh) ;
1065  sou_bh3 = 0.125 * (7.*((lapconf_bh_comp+0.5)/lapconf_bh
1066  * pow(confo_bh/(confo_bh_comp+0.5),6.) - 1.)
1067  * (lapconf_bh_comp+0.5)
1068  / pow(confo_bh_comp+0.5,8.)
1069  + (pow(confo_bh/(confo_bh_comp+0.5),7.)
1070  *pow((lapconf_bh_comp+0.5)/lapconf_bh,2.)
1071  - 1.) / pow(confo_bh_comp+0.5,7))
1072  * taij_quad_comp_bh ;
1073 
1074  sou_bh3.std_spectral_base() ;
1075 
1076  source_bh_volm = sou_bh1 + sou_bh2 + sou_bh3 ;
1077  source_bh_volm.std_spectral_base() ;
1078  source_bh_volm.annule_domain(0) ;
1079 
1080  integ_bh_v = source_bh_volm.integrale() ;
1081 
1082  //-------------------------------------
1083  // Integration over the NS map
1084  //-------------------------------------
1085 
1086  // Volume integral <- dzpuis should be 4
1087  // -------------------------------------
1088  Scalar sou_ns1(mp_ns) ;
1089  sou_ns1 = 0.5 * pow(confo_ns,4.) * (lapconf_ns + confo_ns)
1090  * ener_euler + lapconf_ns * pow(confo_ns,4.) * s_euler ;
1091  sou_ns1.std_spectral_base() ;
1092  sou_ns1.inc_dzpuis(4) ;
1093 
1094  Scalar sou_ns2(mp_ns) ;
1095  sou_ns2 = 0.125 * (7.*(lapconf_ns_auto+0.5)/(confo_ns_auto+0.5)
1096  + 1.) * taij_quad_auto_ns
1097  / pow(confo_ns_auto+0.5, 7.) / qpig ;
1098  sou_ns2.std_spectral_base() ;
1099 
1100  source_ns_volm = sou_ns1 + sou_ns2 ;
1101  source_ns_volm.std_spectral_base() ;
1102 
1103  integ_ns_v = source_ns_volm.integrale() ;
1104 
1105  cout << "integ_bh_s : " << integ_bh_s/ qpig / msol
1106  << " integ_bh_v : "
1107  << integ_bh_v/ qpig / msol
1108  << " integ_ns_v : " << integ_ns_v/ msol << endl ;
1109 
1110  //--------------------
1111  // Komar mass
1112  //--------------------
1113  mkom = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
1114 
1115  cout << "Komar mass (volume) : " << mkom / msol << " [Mo]"
1116  << endl ;
1117 
1118  }
1119 
1120  p_mass_kom_bhns_vol = new double( mkom ) ;
1121 
1122  }
1123 
1124  return *p_mass_kom_bhns_vol ;
1125 
1126 }
1127 
1128 
1129  //-----------------------------------------//
1130  // Total linear momentum //
1131  //-----------------------------------------//
1132 
1134 
1135  // Fundamental constants and units
1136  // -------------------------------
1137  using namespace Unites ;
1138 
1139  if (p_line_mom_bhns == 0x0) { // a new computation is required
1140 
1141  p_line_mom_bhns = new Tbl(3) ;
1142  p_line_mom_bhns->annule_hard() ; // fills the double array with zeros
1143 
1144  if (hole.is_kerrschild()) {
1145 
1146  // Construction of a new grid and a new affine mapping
1147  // ---------------------------------------------------
1148  int nz = (hole.get_mp()).get_mg()->get_nzone() ;
1149  double* bornes = new double[nz+1] ;
1150  double radius = separ + 2. ;
1151 
1152  for (int i=nz-1; i>0; i--) {
1153  bornes[i] = radius ;
1154  radius /= 2. ;
1155  }
1156  bornes[0] = 0. ;
1157  bornes[nz] = __infinity ;
1158 
1159  Map_af mp_aff(*((hole.get_mp()).get_mg()), bornes) ;
1160  mp_aff.set_ori(0.,0.,0.) ;
1161 
1162  delete [] bornes ;
1163 
1164  // Construction of the extrinsic curvature
1165  // ---------------------------------------
1166  Vector shift_bh(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1167  shift_bh.set_etat_qcq() ;
1168 
1169  Vector shift_ns(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1170  shift_ns.set_etat_qcq() ;
1171 
1172  shift_bh.set(1).import(hole.get_shift_auto()(1)) ;
1173  shift_bh.set(2).import(hole.get_shift_auto()(2)) ;
1174  shift_bh.set(3).import(hole.get_shift_auto()(3)) ;
1175 
1176  shift_ns.set(1).import(star.get_shift_auto()(1)) ;
1177  shift_ns.set(2).import(star.get_shift_auto()(2)) ;
1178  shift_ns.set(3).import(star.get_shift_auto()(3)) ;
1179 
1180  Vector shift_tot = shift_bh + shift_ns ;
1181  shift_tot.std_spectral_base() ;
1182  shift_tot.annule(0, nz-2) ;
1183 
1184  Scalar stc(mp_aff) ;
1185  stc = mp_aff.sint ;
1186  stc.std_spectral_base() ;
1187  Scalar ctc(mp_aff) ;
1188  ctc = mp_aff.cost ;
1189  ctc.std_spectral_base() ;
1190  Scalar spc(mp_aff) ;
1191  spc = mp_aff.sinp ;
1192  spc.std_spectral_base() ;
1193  Scalar cpc(mp_aff) ;
1194  cpc = mp_aff.cosp ;
1195  cpc.std_spectral_base() ;
1196 
1197  Vector lc(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1198  lc.set_etat_qcq() ;
1199  lc.set(1) = stc * cpc ;
1200  lc.set(2) = stc * spc ;
1201  lc.set(3) = ctc ;
1202  lc.std_spectral_base() ;
1203 
1204  Vector kijlj(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1205  kijlj.set_etat_qcq() ;
1206 
1207  Scalar rr(mp_aff) ;
1208  rr = mp_aff.r ;
1209  rr.std_spectral_base() ;
1210 
1211  Scalar xbhsr(mp_aff) ;
1212  xbhsr = (hole.get_mp()).get_ori_x() / rr ;
1213  xbhsr.std_spectral_base() ;
1214 
1215  Scalar ybhsr(mp_aff) ;
1216  ybhsr = (hole.get_mp()).get_ori_y() / rr ;
1217  ybhsr.std_spectral_base() ;
1218 
1219  Scalar rl(mp_aff) ;
1220  rl = sqrt(1. - 2.*xbhsr*lc(1) - 2.*ybhsr*lc(2)
1221  + xbhsr*xbhsr + ybhsr*ybhsr) ;
1222  rl.std_spectral_base() ;
1223 
1224  Vector ll(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1225  ll.set_etat_qcq() ;
1226  ll.set(1) = (lc(1) - xbhsr) / rl ;
1227  ll.set(2) = (lc(2) - ybhsr)/ rl ;
1228  ll.set(3) = lc(3) / rl ;
1229  ll.std_spectral_base() ;
1230 
1231  Scalar lcll(mp_aff) ;
1232  lcll = lc(1)*ll(1) + lc(2)*ll(2) + lc(3)*ll(3) ;
1233  lcll.std_spectral_base() ;
1234 
1235  Scalar divshift(mp_aff) ;
1236  divshift = shift_tot(1).deriv(1) + shift_tot(2).deriv(2)
1237  + shift_tot(3).deriv(3) ;
1238  divshift.std_spectral_base() ;
1239 
1240  Scalar llshift(mp_aff) ;
1241  llshift = ll(1)*shift_tot(1) + ll(2)*shift_tot(2)
1242  + ll(3)*shift_tot(3) ;
1243  llshift.std_spectral_base() ;
1244 
1245  Scalar lcshift(mp_aff) ;
1246  lcshift = lc(1)*shift_tot(1) + lc(2)*shift_tot(2)
1247  + lc(3)*shift_tot(3) ;
1248  lcshift.std_spectral_base() ;
1249 
1250  Vector lcdshft(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1251  for (int i=1; i<=3; i++) {
1252  lcdshft.set(i) = lc(1)*(shift_tot(1).deriv(i))
1253  + lc(2)*(shift_tot(2).deriv(i))
1254  + lc(3)*(shift_tot(3).deriv(i)) ;
1255  }
1256  lcdshft.std_spectral_base() ;
1257 
1258  Vector dshift(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1259  for (int i=1; i<=3; i++) {
1260  dshift.set(i) = shift_tot(i).dsdr() ;
1261  }
1262  dshift.std_spectral_base() ;
1263 
1264  Vector lldshft(mp_aff, CON, mp_aff.get_bvect_cart()) ;
1265  for (int i=1; i<=3; i++) {
1266  lldshft.set(i) = ll(1)*(shift_tot(i).deriv(1))
1267  + ll(2)*(shift_tot(i).deriv(2))
1268  + ll(3)*(shift_tot(i).deriv(3)) ;
1269  }
1270  lldshft.std_spectral_base() ;
1271 
1272  double mass = ggrav * hole.get_mass_bh() ;
1273 
1274  Scalar lap_bh2(mp_aff) ;
1275  lap_bh2 = 1./(1.+2.*mass/rl/rr) ;
1276  lap_bh2.std_spectral_base() ;
1277 
1278  Scalar lap2hbh(mp_aff) ;
1279  lap2hbh = lap_bh2 * mass / rl / rr ;
1280  lap2hbh.std_spectral_base() ;
1281 
1282  Scalar omexsr(mp_aff) ;
1283  omexsr = omega * ((hole.get_mp()).get_ori_x() - x_rot)
1284  / rl / rr ;
1285  omexsr.std_spectral_base() ;
1286 
1287  Scalar omeysr(mp_aff) ;
1288  omeysr = omega * ((hole.get_mp()).get_ori_y() - y_rot)
1289  / rl / rr ;
1290  omeysr.std_spectral_base() ;
1291 
1292  Scalar kk(mp_aff) ;
1293  kk = 4.*sqrt(lap_bh2)*lap2hbh*(1.+3.*mass/rl/rr)/3./rl/rr ;
1294  kk.std_spectral_base() ;
1295 
1296  //-----------------------------------------------------------
1297  // Surface integral at infinity : dzpuis should be 2
1298  //-----------------------------------------------------------
1299 
1300  // Source for X component
1301  // ----------------------
1302  Scalar kij_x1(mp_aff) ;
1303  kij_x1 = omexsr*ll(1)*lc(2) - omeysr*(ll(1)*lc(1)+lcll)
1304  + lcll*shift_tot(1)/rl/rr
1305  + ll(1)*lcshift/rl/rr
1306  + (lc(1)-lap_bh2*(9.+14.*mass/rl/rr)*ll(1)*lcll)
1307  *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
1308  kij_x1.std_spectral_base() ;
1309  kij_x1.inc_dzpuis(2) ;
1310 
1311  Scalar kij_x2(mp_aff) ;
1312  kij_x2 = kk * (lc(1) - 2.*lap2hbh*ll(1)*lcll) ;
1313  kij_x2.std_spectral_base() ;
1314  kij_x2.inc_dzpuis(2) ;
1315 
1316  kijlj.set(1) = lcdshft(1) + dshift(1) - 2.*lc(1)*divshift/3.
1317  + 2.*lap2hbh * (-ll(1)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
1318  + ll(3)*lcdshft(3))
1319  - lcll*lldshft(1)
1320  + 2.*ll(1)*lcll*divshift/3.
1321  + kij_x1)
1322  + kij_x2 ;
1323 
1324  // Source for Y component
1325  // ----------------------
1326  Scalar kij_y1(mp_aff) ;
1327  kij_y1 = omexsr*(ll(2)*lc(2)+lcll) - omeysr*ll(2)*lc(1)
1328  + lcll*shift_tot(2)/rl/rr
1329  + ll(2)*lcshift/rl/rr
1330  + (lc(2)-lap_bh2*(9.+14.*mass/rl/rr)*ll(2)*lcll)
1331  *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
1332  kij_y1.std_spectral_base() ;
1333  kij_y1.inc_dzpuis(2) ;
1334 
1335  Scalar kij_y2(mp_aff) ;
1336  kij_y2 = kk * (lc(2) - 2.*lap2hbh*ll(2)*lcll) ;
1337  kij_y2.std_spectral_base() ;
1338  kij_y2.inc_dzpuis(2) ;
1339 
1340  kijlj.set(2) = lcdshft(2) + dshift(2) - 2.*lc(2)*divshift/3.
1341  + 2.*lap2hbh * (-ll(2)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
1342  + ll(3)*lcdshft(3))
1343  - lcll*lldshft(2)
1344  + 2.*ll(2)*lcll*divshift/3.
1345  + kij_y1)
1346  + kij_y2 ;
1347 
1348  // Source for Z component
1349  // ----------------------
1350  Scalar kij_z1(mp_aff) ;
1351  kij_z1 = omexsr*ll(3)*lc(2) - omeysr*ll(3)*lc(1)
1352  + lcll*shift_tot(3)/rl/rr
1353  + ll(3)*lcshift/rl/rr
1354  + (lc(3)-lap_bh2*(9.+14.*mass/rl/rr)*ll(3)*lcll)
1355  *(llshift/rl/rr + omexsr*ll(2) - omeysr*ll(1))/3. ;
1356  kij_z1.std_spectral_base() ;
1357  kij_z1.inc_dzpuis(2) ;
1358 
1359  Scalar kij_z2(mp_aff) ;
1360  kij_z2 = kk * (lc(3) - 2.*lap2hbh*ll(3)*lcll) ;
1361  kij_z2.std_spectral_base() ;
1362  kij_z2.inc_dzpuis(2) ;
1363 
1364  kijlj.set(3) = lcdshft(3) + dshift(3) - 2.*lc(3)*divshift/3.
1365  + 2.*lap2hbh * (-ll(3)*(ll(1)*lcdshft(1) + ll(2)*lcdshft(2)
1366  + ll(3)*lcdshft(3))
1367  - lcll*lldshft(3)
1368  + 2.*ll(3)*lcll*divshift/3.
1369  + kij_z1)
1370  + kij_z2 ;
1371 
1372  kijlj.std_spectral_base() ;
1373 
1374  // X component dzpuis should be 2
1375  // -----------
1376  double lm_x = mp_aff.integrale_surface_infini(kijlj(1)) ;
1377  p_line_mom_bhns->set(0) = 0.25 * lm_x / qpig ;
1378 
1379  // Y component
1380  // -----------
1381  double lm_y = mp_aff.integrale_surface_infini(kijlj(2)) ;
1382  p_line_mom_bhns->set(1) = 0.25 * lm_y / qpig ;
1383 
1384  // Z component
1385  // -----------
1386  double lm_z = mp_aff.integrale_surface_infini(kijlj(3)) ;
1387  p_line_mom_bhns->set(2) = 0.25 * lm_z / qpig ;
1388 
1389  }
1390  else { // Isotropic coordinates with the maximal slicing
1391 
1392  /*
1393  // Method by the sourface integral at infinity
1394  // -------------------------------------------
1395 
1396  const Map& mp_bh = hole.get_mp() ;
1397  Map_af mp_aff(mp_bh) ;
1398 
1399  Scalar st(mp_bh) ;
1400  st = mp_bh.sint ;
1401  st.std_spectral_base() ;
1402  Scalar ct(mp_bh) ;
1403  ct = mp_bh.cost ;
1404  ct.std_spectral_base() ;
1405  Scalar sp(mp_bh) ;
1406  sp = mp_bh.sinp ;
1407  sp.std_spectral_base() ;
1408  Scalar cp(mp_bh) ;
1409  cp = mp_bh.cosp ;
1410  cp.std_spectral_base() ;
1411 
1412  Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1413  ll.set_etat_qcq() ;
1414  ll.set(1) = st * cp ;
1415  ll.set(2) = st * sp ;
1416  ll.set(3) = ct ;
1417  ll.std_spectral_base() ;
1418 
1419  const Scalar& confo_bh = hole.get_confo_tot() ;
1420  const Sym_tensor& taij_tot_rs = hole.get_taij_tot_rs() ;
1421 
1422  Vector kijlj(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1423  kijlj.set_etat_qcq() ;
1424 
1425  for (int i=1; i<=3; i++) {
1426  kijlj.set(i) = (taij_tot_rs(i,1)%ll(1)
1427  + taij_tot_rs(i,2)%ll(2)
1428  + taij_tot_rs(i,3)%ll(3)) / pow(confo_bh,10.) ;
1429  }
1430 
1431  kijlj.std_spectral_base() ;
1432 
1433  // X component
1434  // -----------
1435  double lm_x = mp_aff.integrale_surface_infini(kijlj(1)) ;
1436  p_line_mom_bhns->set(0) = 0.5 * lm_x / qpig ;
1437 
1438  // Y component
1439  // -----------
1440  double lm_y = mp_aff.integrale_surface_infini(kijlj(2)) ;
1441  p_line_mom_bhns->set(1) = 0.5 * lm_y / qpig ;
1442 
1443  // Z component
1444  // -----------
1445  double lm_z = mp_aff.integrale_surface_infini(kijlj(3)) ;
1446  p_line_mom_bhns->set(2) = 0.5 * lm_z / qpig ;
1447  */
1448 
1449  // Method by the volume integral and the surface integral
1450  // at the BH horizon
1451  // ------------------------------------------------------
1452 
1453  const Map& mp_bh = hole.get_mp() ;
1454  Map_af mp_aff(mp_bh) ;
1455  const Map& mp_ns = star.get_mp() ;
1456 
1457  const Sym_tensor& taij = hole.get_taij_tot() ;
1458  const Scalar& confo_ns = star.get_confo_tot() ;
1459  const Scalar& ee = star.get_ener_euler() ;
1460  const Scalar& pp = star.get_press() ;
1461  const Vector& uu = star.get_u_euler() ;
1462 
1463  double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
1464 
1465  Scalar st(mp_bh) ;
1466  st = mp_bh.sint ;
1467  st.std_spectral_base() ;
1468  Scalar ct(mp_bh) ;
1469  ct = mp_bh.cost ;
1470  ct.std_spectral_base() ;
1471  Scalar sp(mp_bh) ;
1472  sp = mp_bh.sinp ;
1473  sp.std_spectral_base() ;
1474  Scalar cp(mp_bh) ;
1475  cp = mp_bh.cosp ;
1476  cp.std_spectral_base() ;
1477 
1478  Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1479  ll.set_etat_qcq() ;
1480  ll.set(1) = st % cp ;
1481  ll.set(2) = st % sp ;
1482  ll.set(3) = ct ;
1483  ll.std_spectral_base() ;
1484 
1485  // X component
1486  // -----------
1487 
1488  //-------------------------------------
1489  // Integration over the BH map
1490  //-------------------------------------
1491 
1492  // Surface integral <- dzpuis should be 0
1493  // --------------------------------------
1494  Scalar sou_bh_sx(mp_bh) ;
1495  sou_bh_sx = taij(1,1) * ll(1) + taij(1,2) * ll(2)
1496  + taij(1,3) * ll(3) ;
1497  sou_bh_sx.std_spectral_base() ;
1498  sou_bh_sx.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1499 
1500  sou_bh_sx.annule_domain(0) ;
1501  sou_bh_sx.raccord(1) ;
1502 
1503  double integ_bh_s_x = mp_aff.integrale_surface(sou_bh_sx,
1504  radius_ah) ;
1505 
1506  //-------------------------------------
1507  // Integration over the NS map
1508  //-------------------------------------
1509 
1510  // Volume integral <- dzpuis should be 4
1511  // -------------------------------------
1512  Scalar sou_ns_vx(mp_ns) ;
1513 
1514  sou_ns_vx = pow(confo_ns, 10.) * (ee + pp) * uu(1) ;
1515 
1516  sou_ns_vx.std_spectral_base() ;
1517  sou_ns_vx.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1518 
1519  double integ_ns_v_x = sou_ns_vx.integrale() ;
1520 
1521  p_line_mom_bhns->set(0) = integ_ns_v_x + 0.5*integ_bh_s_x/qpig ;
1522 
1523  // Y component
1524  // -----------
1525 
1526  //-------------------------------------
1527  // Integration over the BH map
1528  //-------------------------------------
1529 
1530  // Surface integral <- dzpuis should be 0
1531  // --------------------------------------
1532  Scalar sou_bh_sy(mp_bh) ;
1533  sou_bh_sy = taij(2,1) * ll(1) + taij(2,2) * ll(2)
1534  + taij(2,3) * ll(3) ;
1535  sou_bh_sy.std_spectral_base() ;
1536  sou_bh_sy.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1537 
1538  sou_bh_sy.annule_domain(0) ;
1539  sou_bh_sy.raccord(1) ;
1540 
1541  double integ_bh_s_y = mp_aff.integrale_surface(sou_bh_sy,
1542  radius_ah) ;
1543 
1544  //-------------------------------------
1545  // Integration over the NS map
1546  //-------------------------------------
1547 
1548  // Volume integral <- dzpuis should be 4
1549  // -------------------------------------
1550  Scalar sou_ns_vy(mp_ns) ;
1551 
1552  sou_ns_vy = pow(confo_ns, 10.) * (ee + pp) * uu(2) ;
1553 
1554  sou_ns_vy.std_spectral_base() ;
1555  sou_ns_vy.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1556 
1557  double integ_ns_v_y = sou_ns_vy.integrale() ;
1558 
1559  p_line_mom_bhns->set(1) = integ_ns_v_y + 0.5*integ_bh_s_y/qpig ;
1560 
1561 
1562  // Z component
1563  // -----------
1564 
1565  //-------------------------------------
1566  // Integration over the BH map
1567  //-------------------------------------
1568 
1569  // Surface integral <- dzpuis should be 0
1570  // --------------------------------------
1571  Scalar sou_bh_sz(mp_bh) ;
1572  sou_bh_sz = taij(3,1) * ll(1) + taij(3,2) * ll(2)
1573  + taij(3,3) * ll(3) ;
1574  sou_bh_sz.std_spectral_base() ;
1575  sou_bh_sz.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1576 
1577  sou_bh_sz.annule_domain(0) ;
1578  sou_bh_sz.raccord(1) ;
1579 
1580  double integ_bh_s_z = mp_aff.integrale_surface(sou_bh_sz,
1581  radius_ah) ;
1582 
1583  //-------------------------------------
1584  // Integration over the NS map
1585  //-------------------------------------
1586 
1587  // Volume integral <- dzpuis should be 4
1588  // -------------------------------------
1589  Scalar sou_ns_vz(mp_ns) ;
1590 
1591  sou_ns_vz = pow(confo_ns, 10.) * (ee + pp) * uu(3) ;
1592 
1593  sou_ns_vz.std_spectral_base() ;
1594  sou_ns_vz.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1595 
1596  double integ_ns_v_z = sou_ns_vz.integrale() ;
1597 
1598  p_line_mom_bhns->set(2) = integ_ns_v_z + 0.5*integ_bh_s_z/qpig ;
1599 
1600  }
1601 
1602  }
1603 
1604  return *p_line_mom_bhns ;
1605 
1606 }
1607 
1608 
1609  //------------------------------------------//
1610  // Total angular momentum //
1611  //------------------------------------------//
1612 
1614 
1615  // Fundamental constants and units
1616  // -------------------------------
1617  using namespace Unites ;
1618 
1619  if (p_angu_mom_bhns == 0x0) { // a new computation is required
1620 
1621  p_angu_mom_bhns = new Tbl(3) ;
1622  p_angu_mom_bhns->annule_hard() ; // fills the double array with zeros
1623 
1624  double integ_bh_s_x ;
1625  double integ_bh_s_y ;
1626  double integ_bh_s_z ;
1627  double integ_bh_v_x ;
1628  double integ_bh_v_y ;
1629  double integ_bh_v_z ;
1630  double integ_ns_v_x ;
1631  double integ_ns_v_y ;
1632  double integ_ns_v_z ;
1633 
1634  const Map& mp_bh = hole.get_mp() ;
1635  const Map& mp_ns = star.get_mp() ;
1636 
1637  double radius_ah = mp_bh.val_r(1,-1.,M_PI/2.,0.) ;
1638  Map_af mp_aff(mp_bh) ;
1639 
1640  Scalar source_bh_surf(mp_bh) ;
1641  source_bh_surf.set_etat_qcq() ;
1642 
1643  Scalar source_bh_volm(mp_bh) ;
1644  source_bh_volm.set_etat_qcq() ;
1645 
1646  Scalar source_ns_volm(mp_ns) ;
1647  source_ns_volm.set_etat_qcq() ;
1648 
1649  Scalar rr(mp_bh) ;
1650  rr = mp_bh.r ;
1651  rr.std_spectral_base() ;
1652 
1653  Scalar st(mp_bh) ;
1654  st = mp_bh.sint ;
1655  st.std_spectral_base() ;
1656  Scalar ct(mp_bh) ;
1657  ct = mp_bh.cost ;
1658  ct.std_spectral_base() ;
1659  Scalar sp(mp_bh) ;
1660  sp = mp_bh.sinp ;
1661  sp.std_spectral_base() ;
1662  Scalar cp(mp_bh) ;
1663  cp = mp_bh.cosp ;
1664  cp.std_spectral_base() ;
1665 
1666  Vector ll(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1667  ll.set_etat_qcq() ;
1668  ll.set(1) = st % cp ;
1669  ll.set(2) = st % sp ;
1670  ll.set(3) = ct ;
1671  ll.std_spectral_base() ;
1672 
1673  Scalar xx_bh(mp_bh) ;
1674  xx_bh = mp_bh.xa ;
1675  xx_bh.std_spectral_base() ;
1676  Scalar yy_bh(mp_bh) ;
1677  yy_bh = mp_bh.ya ;
1678  yy_bh.std_spectral_base() ;
1679  Scalar zz_bh(mp_bh) ;
1680  zz_bh = mp_bh.za ;
1681  zz_bh.std_spectral_base() ;
1682 
1683  Scalar xx_ns(mp_ns) ;
1684  xx_ns = mp_ns.xa ;
1685  xx_ns.std_spectral_base() ;
1686  Scalar yy_ns(mp_ns) ;
1687  yy_ns = mp_ns.ya ;
1688  yy_ns.std_spectral_base() ;
1689  Scalar zz_ns(mp_ns) ;
1690  zz_ns = mp_ns.za ;
1691  zz_ns.std_spectral_base() ;
1692 
1693  const Scalar& confo_bh = hole.get_confo_tot() ;
1694  const Sym_tensor& taij = hole.get_taij_tot() ;
1695  const Scalar& confo_ns = star.get_confo_tot() ;
1696  const Scalar& ee = star.get_ener_euler() ;
1697  const Scalar& pp = star.get_press() ;
1698  const Vector& uu = star.get_u_euler() ;
1699 
1700  Vector jbh_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1701  jbh_x.set_etat_qcq() ;
1702  for (int i=1; i<=3; i++)
1703  jbh_x.set(i) = yy_bh * taij(3,i) - zz_bh * taij(2,i) ;
1704 
1705  jbh_x.std_spectral_base() ;
1706 
1707  Vector jbh_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1708  jbh_y.set_etat_qcq() ;
1709  for (int i=1; i<=3; i++)
1710  jbh_y.set(i) = zz_bh * taij(1,i) - xx_bh * taij(3,i) ;
1711 
1712  jbh_y.std_spectral_base() ;
1713 
1714  Vector jbh_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1715  jbh_z.set_etat_qcq() ;
1716  for (int i=1; i<=3; i++)
1717  jbh_z.set(i) = xx_bh * taij(2,i) - yy_bh * taij(1,i) ;
1718 
1719  jbh_z.std_spectral_base() ;
1720 
1721  double mass = ggrav * hole.get_mass_bh() ;
1722 
1723  if (hole.is_kerrschild()) {
1724 
1725  double ori_bh = mp_bh.get_ori_x() ;
1726 
1727  Scalar lap_bh2(mp_bh) ;
1728  lap_bh2 = 1./(1.+2.*mass/rr) ;
1729  lap_bh2.std_spectral_base() ;
1730 
1731  Scalar lcnf(mp_bh) ;
1732  lcnf = log(confo_bh) ;
1733  lcnf.std_spectral_base() ;
1734 
1735  Vector jbhsr_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1736  jbhsr_x.set_etat_qcq() ;
1737  for (int i=1; i<=3; i++)
1738  jbhsr_x.set(i) = ll(2)*taij(3,i) - ll(3)*taij(2,i) ;
1739 
1740  jbhsr_x.std_spectral_base() ;
1741 
1742  Vector jbhsr_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1743  jbhsr_y.set_etat_qcq() ;
1744  for (int i=1; i<=3; i++)
1745  jbhsr_y.set(i) = ll(3)*taij(1,i) - (ll(1)+ori_bh/rr)*taij(3,i) ;
1746 
1747  jbhsr_y.std_spectral_base() ;
1748 
1749  Vector jbhsr_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1750  jbhsr_z.set_etat_qcq() ;
1751  for (int i=1; i<=3; i++)
1752  jbhsr_z.set(i) = (ll(1)+ori_bh/rr)*taij(2,i) - ll(2)*taij(1,i) ;
1753 
1754  jbhsr_z.std_spectral_base() ;
1755 
1756  Scalar tmp1(mp_bh) ; // dzpuis = 0
1757  tmp1 = 2. * pow(lap_bh2,2.5) * mass * (1.+3.*mass/rr)
1758  * pow(confo_bh,6.) * ori_bh / 3. / rr / rr ;
1759  tmp1.std_spectral_base() ;
1760  tmp1.annule_domain(0) ;
1761 
1762  Scalar tmp2(mp_bh) ; // dzpuis = 0
1763  tmp2 = 4. * sqrt(lap_bh2) * (1.+3.*mass/rr) * pow(confo_bh,6.) ;
1764  tmp2.std_spectral_base() ;
1765  tmp2.annule_domain(0) ;
1766 
1767  Scalar lltaij(mp_bh) ; // dzpuis = 2
1768  lltaij = ll(1)*(ll(1)*taij(1,1)+ll(2)*taij(1,2)+ll(3)*taij(1,3))
1769  + ll(2)*(ll(1)*taij(2,1)+ll(2)*taij(2,2)+ll(3)*taij(2,3))
1770  + ll(3)*(ll(1)*taij(3,1)+ll(2)*taij(3,2)+ll(3)*taij(3,3)) ;
1771  lltaij.std_spectral_base() ;
1772  lltaij.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1773 
1774  Scalar dlcnf(mp_bh) ; // dzpuis = 2
1775  dlcnf = - 2. * lap_bh2 * tmp2 * mass * lcnf.dsdr() / rr ;
1776  dlcnf.std_spectral_base() ;
1777  dlcnf.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1778  dlcnf.annule_domain(0) ;
1779 
1780  Scalar tmp3(mp_bh) ; // dzpuis = 0
1781  tmp3 = -2.*pow(lap_bh2,2.5)
1782  *(6.+32.*mass/rr+41.*mass*mass/rr/rr+24.*pow(mass,3.)/pow(rr,3.))
1783  *pow(confo_bh,6.)/3./rr
1784  + 3.* lltaij + dlcnf ;
1785  tmp3.std_spectral_base() ;
1786  tmp3.annule_domain(0) ;
1787 
1788  Scalar tmp4(mp_bh) ; // dzpuis = 0
1789  tmp4 = lap_bh2 * mass / rr ;
1790  tmp4.std_spectral_base() ;
1791 
1792  //-------------//
1793  // X component //
1794  //-------------//
1795 
1796  //-------------------------------------
1797  // Integration over the BH map
1798  //-------------------------------------
1799 
1800  // Surface integral <- dzpuis should be 0
1801  // --------------------------------------
1802  source_bh_surf = jbh_x(1)*ll(1) + jbh_x(2)*ll(2) + jbh_x(3)*ll(3) ;
1803 
1804  source_bh_surf.std_spectral_base() ;
1805  source_bh_surf.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1806  source_bh_surf.annule_domain(0) ;
1807  source_bh_surf.raccord(1) ;
1808 
1809  integ_bh_s_x = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
1810 
1811  // Volume integral <- dzpuis should be 4
1812  // -------------------------------------
1813  source_bh_volm = tmp4
1814  * ( jbhsr_x(1)*ll(1) + jbhsr_x(2)*ll(2) + jbhsr_x(3)*ll(3)
1815  + tmp2 * ( ll(2)*lcnf.dsdz() - ll(3)*lcnf.dsdy() ) ) ;
1816 
1817  source_bh_volm.std_spectral_base() ;
1818  source_bh_volm.inc_dzpuis(2) ; // dzpuis : 2 -> 4
1819  source_bh_volm.annule_domain(0) ;
1820 
1821  integ_bh_v_x = source_bh_volm.integrale() ;
1822 
1823  //-------------------------------------
1824  // Integration over the NS map
1825  //-------------------------------------
1826 
1827  // Volume integral <- dzpuis should be 4
1828  // -------------------------------------
1829  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
1830  * (yy_ns*uu(3) - zz_ns*uu(2)) ;
1831 
1832  source_ns_volm.std_spectral_base() ;
1833  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1834 
1835  integ_ns_v_x = source_ns_volm.integrale() ;
1836 
1837  p_angu_mom_bhns->set(0) = integ_ns_v_x
1838  + 0.5 * (integ_bh_s_x + integ_bh_v_x) / qpig ;
1839 
1840  //-------------//
1841  // Y component //
1842  //-------------//
1843 
1844  //-------------------------------------
1845  // Integration over the BH map
1846  //-------------------------------------
1847 
1848  // Surface integral <- dzpuis should be 0
1849  // --------------------------------------
1850  Scalar jbh_y_ll(mp_bh) ;
1851  jbh_y_ll = jbh_y(1)*ll(1) + jbh_y(2)*ll(2) + jbh_y(3)*ll(3) ;
1852  jbh_y_ll.std_spectral_base() ;
1853  jbh_y_ll.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1854 
1855  source_bh_surf = jbh_y_ll - tmp1 * ll(3) ;
1856  source_bh_surf.std_spectral_base() ;
1857  source_bh_surf.annule_domain(0) ;
1858  source_bh_surf.raccord(1) ;
1859 
1860  integ_bh_s_y = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
1861 
1862  // Volume integral <- dzpuis should be 4
1863  // -------------------------------------
1864  Scalar tmp3_llz(mp_bh) ;
1865  tmp3_llz = tmp3 * ll(3) * ori_bh / rr ;
1866  tmp3_llz.std_spectral_base() ;
1867  tmp3_llz.inc_dzpuis(2) ; // dzpuis : 0 -> 2
1868 
1869  source_bh_volm = tmp4
1870  * ( jbhsr_y(1)*ll(1) + jbhsr_y(2)*ll(2) + jbhsr_y(3)*ll(3)
1871  + tmp2 * ( ll(3)*lcnf.dsdx() - (ll(1)+ori_bh/rr)*lcnf.dsdz() )
1872  - tmp3_llz ) ;
1873 
1874  source_bh_volm.std_spectral_base() ;
1875  source_bh_volm.inc_dzpuis(2) ; // dzpuis : 2 -> 4
1876  source_bh_volm.annule_domain(0) ;
1877 
1878  integ_bh_v_y = source_bh_volm.integrale() ;
1879 
1880  //-------------------------------------
1881  // Integration over the NS map
1882  //-------------------------------------
1883 
1884  // Volume integral <- dzpuis should be 4
1885  // -------------------------------------
1886  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
1887  * (zz_ns*uu(1) - xx_ns*uu(3)) ;
1888 
1889  source_ns_volm.std_spectral_base() ;
1890  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1891 
1892  integ_ns_v_y = source_ns_volm.integrale() ;
1893 
1894  p_angu_mom_bhns->set(1) = integ_ns_v_y
1895  + 0.5 * (integ_bh_s_y + integ_bh_v_y) / qpig ;
1896 
1897  //-------------//
1898  // Z component //
1899  //-------------//
1900 
1901  //-------------------------------------
1902  // Integration over the BH map
1903  //-------------------------------------
1904 
1905  // Surface integral <- dzpuis should be 0
1906  // --------------------------------------
1907  Scalar jbh_z_ll(mp_bh) ;
1908  jbh_z_ll = jbh_z(1)*ll(1) + jbh_z(2)*ll(2) + jbh_z(3)*ll(3) ;
1909  jbh_z_ll.std_spectral_base() ;
1910  jbh_z_ll.dec_dzpuis(2) ; // dzpuis : 2 -> 0
1911  source_bh_surf = jbh_z_ll + tmp1 * ll(2) ;
1912 
1913  source_bh_surf.std_spectral_base() ;
1914  source_bh_surf.annule_domain(0) ;
1915  source_bh_surf.raccord(1) ;
1916 
1917  integ_bh_s_z = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
1918 
1919  // Volume integral <- dzpuis should be 4
1920  // -------------------------------------
1921  Scalar tmp3_lly(mp_bh) ;
1922  tmp3_lly = tmp3 * ll(2) * ori_bh / rr ;
1923  tmp3_lly.std_spectral_base() ;
1924  tmp3_lly.inc_dzpuis(2) ; // dzpuis : 0 -> 2
1925 
1926  source_bh_volm = tmp4
1927  * ( jbhsr_z(1)*ll(1) + jbhsr_z(2)*ll(2) + jbhsr_z(3)*ll(3)
1928  + tmp2 * ( (ll(1)+ori_bh/rr)*lcnf.dsdy() - ll(2)*lcnf.dsdx() )
1929  + tmp3_lly ) ;
1930 
1931  source_bh_volm.std_spectral_base() ;
1932  source_bh_volm.inc_dzpuis(2) ; // dzpuis : 2 -> 4
1933  source_bh_volm.annule_domain(0) ;
1934 
1935  integ_bh_v_z = source_bh_volm.integrale() ;
1936 
1937  //-------------------------------------
1938  // Integration over the NS map
1939  //-------------------------------------
1940 
1941  // Volume integral <- dzpuis should be 4
1942  // -------------------------------------
1943  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
1944  * (xx_ns*uu(2) - yy_ns*uu(1)) ;
1945 
1946  source_ns_volm.std_spectral_base() ;
1947  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
1948 
1949  integ_ns_v_z = source_ns_volm.integrale() ;
1950 
1951  p_angu_mom_bhns->set(2) = integ_ns_v_z
1952  + 0.5 * (integ_bh_s_z + integ_bh_v_z) / qpig ;
1953 
1954  }
1955  else { // Isotropic coordinates with the maximal slicing
1956 
1957  // Sets C/M^2 for each case of the lapse boundary condition
1958  // --------------------------------------------------------
1959  double cc ;
1960 
1961  if (hole.has_bc_lapconf_nd()) { // Neumann boundary condition
1962  if (hole.has_bc_lapconf_fs()) { // First condition
1963  // d(\alpha \psi)/dr = 0
1964  // ---------------------
1965  cc = 2. * (sqrt(13.) - 1.) / 3. ;
1966  }
1967  else { // Second condition
1968  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
1969  // -----------------------------------------
1970  cc = 4. / 3. ;
1971  }
1972  }
1973  else { // Dirichlet boundary condition
1974  if (hole.has_bc_lapconf_fs()) { // First condition
1975  // (\alpha \psi) = 1/2
1976  // -------------------
1977  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
1978  abort() ;
1979  }
1980  else { // Second condition
1981  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
1982  // ----------------------------------
1983  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
1984  abort() ;
1985  // cc = 2. * sqrt(2.) ;
1986  }
1987  }
1988 
1989  Scalar r_are(mp_bh) ;
1990  r_are = hole.r_coord(hole.has_bc_lapconf_nd(),
1991  hole.has_bc_lapconf_fs()) ;
1992  r_are.std_spectral_base() ;
1993 
1994  const Scalar& lapconf_bh = hole.get_lapconf_tot() ;
1995  const Scalar& confo_bh = hole.get_confo_tot() ;
1996  const Sym_tensor& taij_rs = hole.get_taij_tot_rs() ;
1997 
1998  Vector jbh_rs_x(mp_bh, CON, mp_bh.get_bvect_cart()) ;
1999  jbh_rs_x.set_etat_qcq() ;
2000  for (int i=1; i<=3; i++)
2001  jbh_rs_x.set(i) = yy_bh * taij_rs(3,i) - zz_bh * taij_rs(2,i) ;
2002 
2003  jbh_rs_x.std_spectral_base() ;
2004 
2005  Vector jbh_rs_y(mp_bh, CON, mp_bh.get_bvect_cart()) ;
2006  jbh_rs_y.set_etat_qcq() ;
2007  for (int i=1; i<=3; i++)
2008  jbh_rs_y.set(i) = zz_bh * taij_rs(1,i) - xx_bh * taij_rs(3,i) ;
2009 
2010  jbh_rs_y.std_spectral_base() ;
2011 
2012  Vector jbh_rs_z(mp_bh, CON, mp_bh.get_bvect_cart()) ;
2013  jbh_rs_z.set_etat_qcq() ;
2014  for (int i=1; i<=3; i++)
2015  jbh_rs_z.set(i) = xx_bh * taij_rs(2,i) - yy_bh * taij_rs(1,i) ;
2016 
2017  jbh_rs_z.std_spectral_base() ;
2018 
2019  Scalar tmp(mp_bh) ;
2020  tmp = - 2. * mass * mass * cc * pow(confo_bh,7.)
2021  * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
2022  / lapconf_bh / pow(r_are*rr,3.) ;
2023  tmp.std_spectral_base() ;
2024 
2025  //-------------//
2026  // X component //
2027  //-------------//
2028 
2029  //-------------------------------------
2030  // Integration over the BH map
2031  //-------------------------------------
2032 
2033  // Surface integral <- dzpuis should be 0
2034  // --------------------------------------
2035  Scalar sou_bh_sx1(mp_bh) ;
2036  sou_bh_sx1 = jbh_rs_x(1)%ll(1) + jbh_rs_x(2)%ll(2)
2037  + jbh_rs_x(3)%ll(3) ;
2038  sou_bh_sx1.std_spectral_base() ;
2039  sou_bh_sx1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
2040 
2041  Scalar sou_bh_sx2(mp_bh) ;
2042  sou_bh_sx2 = tmp * (yy_bh % ll(3) - zz_bh % ll(2)) ;
2043  sou_bh_sx2.std_spectral_base() ;
2044 
2045  source_bh_surf = sou_bh_sx1 + sou_bh_sx2 ;
2046 
2047  source_bh_surf.std_spectral_base() ;
2048  source_bh_surf.annule_domain(0) ;
2049  source_bh_surf.raccord(1) ;
2050 
2051  integ_bh_s_x = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
2052 
2053  //-------------------------------------
2054  // Integration over the NS map
2055  //-------------------------------------
2056 
2057  // Volume integral <- dzpuis should be 4
2058  // -------------------------------------
2059  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
2060  * (yy_ns*uu(3) - zz_ns*uu(2)) ;
2061 
2062  source_ns_volm.std_spectral_base() ;
2063  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
2064 
2065  integ_ns_v_x = source_ns_volm.integrale() ;
2066 
2067  p_angu_mom_bhns->set(0) = integ_ns_v_x + 0.5*integ_bh_s_x/qpig ;
2068 
2069 
2070  //-------------//
2071  // Y component //
2072  //-------------//
2073 
2074  //-------------------------------------
2075  // Integration over the BH map
2076  //-------------------------------------
2077 
2078  // Surface integral <- dzpuis should be 0
2079  // --------------------------------------
2080  Scalar sou_bh_sy1(mp_bh) ;
2081  sou_bh_sy1 = jbh_rs_y(1)%ll(1) + jbh_rs_y(2)%ll(2)
2082  + jbh_rs_y(3)%ll(3) ;
2083  sou_bh_sy1.std_spectral_base() ;
2084  sou_bh_sy1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
2085 
2086  Scalar sou_bh_sy2(mp_bh) ;
2087  sou_bh_sy2 = tmp * (zz_bh % ll(1) - xx_bh % ll(3)) ;
2088  sou_bh_sy2.std_spectral_base() ;
2089 
2090  source_bh_surf = sou_bh_sy1 + sou_bh_sy2 ;
2091 
2092  source_bh_surf.std_spectral_base() ;
2093  source_bh_surf.annule_domain(0) ;
2094  source_bh_surf.raccord(1) ;
2095 
2096  integ_bh_s_y = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
2097 
2098  //-------------------------------------
2099  // Integration over the NS map
2100  //-------------------------------------
2101 
2102  // Volume integral <- dzpuis should be 4
2103  // -------------------------------------
2104  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
2105  * (zz_ns*uu(1) - xx_ns*uu(3)) ;
2106 
2107  source_ns_volm.std_spectral_base() ;
2108  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
2109 
2110  integ_ns_v_y = source_ns_volm.integrale() ;
2111 
2112  p_angu_mom_bhns->set(1) = integ_ns_v_y + 0.5*integ_bh_s_y/qpig ;
2113 
2114 
2115  //-------------//
2116  // Z component //
2117  //-------------//
2118 
2119  //-------------------------------------
2120  // Integration over the BH map
2121  //-------------------------------------
2122 
2123  // Surface integral <- dzpuis should be 0
2124  // --------------------------------------
2125  Scalar sou_bh_sz1(mp_bh) ;
2126  sou_bh_sz1 = jbh_rs_z(1)%ll(1) + jbh_rs_z(2)%ll(2)
2127  + jbh_rs_z(3)%ll(3) ;
2128  sou_bh_sz1.std_spectral_base() ;
2129  sou_bh_sz1.dec_dzpuis(2) ; // dzpuis : 2 -> 0
2130 
2131  Scalar sou_bh_sz2(mp_bh) ;
2132  sou_bh_sz2 = tmp * (xx_bh % ll(2) - yy_bh % ll(1)) ;
2133  sou_bh_sz2.std_spectral_base() ;
2134 
2135  source_bh_surf = sou_bh_sz1 + sou_bh_sz2 ;
2136 
2137  source_bh_surf.std_spectral_base() ;
2138  source_bh_surf.annule_domain(0) ;
2139  source_bh_surf.raccord(1) ;
2140 
2141  integ_bh_s_z = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
2142 
2143  //-------------------------------------
2144  // Integration over the NS map
2145  //-------------------------------------
2146 
2147  // Volume integral <- dzpuis should be 4
2148  // -------------------------------------
2149  source_ns_volm = pow(confo_ns, 10.) * (ee + pp)
2150  * (xx_ns*uu(2) - yy_ns*uu(1)) ;
2151 
2152  source_ns_volm.std_spectral_base() ;
2153  source_ns_volm.inc_dzpuis(4) ; // dzpuis : 0 -> 4
2154 
2155  integ_ns_v_z = source_ns_volm.integrale() ;
2156 
2157  p_angu_mom_bhns->set(2) = integ_ns_v_z + 0.5*integ_bh_s_z/qpig ;
2158 
2159  }
2160 
2161  }
2162 
2163  return *p_angu_mom_bhns ;
2164 
2165 }
2166 
2167 
2168  //-----------------------------------------------//
2169  // Error on the virial theorem //
2170  //-----------------------------------------------//
2171 
2173 
2174  if (p_virial_bhns_surf == 0x0) { // a new computation is required
2175 
2176  double virial = 1. - mass_kom_bhns_surf() / mass_adm_bhns_surf() ;
2177 
2178  p_virial_bhns_surf = new double( virial ) ;
2179 
2180  }
2181 
2182  return *p_virial_bhns_surf ;
2183 
2184 }
2185 
2186 
2188 
2189  if (p_virial_bhns_vol == 0x0) { // a new computation is required
2190 
2191  double virial = 1. - mass_kom_bhns_vol() / mass_adm_bhns_vol() ;
2192 
2193  p_virial_bhns_vol = new double( virial ) ;
2194 
2195  }
2196 
2197  return *p_virial_bhns_vol ;
2198 
2199 }
2200 
2201  //--------------------------------------------------//
2202  // X coordinate of the barycenter //
2203  //--------------------------------------------------//
2204 
2205 double Bin_bhns::xa_barycenter() const {
2206 
2207  using namespace Unites ;
2208 
2209  if (p_xa_barycenter == 0x0) { // a new computation is required
2210 
2211  double mass_b = star.mass_b_bhns(hole.is_kerrschild(),
2212  hole.get_mass_bh(), separ) ;
2213 
2214  const Map& mp = star.get_mp() ;
2215  Scalar xxa(mp) ;
2216  xxa = mp.xa ; // Absolute X coordinate
2217  xxa.std_spectral_base() ;
2218 
2219  Scalar tmp(mp) ;
2220 
2221  if (hole.is_kerrschild()) {
2222 
2223  Scalar xx(mp) ;
2224  xx = mp.x ;
2225  xx.std_spectral_base() ;
2226  Scalar yy(mp) ;
2227  yy = mp.y ;
2228  yy.std_spectral_base() ;
2229  Scalar zz(mp) ;
2230  zz = mp.z ;
2231  zz.std_spectral_base() ;
2232 
2233  double yns = (star.get_mp()).get_ori_y() ;
2234 
2235  Scalar rbh(mp) ;
2236  rbh = sqrt( (xx+separ)*(xx+separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
2237  rbh.std_spectral_base() ;
2238 
2239  double mass = ggrav * hole.get_mass_bh() ;
2240 
2241  tmp = sqrt( 1.+2.*mass/rbh ) ;
2242  tmp.std_spectral_base() ;
2243 
2244  }
2245  else { // Isotropic coordinates with the maximal slicing
2246 
2247  tmp = 1. ;
2248  tmp.std_spectral_base() ;
2249 
2250  }
2251 
2252  Scalar confo = star.get_confo_tot() ;
2253  confo.std_spectral_base() ;
2254 
2255  Scalar g_euler = star.get_gam_euler() ;
2256  g_euler.std_spectral_base() ;
2257 
2258  Scalar nbary = star.get_nbar() ;
2259  nbary.std_spectral_base() ;
2260 
2261  Scalar dens = pow(confo, 6.) * g_euler * nbary * tmp * xxa ;
2262  dens.std_spectral_base() ;
2263  dens.inc_dzpuis(4) ;
2264  assert(dens.get_dzpuis() == 4) ;
2265 
2266  double xa_bary = dens.integrale() / mass_b ;
2267 
2268  p_xa_barycenter = new double( xa_bary ) ;
2269 
2270  }
2271 
2272  return *p_xa_barycenter ;
2273 
2274 }
2275 
2276 
2277  //--------------------------------------------------//
2278  // Y coordinate of the barycenter //
2279  //--------------------------------------------------//
2280 
2281 double Bin_bhns::ya_barycenter() const {
2282 
2283  using namespace Unites ;
2284 
2285  if (p_ya_barycenter == 0x0) { // a new computation is required
2286 
2287  double mass_b = star.mass_b_bhns(hole.is_kerrschild(),
2288  hole.get_mass_bh(), separ) ;
2289 
2290  const Map& mp = star.get_mp() ;
2291  Scalar yya(mp) ;
2292  yya = mp.ya ; // Absolute Y coordinate
2293  yya.std_spectral_base() ;
2294 
2295  Scalar tmp(mp) ;
2296 
2297  if (hole.is_kerrschild()) {
2298 
2299  Scalar xx(mp) ;
2300  xx = mp.x ;
2301  xx.std_spectral_base() ;
2302  Scalar yy(mp) ;
2303  yy = mp.y ;
2304  yy.std_spectral_base() ;
2305  Scalar zz(mp) ;
2306  zz = mp.z ;
2307  zz.std_spectral_base() ;
2308 
2309  double yns = (star.get_mp()).get_ori_y() ;
2310 
2311  Scalar rbh(mp) ;
2312  rbh = sqrt( (xx+separ)*(xx+separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
2313  rbh.std_spectral_base() ;
2314 
2315  double mass = ggrav * hole.get_mass_bh() ;
2316 
2317  tmp = sqrt( 1.+2.*mass/rbh ) ;
2318  tmp.std_spectral_base() ;
2319 
2320  }
2321  else { // Isotropic coordinates with the maximal slicing
2322 
2323  tmp = 1. ;
2324  tmp.std_spectral_base() ;
2325 
2326  }
2327 
2328  Scalar confo = star.get_confo_tot() ;
2329  confo.std_spectral_base() ;
2330 
2331  Scalar g_euler = star.get_gam_euler() ;
2332  g_euler.std_spectral_base() ;
2333 
2334  Scalar nbary = star.get_nbar() ;
2335  nbary.std_spectral_base() ;
2336 
2337  Scalar dens = pow(confo, 6.) * g_euler * nbary * tmp * yya ;
2338  dens.std_spectral_base() ;
2339  dens.inc_dzpuis(4) ;
2340  assert(dens.get_dzpuis() == 4) ;
2341 
2342  double ya_bary = dens.integrale() / mass_b ;
2343 
2344  p_ya_barycenter = new double( ya_bary ) ;
2345 
2346  }
2347 
2348  return *p_ya_barycenter ;
2349 
2350 }
2351 }
double * p_mass_adm_bhns_surf
Total ADM mass of the system calculated by the surface integral at infinity.
Definition: bin_bhns.h:97
double y_rot
Absolute Y coordinate of the rotation axis.
Definition: bin_bhns.h:89
double * p_mass_kom_bhns_surf
Total Komar mass of the system calculated by the surface integral at infinity.
Definition: bin_bhns.h:107
Hole_bhns hole
Black hole.
Definition: bin_bhns.h:72
double mass_kom_bhns_surf() const
Total Komar mass.
double ya_barycenter() const
Absolute coordinate Y of the barycenter of the baryon density.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_bhns.h:80
double * p_mass_adm_bhns_vol
Total ADM mass of the system calculated by the volume integral and the surface integral at the appare...
Definition: bin_bhns.h:102
double * p_virial_bhns_vol
Virial theorem error calculated by the ADM mass and the Komar mass of the volume integral.
Definition: bin_bhns.h:128
const Tbl & line_mom_bhns() const
Total linear momentum.
double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density.
double separ
Absolute orbital separation between two centers of BH and NS.
Definition: bin_bhns.h:83
Tbl * p_angu_mom_bhns
Total angular momentum of the system.
Definition: bin_bhns.h:118
double virial_bhns_surf() const
Estimates the relative error on the virial theorem $|1 - M_{\rm Komar} / M_{\rm ADM}|$.
const Tbl & angu_mom_bhns() const
Total angular momentum.
double virial_bhns_vol() const
Estimates the relative error on the virial theorem $|1 - M_{\rm Komar} / M_{\rm ADM}|$.
double * p_ya_barycenter
Absolute coordinate Y of the barycenter of the baryon density.
Definition: bin_bhns.h:134
double x_rot
Absolute X coordinate of the rotation axis.
Definition: bin_bhns.h:86
double * p_virial_bhns_surf
Virial theorem error calculated by the ADM mass and the Komar mass of the surface integral at infinit...
Definition: bin_bhns.h:123
double * p_mass_kom_bhns_vol
Total Komar mass of the system calculated by the volume integral and the surface integral at the appa...
Definition: bin_bhns.h:112
Star_bhns star
Neutron star.
Definition: bin_bhns.h:75
double * p_xa_barycenter
Absolute coordinate X of the barycenter of the baryon density.
Definition: bin_bhns.h:131
Tbl * p_line_mom_bhns
Total linear momentum of the system.
Definition: bin_bhns.h:115
double mass_adm_bhns_surf() const
Total ADM mass.
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
bool is_kerrschild() const
Returns true for a Kerr-Schild background, false for a Conformally flat one.
Definition: blackhole.h:218
const Map & get_mp() const
Returns the mapping.
Definition: blackhole.h:213
double get_mass_bh() const
Returns the gravitational mass of BH [{\tt m_unit}].
Definition: blackhole.h:221
double integrale() const
Computes the integral over all space of *this .
Definition: cmp_integ.C:55
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: cmp_raccord.C:170
void inc_dzpuis()
Increases by the value of dzpuis and changes accordingly the values of the Cmp in the external compac...
Definition: cmp_r_manip.C:166
const Scalar & get_confo_auto_rs() const
Returns the part of the conformal factor from numerical computation.
Definition: hole_bhns.h:434
const Scalar & get_taij_quad_auto() const
Returns the part of rs generated by the black hole.
Definition: hole_bhns.h:494
const Scalar & get_lapconf_comp() const
Returns the part of the lapconf function generated by the companion star.
Definition: hole_bhns.h:375
bool has_bc_lapconf_nd() const
Returns true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH.
Definition: hole_bhns.h:349
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion star.
Definition: hole_bhns.h:462
const Vector & get_d_lapconf_comp() const
Returns the derivative of the lapconf function generated by the companion star.
Definition: hole_bhns.h:402
const Vector & get_shift_comp() const
Returns the part of the shift vector generated by the companion star.
Definition: hole_bhns.h:413
const Scalar & get_taij_quad_tot_rs() const
Returns the part of rs.
Definition: hole_bhns.h:486
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition: hole_bhns.h:447
const Tensor & get_d_shift_comp() const
Returns the derivative of the shift vector generated by the companion star.
Definition: hole_bhns.h:431
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the black hole.
Definition: hole_bhns.h:408
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the black hole.
Definition: hole_bhns.h:439
const Vector & get_shift_auto_rs() const
Returns the part of the shift vector from numerical computation.
Definition: hole_bhns.h:405
bool has_bc_lapconf_fs() const
Returns true for the first type BC for the lapconf function, false for the second type BH.
Definition: hole_bhns.h:354
const Vector & get_d_confo_auto_rs() const
Returns the derivative of the conformal factor generated by the black hole.
Definition: hole_bhns.h:452
const Vector & get_d_lapconf_auto_rs() const
Returns the derivative of the lapconf function generated by the black hole.
Definition: hole_bhns.h:391
const Scalar & get_taij_quad_tot_rot() const
Returns the part of rot.
Definition: hole_bhns.h:488
const Scalar & get_lapconf_auto_rs() const
Returns the part of the lapconf function from numerical computation.
Definition: hole_bhns.h:365
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition: hole_bhns.h:378
const Sym_tensor & get_taij_tot() const
Returns the total extrinsic curvature tensor.
Definition: hole_bhns.h:468
const Scalar & get_confo_comp() const
Returns the part of the conformal factor generated by the companion star.
Definition: hole_bhns.h:444
const Sym_tensor & get_taij_tot_rs() const
Returns the part of rs of the extrinsic curvature tensor.
Definition: hole_bhns.h:465
const Scalar & get_taij_quad_comp() const
Returns the part of rs generated by the companion star.
Definition: hole_bhns.h:497
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 .
Base class for coordinate mappings.
Definition: map.h:670
Coord cosp
Definition: map.h:724
Coord y
y coordinate centered on the grid
Definition: map.h:727
Coord sint
Definition: map.h:721
Coord ya
Absolute y coordinate.
Definition: map.h:731
Coord r
r coordinate centered on the grid
Definition: map.h:718
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Definition: map.C:253
Coord za
Absolute z coordinate.
Definition: map.h:732
Coord sinp
Definition: map.h:723
Coord z
z coordinate centered on the grid
Definition: map.h:728
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Coord x
x coordinate centered on the grid
Definition: map.h:726
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 xa
Absolute x coordinate.
Definition: map.h:730
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 & dsdz() const
Returns of *this , where .
Definition: scalar_deriv.C:328
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
const Scalar & dsdy() const
Returns of *this , where .
Definition: scalar_deriv.C:297
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 .
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 void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
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...
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition: scalar_import.C:68
const Scalar & get_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
Definition: star_bhns.h:339
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the star.
Definition: star_bhns.h:383
const Vector & get_d_confo_auto() const
Returns the derivative of the conformal factor generated by the star.
Definition: star_bhns.h:396
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition: star_bhns.h:347
const Scalar & get_taij_quad_auto() const
Returns the part of the scalar generated by .
Definition: star_bhns.h:415
const Vector & get_d_lapconf_auto() const
Returns the derivative of the lapse function generated by the star.
Definition: star_bhns.h:356
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion black hole.
Definition: star_bhns.h:401
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the star.
Definition: star_bhns.h:364
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition: star_bhns.h:391
const Scalar & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Definition: star.h:379
const Scalar & get_press() const
Returns the fluid pressure.
Definition: star.h:373
const Scalar & get_gam_euler() const
Returns the Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:382
const Vector & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:385
const Scalar & get_nbar() const
Returns the proper baryon density.
Definition: star.h:367
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: star.h:376
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
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 handling.
Definition: tensor.h:288
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition: tensor.C:671
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.