28char gravastar_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Gravastar/gravastar.C,v 1.4 2016/09/19 15:26:23 j_novak Exp $" ;
40#include "utilitaires.h"
52 :
Star_rot(mpi,nzet_i,true,eos_i), rho_core(rho_core_i)
76 double epsilon = 1.e-12 ;
83 for (
int l=0; l<nz; l++) {
85 for (
int k=0; k<mg->
get_np(l); k++) {
86 for (
int j=0; j<mg->
get_nt(l); j++) {
87 for (
int i=0; i<mg->
get_nr(l); i++) {
99 fact_ent.
set_domain(0) = 1 + epsilon * xi(0) * xi(0) ;
100 fact_ent.
set_domain(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
102 for (
int l=
nzet; l<nz; l++) {
110 cout <<
"Gravastar::equation_of_state: not ready yet for nzet > 2 !"
114 ent_eos = fact_ent * ent_eos ;
130 for (
int l=1; l<
nzet; l++) {
146 for (
int l=1; l<
nzet; l++) {
160 for (
int l=1; l<
nzet; l++) {
191 ost <<
"Number of domains occupied by the star : " <<
nzet << endl ;
193 ost <<
"Equation of state : " << endl ;
197 int l_cr = 0, i_cr = mg->
get_nr(l_cr) - 1, j_cr = mg->
get_nt(l_cr) - 1, k_cr = 0 ;
199 ost << endl <<
"Inner crust enthalpy : " <<
ent.
val_grid_point(l_cr, k_cr, j_cr, i_cr) <<
" c^2" << endl ;
202 <<
" rho_nuc c^2" << endl ;
204 <<
" rho_nuc c^2" << endl ;
211 <<
"Coordinate equatorial radius (phi=0) a1 = "
212 <<
ray_eq()/km <<
" km" << endl ;
213 ost <<
"Coordinate equatorial radius (phi=pi/2) a2 = "
215 ost <<
"Coordinate equatorial radius (phi=pi): "
217 ost <<
"Coordinate polar radius a3 = "
228 if (
omega != __infinity) {
229 ost <<
"Uniformly rotating star" << endl ;
230 ost <<
"-----------------------" << endl ;
232 double freq =
omega / (2.*M_PI) ;
233 ost <<
"Omega : " <<
omega * f_unit
234 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
235 ost <<
"Rotation period : " << 1000. / (freq * f_unit) <<
" ms"
240 ost <<
"Differentially rotating star" << endl ;
241 ost <<
"----------------------------" << endl ;
243 double freq = omega_c / (2.*M_PI) ;
244 ost <<
"Central value of Omega : " << omega_c * f_unit
245 <<
" rad/s f : " << freq * f_unit <<
" Hz" << endl ;
246 ost <<
"Central rotation period : " << 1000. / (freq * f_unit) <<
" ms"
254 if ( (omega_c==0) && (nphi_c==0) ) {
255 ost <<
"Central N^phi : " << nphi_c << endl ;
258 ost <<
"Central N^phi/Omega : " << nphi_c / omega_c << endl ;
261 ost <<
"Error on the virial identity GRV2 : " <<
grv2() << endl ;
262 double xgrv3 =
grv3(&ost) ;
263 ost <<
"Error on the virial identity GRV3 : " << xgrv3 << endl ;
266 double mom_quad_38si =
mom_quad() * rho_unit * (
pow(r_unit,
double(5.))
269 ost <<
"Quadrupole moment Q : " << mom_quad_38si <<
" 10^38 kg m^2"
278 ost <<
"Angular momentum J : "
279 <<
angu_mom()/( qpig / (4* M_PI) * msol*msol) <<
" G M_sol^2 / c"
284if (
omega != __infinity) {
286 double mom_iner_38si = mom_iner * rho_unit * (
pow(r_unit,
double(5.))
288 ost <<
"Moment of inertia: " << mom_iner_38si <<
" 10^38 kg m^2"
292 ost <<
"Ratio T/W : " <<
tsw() << endl ;
293 ost <<
"Circumferential equatorial radius R_circ : "
294 <<
r_circ()/km <<
" km" << endl ;
295 ost <<
"Coordinate equatorial radius r_eq : " <<
ray_eq()/km <<
" km"
297 ost <<
"Flattening r_pole/r_eq : " <<
aplat() << endl ;
301 int lsurf =
nzet - 1;
304 ost <<
"Equatorial value of the velocity U: "
306 ost <<
"Redshift at the equator (forward) : " <<
z_eqf() << endl ;
307 ost <<
"Redshift at the equator (backward): " <<
z_eqb() << endl ;
308 ost <<
"Redshift at the pole : " <<
z_pole() << endl ;
311 ost <<
"Central value of log(N) : "
314 ost <<
"Central value of dzeta=log(AN) : "
317 if ( (omega_c==0) && (nphi_c==0) ) {
318 ost <<
"Central N^phi : " << nphi_c << endl ;
321 ost <<
"Central N^phi/Omega : " << nphi_c / omega_c << endl ;
324 ost <<
" ... w_shift (NB: components in the star Cartesian frame) [c] : "
326 <<
w_shift(1).val_grid_point(0, 0, 0, 0) <<
" "
327 <<
w_shift(2).val_grid_point(0, 0, 0, 0) <<
" "
328 <<
w_shift(3).val_grid_point(0, 0, 0, 0) << endl ;
330 ost <<
"Central value of khi_shift [km c] : "
337 "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
340 ost << diff_a_b(l) <<
" " ;
Equation of state base class.
Cmp nbar_ent(const Cmp &ent, int nzet, int l_min=0, const Param *par=0x0) const
Computes the baryon density field from the log-enthalpy field and extra parameters.
Cmp ener_ent(const Cmp &ent, int nzet, int l_min=0, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy and extra parameters.
Cmp press_ent(const Cmp &ent, int nzet, int l_min=0, const Param *par=0x0) const
Computes the pressure from the log-enthalpy and extra parameters.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
double rho_core
Energy density in gravastar's core.
Gravastar(Map &mp_i, int nzet_i, const Eos &eos_i, const double rho_core_i)
Standard constructor.
void equation_of_state()
Allows to computes the proper baryon and energy density, as well as pressure from the enthalpy,...
double * x
Array of values of at the nr collocation points.
Base class for coordinate mappings.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Tensor field of valence 0 (or component of a tensorial field).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Tbl & set_domain(int l)
Read/write of the value in a given domain.
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Class for isolated rotating stars.
virtual double mom_quad() const
Quadrupole moment.
virtual double z_pole() const
Redshift factor at North pole.
virtual double z_eqb() const
Backward redshift factor at equator.
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Scalar b_car
Square of the metric factor B.
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
double omega
Rotation angular velocity ([f_unit] )
virtual double r_circ() const
Circumferential radius.
Scalar uuu
Norm of u_euler.
virtual double z_eqf() const
Forward redshift factor at equator.
Scalar nphi
Metric coefficient .
virtual double aplat() const
Flatening r_pole/r_eq.
virtual double tsw() const
Ratio T/W.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Scalar dzeta
Metric potential .
Scalar a_car
Square of the metric factor A.
virtual void del_deriv() const
Deletes all the derived quantities.
virtual double grv2() const
Error on the virial identity GRV2.
virtual double angu_mom() const
Angular momentum.
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Scalar ener
Total energy density in the fluid frame.
Scalar logn
Logarithm of the lapse N .
Scalar nn
Lapse function N .
const Eos & eos
Equation of state of the stellar matter.
Scalar nbar
Baryon density in the fluid frame.
double ray_eq() const
Coordinate radius at , [r_unit].
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Scalar press
Fluid pressure.
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Map & mp
Mapping associated with the star.
int nzet
Number of domains of *mp occupied by the star.
double ray_pole() const
Coordinate radius at [r_unit].
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int get_taille() const
Gives the total size (ie dim.taille)
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Cmp pow(const Cmp &, int)
Power .
Standard units of space, time and mass.