73 #include "bin_ns_bh.h" 75 #include "utilitaires.h" 76 #include "graphique.h" 81 double Bin_ns_bh::adm_systeme()
const {
88 auxi_mp.integrale_surface_infini(der_deux) ;
94 double Bin_ns_bh::adm_systeme_volume()
const {
101 Tenseur work_bh(
hole.
mp) ;
102 work_bh.set_etat_qcq() ;
103 for (
int i=0 ; i<3 ; i++) {
104 work_bh.set() = auxi_bh(i, i) ;
105 kk_bh = kk_bh + work_bh ;
109 integ_bh.std_base_scal() ;
120 rad.set(1) =
sin(phi)*
sin(tet) ;
121 rad.set(2) =
cos(tet) ;
123 Cmp integ_hor2 (
hole.
mp) ;
124 integ_hor2.annule_hard() ;
125 integ_hor2.set_dzpuis(2) ;
126 for (
int m=0 ; m<3 ; m++)
127 for (
int n=0 ; n<3 ; n++)
135 Tenseur work_ns(
star.
mp) ;
136 work_ns.set_etat_qcq() ;
137 for (
int i=0 ; i<3 ; i++) {
138 work_ns.set() = auxi_ns(i, i) ;
139 kk_ns = kk_ns + work_ns ;
142 integ_ns.std_base_scal() ;
145 integ_matter.std_base_scal() ;
147 double masse = (integ_bh.integrale()+integ_ns.integrale())/16/M_PI +
150 + integ_matter.integrale()*ggrav ;
155 double Bin_ns_bh::komar_systeme()
const {
162 auxi_mp.integrale_surface_infini(der_deux) ;
169 double Bin_ns_bh::viriel()
const {
170 double adm = adm_systeme() ;
171 double komar = komar_systeme() ;
173 return (adm-komar)/adm ;
176 double Bin_ns_bh::moment_systeme_inf()
const {
184 double* bornes =
new double [nzones+1] ;
186 for (
int i=nzones-1 ; i>0 ; i--) {
187 bornes[i] = courant ;
191 bornes[nzones] = __infinity ;
198 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
199 k_total.set_etat_qcq() ;
201 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
202 shift_un.set_etat_qcq() ;
204 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
205 shift_deux.set_etat_qcq() ;
211 shift_un.change_triad (mapping.get_bvect_cart()) ;
217 shift_deux.change_triad(mapping.get_bvect_cart()) ;
219 Tenseur shift_tot (shift_un+shift_deux) ;
220 shift_tot.set_std_base() ;
221 shift_tot.annule(0, nzones-2) ;
223 shift_tot.inc2_dzpuis() ;
224 shift_tot.dec2_dzpuis() ;
226 Tenseur grad (shift_tot.gradient()) ;
227 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
228 for (
int i=0 ; i<3 ; i++) {
229 k_total.set(i, i) = grad(i, i)-trace()/3. ;
230 for (
int j=i+1 ; j<3 ; j++)
231 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
234 for (
int lig=0 ; lig<3 ; lig++)
235 for (
int col=lig ; col<3 ; col++)
236 k_total.set(lig, col).mult_r_zec() ;
238 Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
239 vecteur_un.set_etat_qcq() ;
240 for (
int i=0 ; i<3 ; i++)
241 vecteur_un.set(i) = k_total(0, i) ;
242 vecteur_un.change_triad (mapping.get_bvect_spher()) ;
243 Cmp integrant_un (vecteur_un(0)) ;
245 Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
246 vecteur_deux.set_etat_qcq() ;
247 for (
int i=0 ; i<3 ; i++)
248 vecteur_deux.set(i) = k_total(1, i) ;
249 vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
250 Cmp integrant_deux (vecteur_deux(0)) ;
253 integrant_un.va = integrant_un.va.mult_st() ;
254 integrant_un.va = integrant_un.va.mult_sp() ;
256 integrant_deux.va = integrant_deux.va.mult_st() ;
257 integrant_deux.va = integrant_deux.va.mult_cp() ;
259 double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
266 double Bin_ns_bh::moment_systeme_hor()
const {
276 xa_bh.std_base_scal() ;
280 ya_bh.std_base_scal() ;
283 vecteur_bh.set_etat_qcq() ;
284 for (
int i=0 ; i<3 ; i++)
287 vecteur_bh.set_std_base() ;
292 integrant_bh.std_base_scal() ;
299 xa_ns.std_base_scal() ;
303 ya_ns.std_base_scal() ;
307 integrant_ns.std_base_scal() ;
309 double moment_ns = integrant_ns.integrale() * ggrav ;
310 return moment_ns + moment_bh ;
314 Tbl Bin_ns_bh::linear_momentum_systeme_inf()
const {
321 double* bornes =
new double [nzones+1] ;
323 for (
int i=nzones-1 ; i>0 ; i--) {
324 bornes[i] = courant ;
328 bornes[nzones] = __infinity ;
335 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
336 k_total.set_etat_qcq() ;
338 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
339 shift_un.set_etat_qcq() ;
341 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
342 shift_deux.set_etat_qcq() ;
348 shift_un.change_triad (mapping.get_bvect_cart()) ;
354 shift_deux.change_triad(mapping.get_bvect_cart()) ;
356 shift_un.set_std_base() ;
357 shift_deux.set_std_base() ;
359 Tenseur shift_tot (shift_un+shift_deux) ;
360 shift_tot.set_std_base() ;
361 shift_tot.annule(0, nzones-2) ;
363 Cmp compy (shift_tot(1)) ;
366 int nr = mapping.get_mg()->get_nr(nzones-1) ;
367 int nt = mapping.get_mg()->get_nt(nzones-1) ;
368 int np = mapping.get_mg()->get_np(nzones-1) ;
369 Tbl val_inf (nt*np) ;
370 val_inf.set_etat_qcq() ;
371 for (
int k=0 ; k<np ; k++)
372 for (
int j=0 ; j<nt ; j++)
373 val_inf.set(k*nt + j) = fabs(compy (nzones-1, k, j, nr-1)) ;
375 Tenseur grad (shift_tot.gradient()) ;
376 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
377 for (
int i=0 ; i<3 ; i++) {
378 k_total.set(i, i) = grad(i, i)-trace()/3. ;
379 for (
int j=i+1 ; j<3 ; j++)
380 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
383 for (
int comp=0 ; comp<3 ; comp++) {
384 Tenseur vecteur (mapping, 1, CON, mapping.get_bvect_cart()) ;
385 vecteur.set_etat_qcq() ;
386 for (
int i=0 ; i<3 ; i++)
387 vecteur.set(i) = k_total(i, comp) ;
388 vecteur.change_triad (mapping.get_bvect_spher()) ;
389 Cmp integrant (vecteur(0)) ;
391 res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
396 double Bin_ns_bh::distance_propre_axe_bh (
const int nr)
const {
401 double pente = -2./x_bh ;
402 double constante = - 1. ;
406 double air_un, theta_un, phi_un ;
407 double air_deux, theta_deux, phi_deux ;
409 double* coloc =
new double[nr] ;
410 double* coef =
new double[nr] ;
411 int* deg =
new int[3] ;
412 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
414 for (
int i=0 ; i<nr ; i++) {
415 ksi = -
cos (M_PI*i/(nr-1)) ;
416 xabs = (ksi+constante)/pente ;
426 cfrcheb(deg, deg, coloc, deg, coef) ;
429 double* som =
new double[nr] ;
431 for (
int i=2 ; i<nr ; i+=2)
432 som[i] = 1./(i+1)-1./(i-1) ;
433 for (
int i=1 ; i<nr ; i+=2)
437 for (
int i=0 ; i<nr ; i++)
438 res += som[i]*coef[i] ;
450 double Bin_ns_bh::distance_propre_axe_ns (
const int nr)
const {
455 double pente = 2./x_ns ;
456 double constante = - 1. ;
460 double air_un, theta_un, phi_un ;
461 double air_deux, theta_deux, phi_deux ;
463 double* coloc =
new double[nr] ;
464 double* coef =
new double[nr] ;
465 int* deg =
new int[3] ;
466 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
468 for (
int i=0 ; i<nr ; i++) {
469 ksi = -
cos (M_PI*i/(nr-1)) ;
470 xabs = (ksi-constante)/pente ;
480 cfrcheb(deg, deg, coloc, deg, coef) ;
483 double* som =
new double[nr] ;
485 for (
int i=2 ; i<nr ; i+=2)
486 som[i] = 1./(i+1)-1./(i-1) ;
487 for (
int i=1 ; i<nr ; i+=2)
491 for (
int i=0 ; i<nr ; i++)
492 res += som[i]*coef[i] ;
504 double Bin_ns_bh::smarr()
const {
510 psiq_t.set_std_base() ;
513 furmet.set_etat_qcq() ;
514 for (
int i=0 ; i< 3 ; i++) {
515 furmet.set(i,i) = 1/psiq_t() ;
516 for(
int j=i+1 ; j<3 ; j++)
517 furmet.set(i,j) = 0 ;
519 Metrique met (furmet,
false) ;
524 Tenseur_sym kij_cov (
manipule (kij, met)) ;
529 aime.set_etat_qcq() ;
534 aime.set_std_base() ;
535 shift = shift - aime ;
540 Tenseur u_i_bas (
manipule(u_euler, met)) ;
549 integ_matter.std_base_scal() ;
551 double matter_term = integ_matter.integrale()*qpig/4/M_PI ;
561 rad.set(1) =
sin(phi)*
sin(tet) ;
562 rad.set(2) =
cos(tet) ;
567 for (
int m=0 ; m<3 ; m++)
568 for (
int n=0 ; n<3 ; n++)
571 temp.std_base_scal() ;
575 integ_hor.std_base_scal() ;
576 integ_hor.raccord(1) ;
580 double m_test = hor_term + matter_term + 2*
omega*moment_systeme_inf() +
Coord xa
Absolute x coordinate.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Map & get_mp() const
Returns the mapping.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Tenseur & get_shift() const
Returns the total shift vector .
double get_ori_x() const
Returns the x coordinate of the origin.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Tenseur confpsi_comp
Part of the conformal factor $$ generated principaly by the companion star.
const Tenseur & get_psi_tot() const
Returns the total .
int get_nzet() const
Returns the number of domains occupied by the star.
double omega_local
local angular velocity
Tenseur psi_auto
Part of generated by the hole.
const Tenseur & get_n_tot() const
Returns the total N .
Tenseur press
Fluid pressure.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Tenseur n_auto
Part of the lapse { N} generated principaly by the star.
Cmp cos(const Cmp &)
Cosine.
Coord tet
coordinate centered on the grid
Tenseur shift_auto
Part of generated by the hole.
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Coord phi
coordinate centered on the grid
const Tenseur & get_press() const
Returns the fluid pressure.
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.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Et_bin_nsbh star
The neutron star.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Map_af & mp
Affine mapping.
const Tenseur & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Map & mp
Mapping associated with the star.
int get_nzone() const
Returns the number of domains.
double get_rayon() const
Returns the radius of the horizon.
const Tenseur & get_nnn() const
Returns the total lapse function N.
double omega
Angular velocity with respect to an asymptotically inertial observer.
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $$.
double rayon
Radius of the horizon in LORENE's units.
Cmp pow(const Cmp &, int)
Power .
Tenseur_sym tkij_auto
Auto .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int nzet
Number of domains of *mp occupied by the star.
Tenseur_sym tkij_tot
Total .
Coord ya
Absolute y coordinate.
Tenseur confpsi_auto
Part of the conformal factor $$ generated principaly by the star.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
const Tenseur_sym & get_tkij_tot() const
Returns the total extrinsic curvature tensor $K^{ij}$ generated by { shift_auto} and { shift_comp}...
Tenseur_sym tkij_tot
Total extrinsic curvature tensor $K^{ij}$ generated by { shift_auto} and { shift_comp}.
Bhole hole
The black hole.
Cmp sin(const Cmp &)
Sine.
Tenseur n_auto
Part of N generated by the hole.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor $K^{ij}$ generated by { shift_auto}.
const Tenseur & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X...
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .