26 char Binary_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Binary/binary.C,v 1.16 2014/10/13 08:52:44 j_novak Exp $" ;
83 #include "utilitaires.h" 84 #include "graphique.h" 97 Map& mp2,
int nzet2,
const Eos& eos2,
int irrot2,
99 : star1(mp1, nzet1, eos1, irrot1, conf_flat),
100 star2(mp2, nzet2, eos2, irrot2, conf_flat)
132 :
star1(mp1, eos1, fich),
133 star2(mp2, eos2, fich)
238 ost <<
"Binary neutron stars" << endl ;
239 ost <<
"=============" << endl ;
241 "Orbital angular velocity : " <<
omega * f_unit <<
" rad/s" << endl ;
243 "Coordinate separation between the two stellar centers : " 246 "Absolute coordinate X of the rotation axis : " <<
x_axe / km
248 ost << endl <<
"Star 1 : " << endl ;
249 ost <<
"====== " << endl ;
250 ost <<
star1 << endl ;
251 ost <<
"Star 2 : " << endl ;
252 ost <<
"====== " << endl ;
253 ost <<
star2 << endl ;
267 if (p_eos_poly != 0x0) {
271 double kappa = p_eos_poly->
get_kap() ;
272 double gamma = p_eos_poly->
get_gam() ; ;
273 double kap_ns2 =
pow( kappa, 0.5 /(gamma-1) ) ;
276 double r_poly = kap_ns2 /
sqrt(ggrav) ;
279 double t_poly = r_poly ;
282 double m_poly = r_poly / ggrav ;
285 double j_poly = r_poly * r_poly / ggrav ;
288 ost << endl <<
"Quantities in polytropic units : " << endl ;
289 ost <<
"==============================" << endl ;
290 ost <<
" ( r_poly = " << r_poly / km <<
" km )" << endl ;
291 ost <<
" d_e_max : " <<
separation() / r_poly << endl ;
295 ost <<
" Omega : " <<
omega * t_poly << endl ;
296 ost <<
" J : " <<
angu_mom()(2) / j_poly << endl ;
297 ost <<
" M_ADM : " <<
mass_adm() / m_poly << endl ;
298 ost <<
" M_Komar : " <<
mass_kom() / m_poly << endl ;
299 ost <<
" M_bar(star 1) : " <<
star1.
mass_b() / m_poly << endl ;
300 ost <<
" M_bar(star 2) : " <<
star2.
mass_b() / m_poly << endl ;
301 ost <<
" R_0(star 1) : " <<
303 ost <<
" R_0(star 2) : " <<
319 cout <<
"distance = " << distance << endl ;
320 double lim_un = distance/2. ;
321 double lim_deux = distance/2. ;
322 double int_un = distance/6. ;
323 double int_deux = distance/6. ;
327 fonction_f_un = 0.5*
pow(
329 fonction_f_un.std_spectral_base();
332 fonction_g_un = 0.5*
pow 334 fonction_g_un.std_spectral_base();
337 fonction_f_deux = 0.5*
pow(
338 cos((
star2.
get_mp().
r-int_deux)*M_PI/2./(lim_deux-int_deux)), 2.)+0.5 ;
339 fonction_f_deux.std_spectral_base();
342 fonction_g_deux = 0.5*
pow(
344 fonction_g_deux.std_spectral_base();
360 double xabs, yabs, zabs, air_un, air_deux, theta, phi ;
362 for (
int l=0 ; l<nz_un ; l++) {
371 for (
int k=0 ; k<np ; k++)
372 for (
int j=0 ; j<nt ; j++)
373 for (
int i=0 ; i<nr ; i++) {
375 xabs = xabs_un (l, k, j, i) ;
376 yabs = yabs_un (l, k, j, i) ;
377 zabs = zabs_un (l, k, j, i) ;
381 (xabs, yabs, zabs, air_un, theta, phi) ;
383 (xabs, yabs, zabs, air_deux, theta, phi) ;
385 if (air_un <= lim_un)
387 decouple_un.set_grid_point(l, k, j, i) = 1 ;
390 decouple_un.set_grid_point(l, k, j, i) =
391 fonction_f_un.val_grid_point(l, k, j, i) ;
393 if (air_deux <= lim_deux)
394 if (air_deux < int_deux)
395 decouple_un.set_grid_point(l, k, j, i) = 0 ;
398 decouple_un.set_grid_point(l, k, j, i) =
399 fonction_g_deux.val_point (air_deux, theta, phi) ;
403 decouple_un.set_grid_point(l, k, j, i) = 0.5 ;
408 for (
int k=0 ; k<np ; k++)
409 for (
int j=0 ; j<nt ; j++)
410 decouple_un.set_grid_point(nz_un-1, k, j, nr)=0.5 ;
413 for (
int l=0 ; l<nz_deux ; l++) {
422 for (
int k=0 ; k<np ; k++)
423 for (
int j=0 ; j<nt ; j++)
424 for (
int i=0 ; i<nr ; i++) {
426 xabs = xabs_deux (l, k, j, i) ;
427 yabs = yabs_deux (l, k, j, i) ;
428 zabs = zabs_deux (l, k, j, i) ;
432 (xabs, yabs, zabs, air_un, theta, phi) ;
434 (xabs, yabs, zabs, air_deux, theta, phi) ;
436 if (air_deux <= lim_deux)
437 if (air_deux < int_deux)
438 decouple_deux.set_grid_point(l, k, j, i) = 1 ;
441 decouple_deux.set_grid_point(l, k, j, i) =
442 fonction_f_deux.val_grid_point(l, k, j, i) ;
444 if (air_un <= lim_un)
446 decouple_deux.set_grid_point(l, k, j, i) = 0 ;
449 decouple_deux.set_grid_point(l, k, j, i) =
450 fonction_g_un.val_point (air_un, theta, phi) ;
454 decouple_deux.set_grid_point(l, k, j, i) = 0.5 ;
459 for (
int k=0 ; k<np ; k++)
460 for (
int j=0 ; j<nt ; j++)
461 decouple_deux.set_grid_point(nz_un-1, k, j, nr)=0.5 ;
464 decouple_un.std_spectral_base() ;
465 decouple_deux.std_spectral_base() ;
470 cout <<
"decouple_un" << endl <<
norme(decouple_un/(nr*nt*np)) << endl ;
471 cout <<
"decouple_deux" << endl <<
norme(decouple_deux/(nr*nt*np))
487 ost <<
"# Grid 1 : " << nz1 <<
"x" 489 <<
" R_out(l) [km] : " ;
490 for (
int l=0; l<nz1; l++) {
491 ost <<
" " << mp1.
val_r(l, 1., M_PI/2, 0) / km ;
495 ost <<
"# VE(M) " << endl ;
498 ost.setf(ios::scientific) ;
507 <<
" M_ADM_vol [M_sol] " 508 <<
" M_Komar [M_sol] " 509 <<
" M_Komar_vol [M_sol] " 510 <<
" J [G M_sol^2/c] " << endl ;
517 ost <<
omega / (2*M_PI)* f_unit ; ost.width(22) ;
518 ost <<
mass_adm() / msol ; ost.width(22) ;
520 ost <<
mass_kom() / msol ; ost.width(22) ;
522 ost <<
angu_mom()(2)/ ( qpig / (4* M_PI) * msol*msol) << endl ;
524 ost <<
"# H_c(1)[c^2] " 525 <<
" e_c(1)[rho_nuc] " 526 <<
" M_B(1) [M_sol] " 529 <<
" a3/a1(1) " << endl ;
539 ost <<
"# H_c(2)[c^2] " 540 <<
" e_c(2)[rho_nuc] " 541 <<
" M_B(2) [M_sol] " 544 <<
" a3/a1(2) " << endl ;
561 double kappa = p_eos_poly->
get_kap() ;
562 double gamma = p_eos_poly->
get_gam() ; ;
563 double kap_ns2 =
pow( kappa, 0.5 /(gamma-1.) ) ;
566 double r_poly = kap_ns2 /
sqrt(ggrav) ;
569 double t_poly = r_poly ;
572 double m_poly = r_poly / ggrav ;
575 double j_poly = r_poly * r_poly / ggrav ;
583 <<
" M_B(2) [poly] " << endl ;
586 ost <<
separation() / r_poly ; ost.width(22) ;
588 ost <<
omega * t_poly ; ost.width(22) ;
589 ost <<
mass_adm() / m_poly ; ost.width(22) ;
590 ost <<
angu_mom()(2) / j_poly ; ost.width(22) ;
610 return sqrt( dx*dx + dy*dy + dz*dz ) ;
Coord xa
Absolute x coordinate.
double mass_kom_vol() const
Total Komar mass (computed by a volume integral)
Map & mp
Mapping associated with the star.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
virtual double mass_b() const
Baryon mass.
Cmp sqrt(const Cmp &)
Square root.
double get_ori_y() const
Returns the y coordinate of the origin.
Star_bin star1
First star of the system.
Standard units of space, time and mass.
Equation of state base class.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Tensor field of valence 0 (or component of a tensorial field).
const Scalar & get_ener() const
Returns the proper total energy density.
ostream & operator>>(ostream &) const
Operator >> (function called by the operator <<).
Base class for coordinate mappings.
double get_ori_x() const
Returns the x coordinate of the origin.
const Map & get_mp() const
Returns the mapping.
void fait_decouple()
Calculates decouple which is used to obtain qq_auto by the formula : qq_auto = decouple * qq...
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
double mass_adm_vol() const
Total ADM mass (computed by a volume integral)
Cmp cos(const Cmp &)
Cosine.
double * p_mass_adm
Total ADM mass of the system.
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
double virial() const
Estimates the relative error on the virial theorem.
void display_poly(ostream &) const
Display in polytropic units.
void sauve(FILE *) const
Save in a file.
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.
double * p_ham_constr
Relative error on the Hamiltonian constraint.
Scalar decouple
Function used to construct the part generated by the star from the total .
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
const Tbl & angu_mom() const
Total angular momentum.
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
double get_kap() const
Returns the pressure coefficient (cf.
Tbl * p_mom_constr
Relative error on the momentum constraint.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
void operator=(const Binary &)
Assignment to another Binary.
double * p_total_ener
Total energy of the system.
Polytropic equation of state (relativistic case).
void del_deriv() const
Deletes all the derived quantities.
int get_nzone() const
Returns the number of domains.
double * p_virial
Virial theorem error.
friend ostream & operator<<(ostream &, const Binary &)
Display.
double ray_eq_pi() const
Coordinate radius at , [r_unit].
virtual void sauve(FILE *) const
Save in a file.
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
double x_axe
Absolute X coordinate of the rotation axis.
Cmp pow(const Cmp &, int)
Power .
double ray_pole() const
Coordinate radius at [r_unit].
Star_bin * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Coord ya
Absolute y coordinate.
void write_global(ostream &) const
Write global quantities in a formatted file.
Tbl * p_angu_mom
Total angular momentum of the system.
const Scalar & get_ent() const
Returns the enthalpy field.
Coord za
Absolute z coordinate.
double mass_kom() const
Total Komar mass.
double ray_eq() const
Coordinate radius at , [r_unit].
double get_ori_z() const
Returns the z coordinate of the origin.
Binary(Map &mp1, int nzet1, const Eos &eos1, int irrot1, Map &mp2, int nzet2, const Eos &eos2, int irrot2, int conf_flat)
Standard constructor.
const Eos & get_eos() const
Returns the equation of state.
Cmp sin(const Cmp &)
Sine.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
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...
double mass_adm() const
Total ADM mass.
Star_bin star2
Second star of the system.
double * p_mass_kom
Total Komar mass of the system.
Coord r
r coordinate centered on the grid
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.