23char bin_ns_bh_glob_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_glob.C,v 1.8 2014/10/13 08:52:43 j_novak Exp $" ;
72#include "utilitaires.h"
78double Bin_ns_bh::adm_systeme()
const {
85 auxi_mp.integrale_surface_infini(der_deux) ;
91double Bin_ns_bh::adm_systeme_volume()
const {
99 work_bh.set_etat_qcq() ;
100 for (
int i=0 ; i<3 ; i++) {
101 work_bh.set() = auxi_bh(i, i) ;
102 kk_bh = kk_bh + work_bh ;
106 integ_bh.std_base_scal() ;
116 rad.set(0) =
cos(phi)*
sin(tet) ;
117 rad.set(1) =
sin(phi)*
sin(tet) ;
118 rad.set(2) =
cos(tet) ;
120 Cmp integ_hor2 (
hole.
mp) ;
121 integ_hor2.annule_hard() ;
122 integ_hor2.set_dzpuis(2) ;
123 for (
int m=0 ; m<3 ; m++)
124 for (
int n=0 ; n<3 ; n++)
132 Tenseur work_ns(
star.
mp) ;
133 work_ns.set_etat_qcq() ;
134 for (
int i=0 ; i<3 ; i++) {
135 work_ns.set() = auxi_ns(i, i) ;
136 kk_ns = kk_ns + work_ns ;
139 integ_ns.std_base_scal() ;
142 integ_matter.std_base_scal() ;
144 double masse = (integ_bh.integrale()+integ_ns.integrale())/16/M_PI +
147 + integ_matter.integrale()*ggrav ;
152double Bin_ns_bh::komar_systeme()
const {
159 auxi_mp.integrale_surface_infini(der_deux) ;
166double Bin_ns_bh::viriel()
const {
167 double adm = adm_systeme() ;
168 double komar = komar_systeme() ;
170 return (adm-komar)/adm ;
173double Bin_ns_bh::moment_systeme_inf()
const {
181 double* bornes =
new double [nzones+1] ;
183 for (
int i=nzones-1 ; i>0 ; i--) {
184 bornes[i] = courant ;
188 bornes[nzones] = __infinity ;
195 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
196 k_total.set_etat_qcq() ;
198 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
199 shift_un.set_etat_qcq() ;
201 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
202 shift_deux.set_etat_qcq() ;
208 shift_un.change_triad (mapping.get_bvect_cart()) ;
214 shift_deux.change_triad(mapping.get_bvect_cart()) ;
216 Tenseur shift_tot (shift_un+shift_deux) ;
217 shift_tot.set_std_base() ;
218 shift_tot.annule(0, nzones-2) ;
220 shift_tot.inc2_dzpuis() ;
221 shift_tot.dec2_dzpuis() ;
223 Tenseur grad (shift_tot.gradient()) ;
224 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
225 for (
int i=0 ; i<3 ; i++) {
226 k_total.set(i, i) = grad(i, i)-trace()/3. ;
227 for (
int j=i+1 ; j<3 ; j++)
228 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
231 for (
int lig=0 ; lig<3 ; lig++)
232 for (
int col=lig ; col<3 ; col++)
233 k_total.set(lig, col).mult_r_zec() ;
235 Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
236 vecteur_un.set_etat_qcq() ;
237 for (
int i=0 ; i<3 ; i++)
238 vecteur_un.set(i) = k_total(0, i) ;
239 vecteur_un.change_triad (mapping.get_bvect_spher()) ;
240 Cmp integrant_un (vecteur_un(0)) ;
242 Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
243 vecteur_deux.set_etat_qcq() ;
244 for (
int i=0 ; i<3 ; i++)
245 vecteur_deux.set(i) = k_total(1, i) ;
246 vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
247 Cmp integrant_deux (vecteur_deux(0)) ;
250 integrant_un.va = integrant_un.va.mult_st() ;
251 integrant_un.va = integrant_un.va.mult_sp() ;
253 integrant_deux.va = integrant_deux.va.mult_st() ;
254 integrant_deux.va = integrant_deux.va.mult_cp() ;
256 double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
263double Bin_ns_bh::moment_systeme_hor()
const {
273 xa_bh.std_base_scal() ;
277 ya_bh.std_base_scal() ;
280 vecteur_bh.set_etat_qcq() ;
281 for (
int i=0 ; i<3 ; i++)
284 vecteur_bh.set_std_base() ;
289 integrant_bh.std_base_scal() ;
296 xa_ns.std_base_scal() ;
300 ya_ns.std_base_scal() ;
304 integrant_ns.std_base_scal() ;
306 double moment_ns = integrant_ns.integrale() * ggrav ;
307 return moment_ns + moment_bh ;
311Tbl Bin_ns_bh::linear_momentum_systeme_inf()
const {
318 double* bornes =
new double [nzones+1] ;
320 for (
int i=nzones-1 ; i>0 ; i--) {
321 bornes[i] = courant ;
325 bornes[nzones] = __infinity ;
332 Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
333 k_total.set_etat_qcq() ;
335 Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
336 shift_un.set_etat_qcq() ;
338 Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
339 shift_deux.set_etat_qcq() ;
345 shift_un.change_triad (mapping.get_bvect_cart()) ;
351 shift_deux.change_triad(mapping.get_bvect_cart()) ;
353 shift_un.set_std_base() ;
354 shift_deux.set_std_base() ;
356 Tenseur shift_tot (shift_un+shift_deux) ;
357 shift_tot.set_std_base() ;
358 shift_tot.annule(0, nzones-2) ;
360 Cmp compy (shift_tot(1)) ;
363 int nr = mapping.get_mg()->get_nr(nzones-1) ;
364 int nt = mapping.get_mg()->get_nt(nzones-1) ;
365 int np = mapping.get_mg()->get_np(nzones-1) ;
366 Tbl val_inf (nt*np) ;
367 val_inf.set_etat_qcq() ;
368 for (
int k=0 ; k<np ; k++)
369 for (
int j=0 ; j<nt ; j++)
370 val_inf.set(k*nt + j) = fabs(compy (nzones-1, k, j, nr-1)) ;
372 Tenseur grad (shift_tot.gradient()) ;
373 Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
374 for (
int i=0 ; i<3 ; i++) {
375 k_total.set(i, i) = grad(i, i)-trace()/3. ;
376 for (
int j=i+1 ; j<3 ; j++)
377 k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
380 for (
int comp=0 ; comp<3 ; comp++) {
381 Tenseur vecteur (mapping, 1, CON, mapping.get_bvect_cart()) ;
382 vecteur.set_etat_qcq() ;
383 for (
int i=0 ; i<3 ; i++)
384 vecteur.set(i) = k_total(i, comp) ;
385 vecteur.change_triad (mapping.get_bvect_spher()) ;
386 Cmp integrant (vecteur(0)) ;
388 res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
393double Bin_ns_bh::distance_propre_axe_bh (
const int nr)
const {
398 double pente = -2./x_bh ;
399 double constante = - 1. ;
403 double air_un, theta_un, phi_un ;
404 double air_deux, theta_deux, phi_deux ;
406 double* coloc =
new double[nr] ;
407 double* coef =
new double[nr] ;
408 int* deg =
new int[3] ;
409 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
411 for (
int i=0 ; i<nr ; i++) {
412 ksi = -
cos (M_PI*i/(nr-1)) ;
413 xabs = (ksi+constante)/pente ;
423 cfrcheb(deg, deg, coloc, deg, coef) ;
426 double* som =
new double[nr] ;
428 for (
int i=2 ; i<nr ; i+=2)
429 som[i] = 1./(i+1)-1./(i-1) ;
430 for (
int i=1 ; i<nr ; i+=2)
434 for (
int i=0 ; i<nr ; i++)
435 res += som[i]*coef[i] ;
447double Bin_ns_bh::distance_propre_axe_ns (
const int nr)
const {
452 double pente = 2./x_ns ;
453 double constante = - 1. ;
457 double air_un, theta_un, phi_un ;
458 double air_deux, theta_deux, phi_deux ;
460 double* coloc =
new double[nr] ;
461 double* coef =
new double[nr] ;
462 int* deg =
new int[3] ;
463 deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
465 for (
int i=0 ; i<nr ; i++) {
466 ksi = -
cos (M_PI*i/(nr-1)) ;
467 xabs = (ksi-constante)/pente ;
477 cfrcheb(deg, deg, coloc, deg, coef) ;
480 double* som =
new double[nr] ;
482 for (
int i=2 ; i<nr ; i+=2)
483 som[i] = 1./(i+1)-1./(i-1) ;
484 for (
int i=1 ; i<nr ; i+=2)
488 for (
int i=0 ; i<nr ; i++)
489 res += som[i]*coef[i] ;
501double Bin_ns_bh::smarr()
const {
507 psiq_t.set_std_base() ;
510 furmet.set_etat_qcq() ;
511 for (
int i=0 ; i< 3 ; i++) {
512 furmet.set(i,i) = 1/psiq_t() ;
513 for(
int j=i+1 ; j<3 ; j++)
514 furmet.set(i,j) = 0 ;
516 Metrique met (furmet,
false) ;
521 Tenseur_sym kij_cov (
manipule (kij, met)) ;
526 aime.set_etat_qcq() ;
531 aime.set_std_base() ;
532 shift = shift - aime ;
537 Tenseur u_i_bas (
manipule(u_euler, met)) ;
546 integ_matter.std_base_scal() ;
548 double matter_term = integ_matter.integrale()*qpig/4/M_PI ;
557 rad.set(0) =
cos(phi)*
sin(tet) ;
558 rad.set(1) =
sin(phi)*
sin(tet) ;
559 rad.set(2) =
cos(tet) ;
564 for (
int m=0 ; m<3 ; m++)
565 for (
int n=0 ; n<3 ; n++)
568 temp.std_base_scal() ;
572 integ_hor.std_base_scal() ;
573 integ_hor.raccord(1) ;
577 double m_test = hor_term + matter_term + 2*
omega*moment_systeme_inf() +
Tenseur psi_auto
Part of generated by the hole.
Tenseur shift_auto
Part of generated by the hole.
double omega_local
local angular velocity
Tenseur_sym tkij_auto
Auto .
Tenseur n_auto
Part of N generated by the hole.
double rayon
Radius of the horizon in LORENE's units.
Map_af & mp
Affine mapping.
const Tenseur & get_psi_tot() const
Returns the total .
Tenseur_sym tkij_tot
Total .
const Tenseur & get_n_tot() const
Returns the total N .
double get_rayon() const
Returns the radius of the horizon.
Et_bin_nsbh star
The neutron star.
Bhole hole
The black hole.
double omega
Angular velocity with respect to an asymptotically inertial observer.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Tenseur confpsi_auto
Part of the conformal factor $\Psi$ generated principaly by the star.
Tenseur confpsi_comp
Part of the conformal factor $\Psi$ generated principaly by the companion star.
Tenseur_sym tkij_tot
Total extrinsic curvature tensor $K^{ij}$ generated by {\tt shift_auto} and {\tt shift_comp}.
const Tenseur_sym & get_tkij_tot() const
Returns the total extrinsic curvature tensor $K^{ij}$ generated by {\tt shift_auto} and {\tt shift_co...
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor $K^{ij}$ generated by {\tt shift_auto}.
Tenseur n_auto
Part of the lapse {\it N} generated principaly by the star.
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $\Psi$.
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
int nzet
Number of domains of *mp occupied by the star.
const Tenseur & get_nnn() const
Returns the total lapse function N.
int get_nzet() const
Returns the number of domains occupied by the star.
const Tenseur & get_shift() const
Returns the total shift vector .
const Map & get_mp() const
Returns the mapping.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
Tenseur press
Fluid pressure.
const Tenseur & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
const Tenseur & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Tenseur ener_euler
Total energy density in the Eulerian frame.
const Tenseur & get_press() const
Returns the fluid pressure.
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
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 .
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord ya
Absolute y coordinate.
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,...
Coord tet
coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
double get_ori_x() const
Returns the x coordinate of the origin.
Coord phi
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_nzone() const
Returns the number of domains.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Cmp sin(const Cmp &)
Sine.
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
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 manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Standard units of space, time and mass.