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