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 $" ;
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]"
182double 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) ;
204 source_bh_surf.set_etat_qcq() ;
206 Scalar source_bh_volm(mp_bh) ;
207 source_bh_volm.set_etat_qcq() ;
209 Scalar source_ns_volm(mp_ns) ;
210 source_ns_volm.set_etat_qcq() ;
214 rr.std_spectral_base() ;
217 st.std_spectral_base() ;
220 ct.std_spectral_base() ;
223 sp.std_spectral_base() ;
226 cp.std_spectral_base() ;
230 ll.set(1) = st % cp ;
231 ll.set(2) = st % sp ;
233 ll.std_spectral_base() ;
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() ;
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() ;
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) ;
254 llshift_auto.std_spectral_base() ;
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 ;
443 Scalar r_are(mp_bh) ;
450 Scalar divshift_zero(divshift) ;
451 divshift_zero.dec_dzpuis(2) ;
453 Scalar lldllsh_zero(lldllsh) ;
454 lldllsh_zero.dec_dzpuis(2) ;
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.) ;
463 source_bh_surf.std_spectral_base() ;
464 source_bh_surf.annule_domain(0) ;
465 source_bh_surf.raccord(1) ;
467 integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
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() ;
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.)
485 source_bh_volm.std_spectral_base() ;
486 source_bh_volm.annule_domain(0) ;
488 integ_bh_v = source_bh_volm.
integrale() ;
496 Scalar sou_ns1(mp_ns) ;
497 sou_ns1 =
pow(confo_ns, 5.) * ener_euler ;
499 sou_ns1.inc_dzpuis(4) ;
501 source_ns_volm = 0.25 * taij_quad_auto_ns
502 /
pow(confo_ns_auto+0.5, 7.) / qpig + sou_ns1 ;
504 source_ns_volm.std_spectral_base() ;
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]"
657double 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) ;
679 source_bh_surf.set_etat_qcq() ;
681 Scalar source_bh_volm(mp_bh) ;
682 source_bh_volm.set_etat_qcq() ;
684 Scalar source_ns_volm(mp_ns) ;
685 source_ns_volm.set_etat_qcq() ;
689 rr.std_spectral_base() ;
692 st.std_spectral_base() ;
695 ct.std_spectral_base() ;
698 sp.std_spectral_base() ;
701 cp.std_spectral_base() ;
705 ll.set(1) = st % cp ;
706 ll.set(2) = st % sp ;
708 ll.std_spectral_base() ;
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() ;
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) ;
737 llshift_auto.std_spectral_base() ;
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 ;
962 Scalar r_are(mp_bh) ;
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) ;
973 lldlapconf_is.std_spectral_base() ;
977 Scalar divshift_zero(divshift) ;
978 divshift_zero.dec_dzpuis(2) ;
980 Scalar lldllsh_zero(lldllsh) ;
981 lldllsh_zero.dec_dzpuis(2) ;
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.) ;
991 sou_bh_s_psi.std_spectral_base() ;
992 sou_bh_s_psi.annule_domain(0) ;
993 sou_bh_s_psi.raccord(1) ;
998 source_bh_surf = sou_bh_s_psi ;
1000 source_bh_surf.std_spectral_base() ;
1001 source_bh_surf.annule_domain(0) ;
1002 source_bh_surf.raccord(1) ;
1007 Scalar sou_bh_s1(mp_bh) ;
1008 sou_bh_s1 = 0.5 * lapconf_bh / rr ;
1010 sou_bh_s1.annule_domain(0) ;
1011 sou_bh_s1.raccord(1) ;
1013 source_bh_surf = sou_bh_s1 + sou_bh_s_psi ;
1015 source_bh_surf.std_spectral_base() ;
1016 source_bh_surf.annule_domain(0) ;
1017 source_bh_surf.raccord(1) ;
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) ;
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 ;
1034 sou_bh_s2.std_spectral_base() ;
1036 source_bh_surf = sou_bh_s1 + sou_bh_s2 + sou_bh_s_psi ;
1038 source_bh_surf.std_spectral_base() ;
1039 source_bh_surf.annule_domain(0) ;
1044 integ_bh_s = mp_aff.integrale_surface(source_bh_surf, radius_ah) ;
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.) ;
1055 sou_bh1.std_spectral_base() ;
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 ;
1062 sou_bh2.std_spectral_base() ;
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 ;
1074 sou_bh3.std_spectral_base() ;
1076 source_bh_volm = sou_bh1 + sou_bh2 + sou_bh3 ;
1077 source_bh_volm.std_spectral_base() ;
1078 source_bh_volm.annule_domain(0) ;
1080 integ_bh_v = source_bh_volm.integrale() ;
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 ;
1092 sou_ns1.inc_dzpuis(4) ;
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() ;
1100 source_ns_volm = sou_ns1 + sou_ns2 ;
1101 source_ns_volm.std_spectral_base() ;
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) ;
1709 for (
int i=1; i<=3; i++)
1710 jbh_y.
set(i) = zz_bh * taij(1,i) - xx_bh * taij(3,i) ;
1716 for (
int i=1; i<=3; i++)
1717 jbh_z.
set(i) = xx_bh * taij(2,i) - yy_bh * taij(1,i) ;
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) ;
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) ;
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) ;
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) ;
1852 jbh_y_ll.std_spectral_base() ;
1853 jbh_y_ll.dec_dzpuis(2) ;
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) ;
1909 jbh_z_ll.std_spectral_base() ;
1910 jbh_z_ll.dec_dzpuis(2) ;
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) ;
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) ;
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) ;
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 ;
double * p_mass_adm_bhns_surf
Total ADM mass of the system calculated by the surface integral at infinity.
double y_rot
Absolute Y coordinate of the rotation axis.
double * p_mass_kom_bhns_surf
Total Komar mass of the system calculated by the surface integral at infinity.
Hole_bhns hole
Black hole.
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.
double * p_mass_adm_bhns_vol
Total ADM mass of the system calculated by the volume integral and the surface integral at the appare...
double * p_virial_bhns_vol
Virial theorem error calculated by the ADM mass and the Komar mass of the volume integral.
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.
Tbl * p_angu_mom_bhns
Total angular momentum of the system.
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.
double x_rot
Absolute X coordinate of the rotation axis.
double * p_virial_bhns_surf
Virial theorem error calculated by the ADM mass and the Komar mass of the surface integral at infinit...
double * p_mass_kom_bhns_vol
Total Komar mass of the system calculated by the volume integral and the surface integral at the appa...
Star_bhns star
Neutron star.
double * p_xa_barycenter
Absolute coordinate X of the barycenter of the baryon density.
Tbl * p_line_mom_bhns
Total linear momentum of the system.
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.
const Map & get_mp() const
Returns the mapping.
double get_mass_bh() const
Returns the gravitational mass of BH [{\tt m_unit}].
double integrale() const
Computes the integral over all space of *this .
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.
const Scalar & get_taij_quad_tot_rot() const
Returns the part of rot.
const Vector & get_d_confo_auto_rs() const
Returns the derivative of the conformal factor generated by the black hole.
const Scalar & get_taij_quad_comp() const
Returns the part of rs generated by the companion star.
const Vector & get_d_lapconf_auto_rs() const
Returns the derivative of the lapconf function generated by the black hole.
const Scalar & get_lapconf_comp() const
Returns the part of the lapconf function generated by the companion star.
const Tensor & get_d_shift_comp() const
Returns the derivative of the shift vector generated by the companion star.
const Scalar & get_confo_auto_rs() const
Returns the part of the conformal factor from numerical computation.
const Scalar & get_lapconf_auto_rs() const
Returns the part of the lapconf function from numerical computation.
const Scalar & get_taij_quad_tot_rs() const
Returns the part of rs.
const Vector & get_shift_auto_rs() const
Returns the part of the shift vector from numerical computation.
const Vector & get_d_lapconf_comp() const
Returns the derivative of the lapconf function generated by the companion star.
const Sym_tensor & get_taij_tot() const
Returns the total extrinsic curvature tensor.
bool has_bc_lapconf_fs() const
Returns true for the first type BC for the lapconf function, false for the second type BH.
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the black hole.
const Scalar & get_confo_tot() const
Returns the total conformal factor.
const Scalar & get_confo_comp() const
Returns the part of the conformal factor generated by the companion star.
const Scalar & get_taij_quad_auto() const
Returns the part of rs generated by the black hole.
const Vector & get_shift_comp() const
Returns the part of the shift vector generated by the companion star.
const Sym_tensor & get_taij_tot_rs() const
Returns the part of rs of the extrinsic curvature tensor.
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the black hole.
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion star.
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.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord y
y coordinate centered on the grid
Coord ya
Absolute y coordinate.
Coord r
r coordinate centered on the grid
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Coord za
Absolute z coordinate.
Coord z
z coordinate centered on the grid
double get_ori_x() const
Returns the x coordinate of the origin.
Coord x
x coordinate centered on the grid
Coord xa
Absolute x coordinate.
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).
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.
const Scalar & dsdy() const
Returns of *this , where .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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.
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.
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 Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the star.
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion black hole.
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the star.
const Scalar & get_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
const Scalar & get_taij_quad_auto() const
Returns the part of the scalar generated by .
const Vector & get_d_confo_auto() const
Returns the derivative of the conformal factor generated by the star.
const Vector & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
const Scalar & get_gam_euler() const
Returns the Lorentz factor between the fluid and Eulerian observers.
const Map & get_mp() const
Returns the mapping.
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
const Scalar & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
const Scalar & get_press() const
Returns the fluid pressure.
const Scalar & get_nbar() const
Returns the proper baryon density.
Class intended to describe valence-2 symmetric tensors.
void annule_hard()
Sets the Tbl to zero in a hard way.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Scalar & set(int)
Read/write access to a component.
Cmp sqrt(const Cmp &)
Square root.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Standard units of space, time and mass.