85 #include "utilitaires.h" 98 :
Star(mpi, nzet_i, eos_i),
110 w_shift(mpi, CON, mp.get_bvect_cart()),
112 tkij(mpi, COV, mp.get_bvect_cart()),
119 ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
169 relativistic(et.relativistic),
183 khi_shift(et.khi_shift),
186 ssjm1_nuf(et.ssjm1_nuf),
187 ssjm1_nuq(et.ssjm1_nuq),
188 ssjm1_dzeta(et.ssjm1_dzeta),
189 ssjm1_tggg(et.ssjm1_tggg),
190 ssjm1_khi(et.ssjm1_khi),
191 ssjm1_wshift(et.ssjm1_wshift)
201 :
Star(mpi, eos_i, fich),
212 w_shift(mpi, CON, mp.get_bvect_cart()),
214 tkij(mpi, COV, mp.get_bvect_cart()),
221 ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
228 size_t nread = fread(&
relativistic,
sizeof(
bool), 1, fich) ;
230 cerr <<
"Star_rot::Star_rot(FILE*): Problem reading the data file." << endl ;
455 if (
omega != __infinity) {
456 ost <<
"Uniformly rotating star" << endl ;
457 ost <<
"-----------------------" << endl ;
459 double freq =
omega / (2.*M_PI) ;
460 ost <<
"Omega : " <<
omega * f_unit
461 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
462 ost <<
"Rotation period : " << 1000. / (freq * f_unit) <<
" ms" 467 ost <<
"Differentially rotating star" << endl ;
468 ost <<
"----------------------------" << endl ;
470 double freq = omega_c / (2.*M_PI) ;
471 ost <<
"Central value of Omega : " << omega_c * f_unit
472 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
473 ost <<
"Central rotation period : " << 1000. / (freq * f_unit) <<
" ms" 477 ost <<
"Relativistic star" << endl ;
480 ost <<
"Newtonian star" << endl ;
483 ost <<
"Compactness G M_g /(c^2 R_circ) : " << compact << endl ;
486 if ( (omega_c==0) && (nphi_c==0) ) {
487 ost <<
"Central N^phi : " << nphi_c << endl ;
490 ost <<
"Central N^phi/Omega : " << nphi_c / omega_c << endl ;
493 ost <<
"Error on the virial identity GRV2 : " <<
grv2() << endl ;
494 double xgrv3 =
grv3(&ost) ;
495 ost <<
"Error on the virial identity GRV3 : " << xgrv3 << endl ;
497 double mom_quad_38si =
mom_quad() * rho_unit * (
pow(r_unit,
double(5.))
499 ost <<
"Quadrupole moment Q : " << mom_quad_38si <<
" 10^38 kg m^2" 501 ost <<
"Q / (M R_circ^2) : " 503 ost <<
"c^4 Q / (G^2 M^3) : " 507 ost <<
"Angular momentum J : " 508 <<
angu_mom()/( qpig / (4* M_PI) * msol*msol) <<
" G M_sol^2 / c" 510 ost <<
"c J / (G M^2) : " 513 if (
omega != __infinity) {
515 double mom_iner_38si = mom_iner * rho_unit * (
pow(r_unit,
double(5.))
517 ost <<
"Moment of inertia: " << mom_iner_38si <<
" 10^38 kg m^2" 521 ost <<
"Ratio T/W : " <<
tsw() << endl ;
522 ost <<
"Circumferential equatorial radius R_circ : " 523 <<
r_circ()/km <<
" km" << endl ;
525 ost <<
"Circumferential polar (meridional) radius R_circ_merid : " 527 ost <<
"Flattening R_circ_merid / R_circ_eq : " 529 ost <<
"Surface area : " <<
area()/(km*km) <<
" km^2" << endl ;
530 ost <<
"Mean radius : " <<
mean_radius()/km <<
" km" << endl ;
533 "Skipping polar radius / surface statements due to number of points in phi direction np > 1" 536 ost <<
"Coordinate equatorial radius r_eq : " <<
ray_eq()/km <<
" km" 538 ost <<
"Flattening r_pole/r_eq : " <<
aplat() << endl ;
539 int lsurf =
nzet - 1;
542 ost <<
"Equatorial value of the velocity U: " 544 ost <<
"Redshift at the equator (forward) : " <<
z_eqf() << endl ;
545 ost <<
"Redshift at the equator (backward): " <<
z_eqb() << endl ;
546 ost <<
"Redshift at the pole : " <<
z_pole() << endl ;
549 ost <<
"Central value of log(N) : " 552 ost <<
"Central value of dzeta=log(AN) : " 555 if ( (omega_c==0) && (nphi_c==0) ) {
556 ost <<
"Central N^phi : " << nphi_c << endl ;
559 ost <<
"Central N^phi/Omega : " << nphi_c / omega_c << endl ;
562 ost <<
" ... w_shift (NB: components in the star Cartesian frame) [c] : " 564 <<
w_shift(1).val_grid_point(0, 0, 0, 0) <<
" " 565 <<
w_shift(2).val_grid_point(0, 0, 0, 0) <<
" " 566 <<
w_shift(3).val_grid_point(0, 0, 0, 0) << endl ;
568 ost <<
"Central value of khi_shift [km c] : " 575 "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : " 578 ost << diff_a_b(l) <<
" " ;
585 double r_grav = ggrav *
mass_g() ;
586 double r_isco_appr = 6. * r_grav * ( 1. -
pow(2./3.,1.5) * jdimless ) ;
591 double f_isco_appr = ( 1. + 11. /6. /
sqrt(6.) * jdimless ) / r_grav /
592 (12. * M_PI ) /
sqrt(6.) ;
594 ost << endl <<
"Innermost stable circular orbit (ISCO) : " << endl ;
595 double xr_isco =
r_isco(&ost) ;
596 ost <<
" circumferential radius r_isco = " 597 << xr_isco / km <<
" km" << endl ;
598 ost <<
" (approx. 6M + 1st order in j : " 599 << r_isco_appr / km <<
" km)" << endl ;
600 ost <<
" (approx. 6M : " 601 << 6. * r_grav / km <<
" km)" << endl ;
602 ost <<
" orbital frequency f_isco = " 603 <<
f_isco() * f_unit <<
" Hz" << endl ;
604 ost <<
" (approx. 1st order in j : " 605 << f_isco_appr * f_unit <<
" Hz)" << endl ;
618 double freq = omega_c / (2.*M_PI) ;
619 ost <<
"Central Omega : " << omega_c * f_unit
620 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
621 ost <<
"Rotation period : " << 1000. / (freq * f_unit) <<
" ms" 623 ost << endl <<
"Central enthalpy : " <<
ent.
val_grid_point(0,0,0,0) <<
" c^2" << endl ;
625 <<
" x 0.1 fm^-3" << endl ;
627 <<
" rho_nuc c^2" << endl ;
629 <<
" rho_nuc c^2" << endl ;
631 ost <<
"Central value of log(N) : " 634 ost <<
"Central value of dzeta=log(AN) : " 640 if ( (omega_c==0) && (nphi_c==0) ) {
641 ost <<
"Central N^phi : " << nphi_c << endl ;
644 ost <<
"Central N^phi/Omega : " << nphi_c / omega_c
649 int lsurf =
nzet - 1;
652 ost <<
"Equatorial value of the velocity U: " 656 <<
"Coordinate equatorial radius r_eq = " 657 <<
ray_eq()/km <<
" km" << endl ;
658 ost <<
"Flattening r_pole/r_eq : " <<
aplat() << endl ;
660 ost <<
"Flattening r_circ_pole/r_circ_eq : " 681 if (p_eos_poly != 0x0) {
683 double kappa = p_eos_poly->
get_kap() ;
686 double kap_ns2 =
pow( kappa, 0.5 /(p_eos_poly->
get_gam()-1) ) ;
689 double r_poly = kap_ns2 /
sqrt(ggrav) ;
692 double t_poly = r_poly ;
695 double m_poly = r_poly / ggrav ;
698 double j_poly = r_poly * r_poly / ggrav ;
701 double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
704 ost << endl <<
"Quantities in polytropic units : " << endl ;
705 ost <<
"==============================" << endl ;
706 ost <<
" ( r_poly = " << r_poly / km <<
" km )" << endl ;
709 ost <<
" Omega_c : " <<
get_omega_c() * t_poly << endl ;
710 ost <<
" P_c : " << 2.*M_PI /
get_omega_c() / t_poly << endl ;
711 ost <<
" M_bar : " <<
mass_b() / m_poly << endl ;
712 ost <<
" M : " <<
mass_g() / m_poly << endl ;
713 ost <<
" J : " <<
angu_mom() / j_poly << endl ;
714 ost <<
" r_eq : " <<
ray_eq() / r_poly << endl ;
715 ost <<
" R_circ_eq : " <<
r_circ() / r_poly << endl ;
716 ost <<
" R_circ_merid : " <<
r_circ_merid() / r_poly << endl ;
717 ost <<
" R_mean : " <<
mean_radius() / r_poly << endl ;
754 xk.
set(1).set_dzpuis(-1) ;
757 xk.
set(2).set_dzpuis(-1) ;
760 xk.
set(3).set_dzpuis(-1) ;
761 xk.std_spectral_base() ;
768 double lambda = double(1) / double(3) ;
775 + (lambda/2./(lambda+1)) * (d_khi + x_d_w) ;
781 if ( (
beta(1).get_etat() == ETATZERO) && (
beta(2).get_etat() == ETATZERO) ) {
virtual double mass_g() const
Gravitational mass.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Scalar dzeta
Metric potential .
virtual double mean_radius() const
Mean star radius from the area .
Scalar ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta ...
Scalar a_car
Square of the metric factor A.
Map & mp
Mapping associated with the star.
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
double * p_r_isco
Circumferential radius of the ISCO.
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
virtual double r_circ_merid() const
Crcumferential meridional radius.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
double * p_f_eq
Orbital frequency at the equator.
Cmp sqrt(const Cmp &)
Square root.
double * p_r_circ
Circumferential radius (equator)
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Scalar tnphi
Component of the shift vector.
double * p_mom_quad
Quadrupole moment.
Scalar bbb
Metric factor B.
const Eos & eos
Equation of state of the stellar matter.
Standard units of space, time and mass.
Equation of state base class.
virtual void comp_p_from_cartesian(const Scalar &v_x, const Scalar &v_y, Scalar &v_p) const =0
Computes the Spherical component (with respect to bvect_spher ) of a vector given by its cartesian c...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual double z_pole() const
Redshift factor at North pole.
double & set(int i)
Read/write of a particular element (index i) (1D case)
double * p_f_isco
Orbital frequency of the ISCO.
Tensor field of valence 0 (or component of a tensorial field).
virtual double grv2() const
Error on the virial identity GRV2.
Base class for coordinate mappings.
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
virtual void del_deriv() const
Deletes all the derived quantities.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
virtual void sauve(FILE *) const
Save in a binary file.
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
double omega
Rotation angular velocity ([f_unit] )
Scalar nbar
Baryon density in the fluid frame.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Tensor field of valence 1.
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
virtual double area() const
Integrated surface area in .
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Class for isolated rotating stars.
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg ...
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
void operator=(const Star_rot &)
Assignment to another Star_rot.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
void operator=(const Star &)
Assignment to another Star.
int nzet
Number of domains of *mp occupied by the star.
virtual ~Star_rot()
Destructor.
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Scalar nphi
Metric coefficient .
double * p_z_eqf
Forward redshift factor at equator.
virtual void sauve(FILE *) const
Save in a file.
Scalar ener
Total energy density in the fluid frame.
double * p_grv2
Error on the virial identity GRV2.
double * p_area
Integrated surface area.
double * p_r_circ_merid
Circumferential radius (meridian)
Scalar press
Fluid pressure.
virtual double mass_b() const
Baryon mass.
double * p_espec_isco
Specific energy of a particle on the ISCO.
double get_kap() const
Returns the pressure coefficient (cf.
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
virtual double mom_quad() const
Quadrupole moment.
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
virtual double angu_mom() const
Angular momentum.
double * p_z_pole
Redshift factor at North pole.
Polytropic equation of state (relativistic case).
virtual void display_poly(ostream &) const
Display in polytropic units.
Star_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
int get_nzone() const
Returns the number of domains.
double * p_grv3
Error on the virial identity GRV3.
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
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.
virtual double r_circ() const
Circumferential equatorial radius.
double * p_angu_mom
Angular momentum.
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Scalar b_car
Square of the metric factor B.
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.
double * p_aplat
Flatening r_pole/r_eq.
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
virtual void del_deriv() const
Deletes all the derived quantities.
virtual double z_eqb() const
Backward redshift factor at equator.
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Scalar nn
Lapse function N .
Coord y
y coordinate centered on the grid
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
int get_taille() const
Gives the total size (ie dim.taille)
Scalar uuu
Norm of u_euler.
double ray_eq() const
Coordinate radius at , [r_unit].
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Coord x
x coordinate centered on the grid
void div_rsint()
Division by everywhere; dzpuis is not changed.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
virtual void sauve(FILE *) const
Save in a file.
virtual void sauve(FILE *) const
Save in a file.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
Scalar tggg
Metric potential .
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift ...
virtual double aplat() const
Flatening r_pole/r_eq.
Tbl * p_surf_grav
Surface gravity (along the theta direction)
Coord z
z coordinate centered on the grid
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
virtual double z_eqf() const
Forward redshift factor at equator.
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
virtual double tsw() const
Ratio T/W.
double * p_z_eqb
Backward redshift factor at equator.