28 char bin_bhns_global_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_global.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
58 #include "blackhole.h" 60 #include "utilitaires.h" 92 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
110 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
122 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
128 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
143 Scalar lldconf_iso = confo_bh_auto_rs.
dsdr() ;
147 anoth = 0.5 *
sqrt(r_are)
148 * (
sqrt(1. -2.*mass/r_are/rr + cc*cc*
pow(mass/r_are/rr,4.))
164 cout <<
"ADM mass (surface) : " << madm / msol <<
" [Mo]" 182 double Bin_bhns::mass_adm_bhns_vol()
const {
198 double radius_ah = mp_bh.
val_r(1,-1.,M_PI/2.,0.) ;
203 Scalar source_bh_surf(mp_bh) ;
206 Scalar source_bh_volm(mp_bh) ;
209 Scalar source_ns_volm(mp_ns) ;
230 ll.set(1) = st % cp ;
231 ll.
set(2) = st % sp ;
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) ;
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)) ;
251 Scalar llshift_auto(mp_bh) ;
252 llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
253 + ll(3)%shift_auto_rs(3) ;
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)) ;
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) ;
292 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
419 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
431 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
437 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
450 Scalar divshift_zero(divshift) ;
453 Scalar lldllsh_zero(lldllsh) ;
456 source_bh_surf = confo_bh / rr
457 -
pow(confo_bh, 4.) * (divshift_zero - 3.*lldllsh_zero)
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.) ;
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.) ;
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.)
488 integ_bh_v = source_bh_volm.
integrale() ;
497 sou_ns1 =
pow(confo_ns, 5.) * ener_euler ;
501 source_ns_volm = 0.25 * taij_quad_auto_ns
502 /
pow(confo_ns_auto+0.5, 7.) / qpig + sou_ns1 ;
506 integ_ns_v = source_ns_volm.
integrale() ;
508 cout <<
"integ_bh_s : " << integ_bh_s/ qpig / msol
510 << integ_bh_v/ qpig / msol
511 <<
" integ_ns_v : " << integ_ns_v/ msol << endl ;
516 madm = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
518 cout <<
"ADM mass (volume) : " << madm / msol <<
" [Mo]" 561 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
579 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
591 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
597 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
614 lldlap_bh = lapconf_bh_auto_rs.
dsdr()
615 - confo_bh_auto_rs.
dsdr() ;
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 ;
631 lldlap_ns = lapconf_ns_auto.
dsdr() - confo_ns_auto.
dsdr() ;
638 cout <<
"Komar mass (surface) : " << mkom / msol <<
" [Mo]" 657 double Bin_bhns::mass_kom_bhns_vol()
const {
673 double radius_ah = mp_bh.
val_r(1,-1.,M_PI/2.,0.) ;
678 Scalar source_bh_surf(mp_bh) ;
681 Scalar source_bh_volm(mp_bh) ;
684 Scalar source_ns_volm(mp_ns) ;
705 ll.set(1) = st % cp ;
706 ll.
set(2) = st % sp ;
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) ;
734 Scalar llshift_auto(mp_bh) ;
735 llshift_auto = ll(1)%shift_auto_rs(1) + ll(2)%shift_auto_rs(2)
736 + ll(3)%shift_auto_rs(3) ;
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)) ;
801 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
938 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
950 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
956 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
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) ;
977 Scalar divshift_zero(divshift) ;
980 Scalar lldllsh_zero(lldllsh) ;
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)
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.) ;
998 source_bh_surf = sou_bh_s_psi ;
1007 Scalar sou_bh_s1(mp_bh) ;
1008 sou_bh_s1 = 0.5 * lapconf_bh / rr ;
1013 source_bh_surf = sou_bh_s1 + sou_bh_s_psi ;
1023 Scalar sou_bh_s1(mp_bh) ;
1024 sou_bh_s1 = lldlapconf_is ;
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 ;
1036 source_bh_surf = sou_bh_s1 + sou_bh_s2 + sou_bh_s_psi ;
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.) ;
1059 sou_bh2 = 0.125 * (7.*lapconf_bh/
pow(confo_bh,8.)
1060 + 1./
pow(confo_bh,7.)) * taij_quad_auto_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 ;
1076 source_bh_volm = sou_bh1 + sou_bh2 + sou_bh3 ;
1080 integ_bh_v = source_bh_volm.
integrale() ;
1089 sou_ns1 = 0.5 *
pow(confo_ns,4.) * (lapconf_ns + confo_ns)
1090 * ener_euler + lapconf_ns *
pow(confo_ns,4.) * s_euler ;
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 ;
1100 source_ns_volm = sou_ns1 + sou_ns2 ;
1103 integ_ns_v = source_ns_volm.
integrale() ;
1105 cout <<
"integ_bh_s : " << integ_bh_s/ qpig / msol
1107 << integ_bh_v/ qpig / msol
1108 <<
" integ_ns_v : " << integ_ns_v/ msol << endl ;
1113 mkom = (integ_bh_s + integ_bh_v) / qpig + integ_ns_v ;
1115 cout <<
"Komar mass (volume) : " << mkom / msol <<
" [Mo]" 1148 int nz = (
hole.
get_mp()).get_mg()->get_nzone() ;
1149 double* bornes =
new double[nz+1] ;
1150 double radius =
separ + 2. ;
1152 for (
int i=nz-1; i>0; i--) {
1153 bornes[i] = radius ;
1157 bornes[nz] = __infinity ;
1180 Vector shift_tot = shift_bh + shift_ns ;
1182 shift_tot.
annule(0, nz-2) ;
1199 lc.set(1) = stc * cpc ;
1200 lc.
set(2) = stc * spc ;
1220 rl =
sqrt(1. - 2.*xbhsr*lc(1) - 2.*ybhsr*lc(2)
1221 + xbhsr*xbhsr + ybhsr*ybhsr) ;
1226 ll.set(1) = (lc(1) - xbhsr) / rl ;
1227 ll.
set(2) = (lc(2) - ybhsr)/ rl ;
1228 ll.
set(3) = lc(3) / rl ;
1232 lcll = lc(1)*ll(1) + lc(2)*ll(2) + lc(3)*ll(3) ;
1235 Scalar divshift(mp_aff) ;
1236 divshift = shift_tot(1).deriv(1) + shift_tot(2).deriv(2)
1237 + shift_tot(3).deriv(3) ;
1241 llshift = ll(1)*shift_tot(1) + ll(2)*shift_tot(2)
1242 + ll(3)*shift_tot(3) ;
1246 lcshift = lc(1)*shift_tot(1) + lc(2)*shift_tot(2)
1247 + lc(3)*shift_tot(3) ;
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)) ;
1259 for (
int i=1; i<=3; i++) {
1260 dshift.
set(i) = shift_tot(i).dsdr() ;
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)) ;
1275 lap_bh2 = 1./(1.+2.*mass/rl/rr) ;
1279 lap2hbh = lap_bh2 * mass / rl / rr ;
1293 kk = 4.*
sqrt(lap_bh2)*lap2hbh*(1.+3.*mass/rl/rr)/3./rl/rr ;
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. ;
1312 kij_x2 = kk * (lc(1) - 2.*lap2hbh*ll(1)*lcll) ;
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)
1320 + 2.*ll(1)*lcll*divshift/3.
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. ;
1333 kij_y1.inc_dzpuis(2) ;
1336 kij_y2 = kk * (lc(2) - 2.*lap2hbh*ll(2)*lcll) ;
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)
1344 + 2.*ll(2)*lcll*divshift/3.
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. ;
1357 kij_z1.inc_dzpuis(2) ;
1360 kij_z2 = kk * (lc(3) - 2.*lap2hbh*ll(3)*lcll) ;
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)
1368 + 2.*ll(3)*lcll*divshift/3.
1463 double radius_ah = mp_bh.
val_r(1,-1.,M_PI/2.,0.) ;
1480 ll.set(1) = st % cp ;
1481 ll.
set(2) = st % sp ;
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) ;
1512 Scalar sou_ns_vx(mp_ns) ;
1514 sou_ns_vx =
pow(confo_ns, 10.) * (ee + pp) * uu(1) ;
1519 double integ_ns_v_x = sou_ns_vx.
integrale() ;
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) ;
1550 Scalar sou_ns_vy(mp_ns) ;
1552 sou_ns_vy =
pow(confo_ns, 10.) * (ee + pp) * uu(2) ;
1557 double integ_ns_v_y = sou_ns_vy.
integrale() ;
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) ;
1589 Scalar sou_ns_vz(mp_ns) ;
1591 sou_ns_vz =
pow(confo_ns, 10.) * (ee + pp) * uu(3) ;
1596 double integ_ns_v_z = sou_ns_vz.
integrale() ;
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 ;
1637 double radius_ah = mp_bh.
val_r(1,-1.,M_PI/2.,0.) ;
1640 Scalar source_bh_surf(mp_bh) ;
1643 Scalar source_bh_volm(mp_bh) ;
1646 Scalar source_ns_volm(mp_ns) ;
1668 ll.set(1) = st % cp ;
1669 ll.
set(2) = st % sp ;
1702 for (
int i=1; i<=3; i++)
1703 jbh_x.set(i) = yy_bh * taij(3,i) - zz_bh * taij(2,i) ;
1705 jbh_x.std_spectral_base() ;
1709 for (
int i=1; i<=3; i++)
1710 jbh_y.set(i) = zz_bh * taij(1,i) - xx_bh * taij(3,i) ;
1712 jbh_y.std_spectral_base() ;
1716 for (
int i=1; i<=3; i++)
1717 jbh_z.set(i) = xx_bh * taij(2,i) - yy_bh * taij(1,i) ;
1719 jbh_z.std_spectral_base() ;
1728 lap_bh2 = 1./(1.+2.*mass/rr) ;
1732 lcnf =
log(confo_bh) ;
1737 for (
int i=1; i<=3; i++)
1738 jbhsr_x.set(i) = ll(2)*taij(3,i) - ll(3)*taij(2,i) ;
1740 jbhsr_x.std_spectral_base() ;
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) ;
1747 jbhsr_y.std_spectral_base() ;
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) ;
1754 jbhsr_z.std_spectral_base() ;
1757 tmp1 = 2. *
pow(lap_bh2,2.5) * mass * (1.+3.*mass/rr)
1758 *
pow(confo_bh,6.) * ori_bh / 3. / rr / rr ;
1763 tmp2 = 4. *
sqrt(lap_bh2) * (1.+3.*mass/rr) *
pow(confo_bh,6.) ;
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)) ;
1775 dlcnf = - 2. * lap_bh2 * tmp2 * mass * lcnf.
dsdr() / rr ;
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 ;
1789 tmp4 = lap_bh2 * mass / rr ;
1802 source_bh_surf = jbh_x(1)*ll(1) + jbh_x(2)*ll(2) + jbh_x(3)*ll(3) ;
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() ) ) ;
1821 integ_bh_v_x = source_bh_volm.
integrale() ;
1829 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
1830 * (yy_ns*uu(3) - zz_ns*uu(2)) ;
1835 integ_ns_v_x = source_ns_volm.
integrale() ;
1838 + 0.5 * (integ_bh_s_x + integ_bh_v_x) / qpig ;
1851 jbh_y_ll = jbh_y(1)*ll(1) + jbh_y(2)*ll(2) + jbh_y(3)*ll(3) ;
1855 source_bh_surf = jbh_y_ll - tmp1 * ll(3) ;
1865 tmp3_llz = tmp3 * ll(3) * ori_bh / rr ;
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() )
1878 integ_bh_v_y = source_bh_volm.
integrale() ;
1886 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
1887 * (zz_ns*uu(1) - xx_ns*uu(3)) ;
1892 integ_ns_v_y = source_ns_volm.
integrale() ;
1895 + 0.5 * (integ_bh_s_y + integ_bh_v_y) / qpig ;
1908 jbh_z_ll = jbh_z(1)*ll(1) + jbh_z(2)*ll(2) + jbh_z(3)*ll(3) ;
1911 source_bh_surf = jbh_z_ll + tmp1 * ll(2) ;
1922 tmp3_lly = tmp3 * ll(2) * ori_bh / rr ;
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() )
1935 integ_bh_v_z = source_bh_volm.
integrale() ;
1943 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
1944 * (xx_ns*uu(2) - yy_ns*uu(1)) ;
1949 integ_ns_v_z = source_ns_volm.
integrale() ;
1952 + 0.5 * (integ_bh_s_z + integ_bh_v_z) / qpig ;
1965 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
1977 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
1983 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
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) ;
2003 jbh_rs_x.std_spectral_base() ;
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) ;
2010 jbh_rs_y.std_spectral_base() ;
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) ;
2017 jbh_rs_z.std_spectral_base() ;
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.) ;
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) ;
2041 Scalar sou_bh_sx2(mp_bh) ;
2042 sou_bh_sx2 = tmp * (yy_bh % ll(3) - zz_bh % ll(2)) ;
2045 source_bh_surf = sou_bh_sx1 + sou_bh_sx2 ;
2059 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
2060 * (yy_ns*uu(3) - zz_ns*uu(2)) ;
2065 integ_ns_v_x = source_ns_volm.
integrale() ;
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) ;
2086 Scalar sou_bh_sy2(mp_bh) ;
2087 sou_bh_sy2 = tmp * (zz_bh % ll(1) - xx_bh % ll(3)) ;
2090 source_bh_surf = sou_bh_sy1 + sou_bh_sy2 ;
2104 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
2105 * (zz_ns*uu(1) - xx_ns*uu(3)) ;
2110 integ_ns_v_y = source_ns_volm.
integrale() ;
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) ;
2131 Scalar sou_bh_sz2(mp_bh) ;
2132 sou_bh_sz2 = tmp * (xx_bh % ll(2) - yy_bh % ll(1)) ;
2135 source_bh_surf = sou_bh_sz1 + sou_bh_sz2 ;
2149 source_ns_volm =
pow(confo_ns, 10.) * (ee + pp)
2150 * (xx_ns*uu(2) - yy_ns*uu(1)) ;
2155 integ_ns_v_z = source_ns_volm.
integrale() ;
2191 double virial = 1. - mass_kom_bhns_vol() / mass_adm_bhns_vol() ;
2236 rbh =
sqrt( (xx+
separ)*(xx+
separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
2241 tmp =
sqrt( 1.+2.*mass/rbh ) ;
2261 Scalar dens =
pow(confo, 6.) * g_euler * nbary * tmp * xxa ;
2266 double xa_bary = dens.
integrale() / mass_b ;
2312 rbh =
sqrt( (xx+
separ)*(xx+
separ) + (yy+yns)*(yy+yns) + zz*zz ) ;
2317 tmp =
sqrt( 1.+2.*mass/rbh ) ;
2337 Scalar dens =
pow(confo, 6.) * g_euler * nbary * tmp * yya ;
2342 double ya_bary = dens.
integrale() / mass_b ;
Coord xa
Absolute x coordinate.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
double get_mass_bh() const
Returns the gravitational mass of BH [{ m_unit}].
Cmp log(const Cmp &)
Neperian logarithm.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Hole_bhns hole
Black hole.
const Tbl & line_mom_bhns() const
Total linear momentum.
const Scalar & get_press() const
Returns the fluid pressure.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
double ya_barycenter() const
Absolute coordinate Y of the barycenter of the baryon density.
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the black hole.
const Vector & get_d_confo_auto_rs() const
Returns the derivative of the conformal factor generated by the black hole.
const Vector & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Cmp sqrt(const Cmp &)
Square root.
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Standard units of space, time and mass.
double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density.
bool is_kerrschild() const
Returns true for a Kerr-Schild background, false for a Conformally flat one.
double & set(int i)
Read/write of a particular element (index i) (1D case)
const Scalar & dsdz() const
Returns of *this , where .
Tensor field of valence 0 (or component of a tensorial field).
Base class for coordinate mappings.
double get_ori_x() const
Returns the x coordinate of the origin.
double integrale() const
Computes the integral over all space of *this .
Tbl * p_angu_mom_bhns
Total angular momentum of the system.
Tbl * p_line_mom_bhns
Total linear momentum of the system.
bool has_bc_lapconf_fs() const
Returns true for the first type BC for the lapconf function, false for the second type BH...
double * p_mass_kom_bhns_surf
Total Komar mass of the system calculated by the surface integral at infinity.
const Scalar & get_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
const Scalar & get_confo_comp() const
Returns the part of the conformal factor generated by the companion star.
const Scalar & dsdx() const
Returns of *this , where .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
bool has_bc_lapconf_nd() const
Returns true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH...
double virial_bhns_vol() const
Estimates the relative error on the virial theorem $|1 - M_{ Komar} / M_{ ADM}|$. ...
const Map & get_mp() const
Returns the mapping.
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Vector & get_d_confo_auto() const
Returns the derivative of the conformal factor generated by the star.
const Sym_tensor & get_taij_tot_rs() const
Returns the part of rs of the extrinsic curvature tensor.
Tensor field of valence 1.
double mass_kom_bhns_surf() const
Total Komar mass.
double * p_mass_kom_bhns_vol
Total Komar mass of the system calculated by the volume integral and the surface integral at the appa...
double separ
Absolute orbital separation between two centers of BH and NS.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
double * p_xa_barycenter
Absolute coordinate X of the barycenter of the baryon density.
double omega
Angular velocity with respect to an asymptotically inertial observer.
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.
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the star.
const Scalar & get_lapconf_comp() const
Returns the part of the lapconf function generated by the companion star.
const Scalar & get_taij_quad_tot_rot() const
Returns the part of rot.
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 & dsdy() const
Returns of *this , where .
int get_dzpuis() const
Returns dzpuis.
const Vector & get_shift_auto_rs() const
Returns the part of the shift vector from numerical computation.
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Star_bhns star
Neutron star.
const Scalar & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
const Sym_tensor & get_taij_tot() const
Returns the total extrinsic curvature tensor.
const Map & get_mp() const
Returns the mapping.
double mass_adm_bhns_surf() const
Total ADM mass.
const Scalar & get_taij_quad_tot_rs() const
Returns the part of rs.
double * p_mass_adm_bhns_vol
Total ADM mass of the system calculated by the volume integral and the surface integral at the appare...
Cmp pow(const Cmp &, int)
Power .
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the star.
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the black hole.
const Scalar & get_taij_quad_auto() const
Returns the part of rs generated by the black hole.
const Vector & get_d_lapconf_auto() const
Returns the derivative of the lapse function generated by the star.
const Scalar & get_confo_tot() const
Returns the total conformal factor.
const Vector & get_d_lapconf_comp() const
Returns the derivative of the lapconf function generated by the companion star.
Coord ya
Absolute y coordinate.
const Tbl & angu_mom_bhns() const
Total angular momentum.
const Scalar & get_confo_auto_rs() const
Returns the part of the conformal factor from numerical computation.
const Scalar & get_nbar() const
Returns the proper baryon density.
double * p_virial_bhns_vol
Virial theorem error calculated by the ADM mass and the Komar mass of the volume integral.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
double * p_ya_barycenter
Absolute coordinate Y of the barycenter of the baryon density.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Coord y
y coordinate centered on the grid
Coord za
Absolute z coordinate.
const Scalar & get_taij_quad_comp() const
Returns the part of rs generated by the companion star.
const Scalar & dsdr() const
Returns of *this .
double * p_mass_adm_bhns_surf
Total ADM mass of the system calculated by the surface integral at infinity.
Coord x
x coordinate centered on the grid
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion black hole.
double virial_bhns_surf() const
Estimates the relative error on the virial theorem $|1 - M_{ Komar} / M_{ ADM}|$. ...
Scalar & set(int)
Read/write access to a component.
const Scalar & get_gam_euler() const
Returns the Lorentz factor between the fluid and Eulerian observers.
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
double x_rot
Absolute X coordinate of the rotation axis.
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
void annule_hard()
Sets the Tbl to zero in a hard way.
const Vector & get_d_lapconf_auto_rs() const
Returns the derivative of the lapconf function generated by the black hole.
const Scalar & get_lapconf_auto_rs() const
Returns the part of the lapconf function from numerical computation.
Coord z
z coordinate centered on the grid
const Vector & get_shift_comp() const
Returns the part of the shift vector generated by the companion star.
const Scalar & get_taij_quad_auto() const
Returns the part of the scalar generated by .
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...
double * p_virial_bhns_surf
Virial theorem error calculated by the ADM mass and the Komar mass of the surface integral at infinit...
Class intended to describe valence-2 symmetric tensors.
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion star.
double y_rot
Absolute Y coordinate of the rotation axis.
Coord r
r coordinate centered on the grid
const Tensor & get_d_shift_comp() const
Returns the derivative of the shift vector generated by the companion star.