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
28char 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
67namespace 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 ;
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
182double 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
657double 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) ;
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
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
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.
void inc_dzpuis()
Increases by the value of dzpuis and changes accordingly the values of the Cmp in the external compac...
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 Scalar & get_taij_quad_tot_rot() const
Returns the part of rot.
Definition hole_bhns.h:488
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 Scalar & get_taij_quad_comp() const
Returns the part of rs generated by the companion star.
Definition hole_bhns.h:497
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_lapconf_comp() const
Returns the part of the lapconf function generated by the companion star.
Definition hole_bhns.h:375
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 Scalar & get_confo_auto_rs() const
Returns the part of the conformal factor from numerical computation.
Definition hole_bhns.h:434
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_taij_quad_tot_rs() const
Returns the part of rs.
Definition hole_bhns.h:486
const Vector & get_shift_auto_rs() const
Returns the part of the shift vector from numerical computation.
Definition hole_bhns.h:405
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 Sym_tensor & get_taij_tot() const
Returns the total extrinsic curvature tensor.
Definition hole_bhns.h:468
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 Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the black hole.
Definition hole_bhns.h:439
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition hole_bhns.h:447
const Scalar & get_confo_comp() const
Returns the part of the conformal factor generated by the companion star.
Definition hole_bhns.h:444
const Scalar & get_taij_quad_auto() const
Returns the part of rs generated by the black hole.
Definition hole_bhns.h:494
const Vector & get_shift_comp() const
Returns the part of the shift vector generated by the companion star.
Definition hole_bhns.h:413
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_lapconf_tot() const
Returns the total lapconf function.
Definition hole_bhns.h:378
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the black hole.
Definition hole_bhns.h:408
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion star.
Definition hole_bhns.h:462
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
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 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
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 .
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
const Scalar & dsdy() const
Returns of *this , where .
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 .
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 .
const Scalar & dsdr() const
Returns of *this .
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.
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition star_bhns.h:347
const Vector & get_d_lapconf_auto() const
Returns the derivative of the lapse function generated by the star.
Definition star_bhns.h:356
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition star_bhns.h:391
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_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_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
Definition star_bhns.h:339
const Scalar & get_taij_quad_auto() const
Returns the part of the scalar generated by .
Definition star_bhns.h:415
const Vector & get_d_confo_auto() const
Returns the derivative of the conformal factor generated by the star.
Definition star_bhns.h:396
const Vector & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition star.h:385
const Scalar & get_gam_euler() const
Returns the Lorentz factor between the fluid and Eulerian observers.
Definition star.h:382
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition star.h:376
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_nbar() const
Returns the proper baryon density.
Definition star.h:367
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.