61 #include "hole_bhns.h" 68 #include "utilitaires.h" 73 int filter_r,
int filter_r_s,
int filter_p_s,
74 double x_rot,
double y_rot,
double precis,
75 double omega_orb,
double resize_bh,
76 const Tbl& fact_resize,
Tbl& diff) {
91 double rr_in_1 =
mp.
val_r(1, -1., M_PI/2, 0.) ;
215 double rr_out_nm3 =
mp.
val_r(nz-3, 1., M_PI/2., 0.) ;
216 mp.
resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(1)) ;
219 double rr_out_nm4 =
mp.
val_r(nz-4, 1., M_PI/2., 0.) ;
220 mp.
resize(nz-4, rr_in_1/rr_out_nm4 * fact_resize(0)) ;
225 double rr_out_nm2 =
mp.
val_r(nz-2, 1., M_PI/2., 0.) ;
226 mp.
resize(nz-2, 2. * rr_in_1 * fact_resize(1) / rr_out_nm2) ;
231 double rr_out_1 =
mp.
val_r(1, 1., M_PI/2., 0.) ;
232 mp.
resize(1, rr_in_1/rr_out_1 * resize_bh) ;
237 double rr_out_nm4_new =
mp.
val_r(nz-4, 1., M_PI/2., 0.) ;
239 for (
int i=1; i<nz-5; i++) {
241 double rr_out_i =
mp.
val_r(i, 1., M_PI/2., 0.) ;
243 double rr_mid = rr_out_i
244 + (rr_out_nm4_new - rr_out_i) /
double(nz - 4 - i) ;
246 double rr_2timesi = 2. * rr_out_i ;
248 if (rr_2timesi < rr_mid) {
250 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
251 mp.
resize(i+1, rr_2timesi / rr_out_ip1) ;
256 double rr_out_ip1 =
mp.
val_r(i+1, 1., M_PI/2., 0.) ;
257 mp.
resize(i+1, rr_mid / rr_out_ip1) ;
283 double& diff_lapconf = diff.
set(0) ;
284 double& diff_confo = diff.
set(1) ;
285 double& diff_shift_x = diff.
set(2) ;
286 double& diff_shift_y = diff.
set(3) ;
287 double& diff_shift_z = diff.
set(4) ;
303 double mass = ggrav *
mass_bh ;
323 ll.set(1) = st % cp ;
324 ll.
set(2) = st % sp ;
329 for (
int i=1; i<=3; i++) {
341 for (
int mer_bh=0; mer_bh<mermax_bh; mer_bh++) {
343 cout <<
"--------------------------------------------------" << endl ;
344 cout <<
"step: " << mer_bh << endl ;
345 cout <<
"diff_lapconf = " << diff_lapconf << endl ;
346 cout <<
"diff_confo = " << diff_confo << endl ;
347 cout <<
"diff_shift : x = " << diff_shift_x
348 <<
" y = " << diff_shift_y <<
" z = " << diff_shift_z << endl ;
352 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
366 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
378 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
384 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
395 lapbh_iso =
sqrt(1. - 2.*mass/r_are/rr
396 + cc*cc*
pow(mass/r_are/rr,4.)) ;
402 psibh_iso =
sqrt(r_are) ;
408 dlapbh_iso = mass/r_are/rr - 2.*cc*cc*
pow(mass/r_are/rr,4.) ;
437 tmpl3 = 5.25 * cc * cc *
pow(mass,4.) * lapbh_iso
446 source_lapconf = tmpl1 + tmpl2 + tmpl3 ;
458 if (source_lapconf.
get_etat() != ETATZERO) {
459 source_lapconf.
filtre(filter_r) ;
465 lapconf_m1.set_etat_qcq() ;
498 Scalar tmpc2 = 0.75 * cc * cc *
pow(mass,4.)
517 source_confo = tmpc1 + tmpc2 + tmpc3 ;
529 if (source_confo.
get_etat() != ETATZERO) {
530 source_confo.
filtre(filter_r) ;
534 bc_conf =
bc_confo(omega_orb, x_rot, y_rot) ;
559 for (
int i=1; i<=3; i++) {
573 for (
int i=1; i<=3; i++) {
579 for (
int i=1; i<=3; i++) {
580 tmps2.set(i) = 2. * psibh_iso
581 * (dlapbh_iso + 0.5*(lapbh_iso - 1.)
586 tmps2.std_spectral_base() ;
587 tmps2.annule_domain(0) ;
588 for (
int i=1; i<=3; i++) {
589 tmps2.
set(i).raccord(1) ;
591 for (
int i=1; i<=3; i++) {
592 (tmps2.set(i)).inc_dzpuis(2) ;
597 for (
int i=1; i<=3; i++) {
598 tmps3.set(i) = 2. * cc * mass * mass * lapbh_iso
599 * (dlappsi(i) - 3.*ll(i)*(ll(1)%dlappsi(1)
604 tmps3.std_spectral_base() ;
605 tmps3.annule_domain(0) ;
606 for (
int i=1; i<=3; i++) {
607 tmps3.set(i).raccord(1) ;
609 for (
int i=1; i<=3; i++) {
610 (tmps3.set(i)).inc_dzpuis(2) ;
613 Vector tmps4 = 4. * cc * mass * mass
614 * (dlapbh_iso * (1. - psibh_iso*lapbh_iso/
lapconf_tot)
615 + 0.5 * lapbh_iso * (lapbh_iso - 1.)
618 * ll / rr /
pow(r_are*rr,3.) ;
621 for (
int i=1; i<=3; i++) {
624 for (
int i=1; i<=3; i++) {
625 (tmps4.
set(i)).inc_dzpuis(4) ;
629 for (
int i=1; i<=3; i++) {
642 for (
int i=1; i<=3; i++) {
646 source_shift = tmps1 + tmps2 + tmps3 + tmps4 + tmps5 ;
648 source_shift.annule_domain(0) ;
650 for (
int i=1; i<=3; i++) {
651 source_shift.set(i).raccord(1) ;
654 if (filter_r_s != 0) {
655 for (
int i=1; i<=3; i++) {
656 if (source_shift(i).get_etat() != ETATZERO)
657 source_shift.set(i).filtre(filter_r_s) ;
661 if (filter_p_s != 0) {
662 for (
int i=1; i<=3; i++) {
663 if (source_shift(i).get_etat() != ETATZERO) {
664 source_shift.set(i).filtre_phi(filter_p_s, nz-1) ;
687 for (
int i=0; i<3; i++) {
688 source_p.set(i) =
Cmp(source_shift(i+1)) ;
694 for (
int i=0; i<3; i++) {
703 poisson_vect_frontiere(1./3., source_p, resu_p,
704 bc_shif_x, bc_shif_y, bc_shif_z,
707 for (
int i=1; i<=3; i++) {
713 for (
int i=1; i<=3; i++) {
720 for (
int i=1; i<=3; i++) {
733 tdiff_lapconf.
set(0) = 0. ;
734 cout <<
"Relative difference in the lapconf function : " << endl ;
735 for (
int l=0; l<nz; l++) {
736 cout << tdiff_lapconf(l) <<
" " ;
740 diff_lapconf = tdiff_lapconf(1) ;
741 for (
int l=2; l<nz; l++) {
742 diff_lapconf += tdiff_lapconf(l) ;
747 tdiff_confo.
set(0) = 0. ;
748 cout <<
"Relative difference in the conformal factor : " << endl ;
749 for (
int l=0; l<nz; l++) {
750 cout << tdiff_confo(l) <<
" " ;
754 diff_confo = tdiff_confo(1) ;
755 for (
int l=2; l<nz; l++) {
756 diff_confo += tdiff_confo(l) ;
761 tdiff_shift_x.
set(0) = 0. ;
762 cout <<
"Relative difference in the shift vector (x) : " << endl ;
763 for (
int l=0; l<nz; l++) {
764 cout << tdiff_shift_x(l) <<
" " ;
768 diff_shift_x = tdiff_shift_x(1) ;
769 for (
int l=2; l<nz; l++) {
770 diff_shift_x += tdiff_shift_x(l) ;
775 tdiff_shift_y.
set(0) = 0. ;
776 cout <<
"Relative difference in the shift vector (y) : " << endl ;
777 for (
int l=0; l<nz; l++) {
778 cout << tdiff_shift_y(l) <<
" " ;
782 diff_shift_y = tdiff_shift_y(1) ;
783 for (
int l=2; l<nz; l++) {
784 diff_shift_y += tdiff_shift_y(l) ;
789 tdiff_shift_z.
set(0) = 0. ;
790 cout <<
"Relative difference in the shift vector (z) : " << endl ;
791 for (
int l=0; l<nz; l++) {
792 cout << tdiff_shift_z(l) <<
" " ;
796 diff_shift_z = tdiff_shift_z(1) ;
797 for (
int l=2; l<nz; l++) {
798 diff_shift_z += tdiff_shift_z(l) ;
Map & mp
Mapping associated with the black hole.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Scalar lapconf_auto_bh
Part of the lapconf function from the analytic background.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Scalar taij_quad_comp
Part of the scalar from the companion star.
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Cmp sqrt(const Cmp &)
Square root.
Scalar confo_auto
Conformal factor generated by the black hole.
double mass_bh
Gravitational mass of BH.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Vector shift_auto
Shift vector generated by the black hole.
Sym_tensor taij_comp
Part of the extrinsic curvature tensor generated by the companion star.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Values and coefficients of a (real-value) function.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tensor field of valence 1.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Sym_tensor taij_auto_rs
Part of the extrinsic curvature tensor numericalty computed for the black hole.
Scalar lapconf_comp
Lapconf function generated by the companion star.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
const Valeur bc_shift_x(double ome_orb, double y_rot) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const Valeur bc_shift_y(double ome_orb, double x_rot) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Scalar confo_auto_rs
Part of the conformal factor from the numerical computation.
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.
Vector d_lapconf_comp
Derivative of the lapconf function generated by the companion star.
Vector shift_auto_bh
Part of the shift vector from the analytic background.
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 Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Scalar confo_comp
Conformal factor generated by the companion star.
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
Scalar taij_quad_auto
Part of the scalar from the black hole.
Vector shift_auto_rs
Part of the shift vector from the numerical computation.
int get_nzone() const
Returns the number of domains.
const Scalar & deriv(int i) const
Returns of *this , where .
Cmp pow(const Cmp &, int)
Power .
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Scalar lapconf_auto
Lapconf function generated by the black hole.
Scalar confo_tot
Total conformal factor.
Tbl & set(int l)
Read/write of the value in a given domain.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Vector d_confo_comp
Derivative of the conformal factor generated by the companion star.
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Scalar lapconf_auto_rs
Part of the lapconf function from the numerical computation.
bool bc_lapconf_fs
true for the first type BC for the lapconf function, false for the second type BH ...
bool bc_lapconf_nd
true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH ...
Scalar & set(int)
Read/write access to a component.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
const Valeur bc_lapconf() const
Boundary condition on the apparent horizon of the black hole for the lapconf function: 2-D Valeur...
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur...
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Scalar lapconf_tot
Total lapconf function.
Coord r
r coordinate centered on the grid
void equilibrium_bhns(int mer, int mermax_bh, int filter_r, int filter_r_s, int filter_p_s, double x_rot, double y_rot, double precis, double omega_orb, double resize_bh, const Tbl &fact_resize, Tbl &diff)
Computes a black-hole part in a black hole-neutron star binary by giving boundary conditions on the a...
Scalar confo_auto_bh
Part of the conformal factor from the analytic background.