LORENE
gravastar.C
1/*
2 * Methods of class Gravastar
3 *
4 * (see file gravastar.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2010 Frederic Vincent
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char gravastar_C[] = "$Header: /cvsroot/Lorene/C++/Source/Gravastar/gravastar.C,v 1.4 2016/09/19 15:26:23 j_novak Exp $" ;
29
30
31// C++ headers
32//#include <>
33
34// C headers
35#include <cmath>
36
37// Lorene headers
38#include "gravastar.h"
39#include "eos.h"
40#include "utilitaires.h"
41#include "param.h"
42#include "unites.h"
43#include "nbr_spx.h"
44
45
46 //---------------------//
47 // Constructor //
48 //---------------------//
49
50namespace Lorene {
51Gravastar::Gravastar(Map& mpi, int nzet_i, const Eos& eos_i, const double rho_core_i)
52 : Star_rot(mpi,nzet_i,true,eos_i), rho_core(rho_core_i)
53{}
54
55 //------------//
56 // Destructor //
57 //------------//
58
62
63 //----------------------//
64 // Eq. of state //
65 //----------------------//
66
68
69 Scalar ent_eos = ent ;
70
71
72 // Slight rescale of the enthalpy field in case of 2 domains inside the
73 // star
74
75
76 double epsilon = 1.e-12 ;
77
78 const Mg3d* mg = mp.get_mg() ;
79 int nz = mg->get_nzone() ;
80
81 Mtbl xi(mg) ;
82 xi.set_etat_qcq() ;
83 for (int l=0; l<nz; l++) {
84 xi.t[l]->set_etat_qcq() ;
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++) {
88 xi.set(l,k,j,i) =
89 mg->get_grille3d(l)->x[i] ;
90 }
91 }
92 }
93
94 }
95
96 Scalar fact_ent(mp) ;
97 fact_ent.allocate_all() ;
98
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) ;
101
102 for (int l=nzet; l<nz; l++) {
103 fact_ent.set_domain(l) = 1 ;
104 }
105
106 if (nzet > 1) {
107
108 if (nzet > 2) {
109
110 cout << "Gravastar::equation_of_state: not ready yet for nzet > 2 !"
111 << endl ;
112 }
113
114 ent_eos = fact_ent * ent_eos ;
115 ent_eos.std_spectral_base() ;
116 }
117
118
119
120 /*
121 Call to the Eos
122 There is no Eos inside the core where p=-rho=cst. Thus the most internal domain in treated separately.
123 Inside the crust, the Eos is called domain by domain.
124 */
125
126 Scalar tempo(mp) ;
127
128 //nbar.set_etat_qcq() ;
129 nbar = 0 ;
130 for (int l=1; l<nzet; l++) {
131
132 Param par ; // Paramater for multi-domain equation of state
133 par.add_int(l) ;
134
135
136
137 tempo = eos.nbar_ent(ent_eos, 1, l, &par) ;
138 //cout << "tempo eos= " << tempo << endl;
139 nbar = nbar + tempo ;
140
141 }
142
143 //ener.set_etat_qcq() ;
144 ener = rho_core ;
145 ener.annule(1,nz-1);
146 for (int l=1; l<nzet; l++) {
147
148 Param par ; // Paramater for multi-domain equation of state
149 par.add_int(l) ;
150
151 tempo = eos.ener_ent(ent_eos, 1, l, &par) ;
152
153 ener = ener + tempo ;
154
155 }
156
157 //press.set_etat_qcq() ;
158 press = -rho_core ;
159 press.annule(1,nz-1);
160 for (int l=1; l<nzet; l++) {
161
162 Param par ; // Paramater for multi-domain equation of state
163 par.add_int(l) ;
164
165 tempo = eos.press_ent(ent_eos, 1, l, &par) ;
166
167 press = press + tempo ;
168
169 }
170
171
172 // Set the bases for spectral expansion
176
177 // The derived quantities are obsolete
178 del_deriv() ;
179
180}
181
182// Printing
183// --------
184
185ostream& Gravastar::operator>>(ostream& ost) const {
186
187 using namespace Unites ;
188
189 ost << endl ;
190
191 ost << "Number of domains occupied by the star : " << nzet << endl ;
192
193 ost << "Equation of state : " << endl ;
194 ost << eos << endl ;
195
196 const Mg3d* mg = mp.get_mg() ;
197 int l_cr = 0, i_cr = mg->get_nr(l_cr) - 1, j_cr = mg->get_nt(l_cr) - 1, k_cr = 0 ;
198
199 ost << endl << "Inner crust enthalpy : " << ent.val_grid_point(l_cr, k_cr, j_cr, i_cr) << " c^2" << endl ;
200
201 ost << "Central proper energy density : " << ener.val_grid_point(0,0,0,0)
202 << " rho_nuc c^2" << endl ;
203 ost << "Central pressure : " << press.val_grid_point(0,0,0,0)
204 << " rho_nuc c^2" << endl ;
205
206 ost << endl ;
207 ost << "Central lapse N : " << nn.val_grid_point(0,0,0,0) << endl ;
208 // ost << "Central value of lnq : " << lnq.val_grid_point(0,0,0,0) << endl ;
209
210 ost << endl
211 << "Coordinate equatorial radius (phi=0) a1 = "
212 << ray_eq()/km << " km" << endl ;
213 ost << "Coordinate equatorial radius (phi=pi/2) a2 = "
214 << ray_eq_pis2()/km << " km" << endl ;
215 ost << "Coordinate equatorial radius (phi=pi): "
216 << ray_eq_pi()/km << " km" << endl ;
217 ost << "Coordinate polar radius a3 = "
218 << ray_pole()/km << " km" << endl ;
219 ost << "Axis ratio a2/a1 = " << ray_eq_pis2() / ray_eq()
220 << " a3/a1 = " << ray_pole() / ray_eq() << endl ;
221
223
224 double omega_c = get_omega_c() ;
225
226 ost << endl ;
227
228 if (omega != __infinity) {
229 ost << "Uniformly rotating star" << endl ;
230 ost << "-----------------------" << endl ;
231
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"
236 << endl ;
237
238 }
239 else {
240 ost << "Differentially rotating star" << endl ;
241 ost << "----------------------------" << endl ;
242
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"
247 << endl ;
248 }
249
250 // double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
251 //ost << "Compactness G M_g /(c^2 R_circ) : " << compact << endl ;
252
253 double nphi_c = nphi.val_grid_point(0, 0, 0, 0) ;
254 if ( (omega_c==0) && (nphi_c==0) ) {
255 ost << "Central N^phi : " << nphi_c << endl ;
256 }
257 else{
258 ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
259 }
260
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 ;
264
265 // cout << "ICI" << endl;
266 double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
267 / double(1.e38) ) ;
268//cout << "LA" << endl;
269 ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
270 << endl ;
271
272 /*ost << "Q / (M R_circ^2) : "
273 << mom_quad() / ( mass_g() * pow( r_circ(), 2. ) ) << endl ;
274 ost << "c^4 Q / (G^2 M^3) : "
275 << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
276 << endl ; */
277
278 ost << "Angular momentum J : "
279 << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
280 << endl ;
281 /* ost << "c J / (G M^2) : "
282 << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ; */
283
284if (omega != __infinity) {
285 double mom_iner = angu_mom() / omega ;
286 double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
287 / double(1.e38) ) ;
288 ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
289 << endl ;
290 }
291
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"
296 << endl ;
297 ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
298
299 //cout << "ICI" << endl;
300
301 int lsurf = nzet - 1;
302 int nt = mp.get_mg()->get_nt(lsurf) ;
303 int nr = mp.get_mg()->get_nr(lsurf) ;
304 ost << "Equatorial value of the velocity U: "
305 << uuu.val_grid_point(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
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 ;
309
310
311 ost << "Central value of log(N) : "
312 << logn.val_grid_point(0, 0, 0, 0) << endl ;
313
314 ost << "Central value of dzeta=log(AN) : "
315 << dzeta.val_grid_point(0, 0, 0, 0) << endl ;
316
317 if ( (omega_c==0) && (nphi_c==0) ) {
318 ost << "Central N^phi : " << nphi_c << endl ;
319 }
320 else{
321 ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
322 }
323
324 ost << " ... w_shift (NB: components in the star Cartesian frame) [c] : "
325 << endl
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 ;
329
330 ost << "Central value of khi_shift [km c] : "
331 << khi_shift.val_grid_point(0, 0, 0, 0) / km << endl ;
332
333 ost << "Central value of B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
334
335 Tbl diff_a_b = diffrel( a_car, b_car ) ;
336 ost <<
337 "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
338 << endl ;
339 for (int l=0; l<diff_a_b.get_taille(); l++) {
340 ost << diff_a_b(l) << " " ;
341 }
342 ost << endl ;
343
344 // Approximate formula for R_isco = 6 R_g (1-(2/3)^1.5 j )
345 // up to the first order in j
346
347 /*double jdimless = angu_mom() / ( ggrav * pow(mass_g(), 2.) ) ;
348 double r_grav = ggrav * mass_g() ;
349 double r_isco_appr = 6. * r_grav * ( 1. - pow(2./3.,1.5) * jdimless ) ;*/
350
351 // Approximate formula for the ISCO frequency
352 // freq_ms = 6^{-1.5}/2pi/R_g (1+11*6^(-1.5) j )
353 // up to the first order in j
354
355 /* double f_isco_appr = ( 1. + 11. /6. /sqrt(6.) * jdimless ) / r_grav /
356 (12. * M_PI ) / sqrt(6.) ;*/
357
358 /* ost << endl << "Innermost stable circular orbit (ISCO) : " << endl ;
359 double xr_isco = r_isco(&ost) ;
360 ost <<" circumferential radius r_isco = "
361 << xr_isco / km << " km" << endl ;
362 ost <<" (approx. 6M + 1st order in j : "
363 << r_isco_appr / km << " km)" << endl ;
364 ost <<" (approx. 6M : "
365 << 6. * r_grav / km << " km)" << endl ;
366 ost <<" orbital frequency f_isco = "
367 << f_isco() * f_unit << " Hz" << endl ;
368 ost <<" (approx. 1st order in j : "
369 << f_isco_appr * f_unit << " Hz)" << endl ;*/
370
371
372 return ost ;
373
374}
375}
Equation of state base class.
Definition eos.h:190
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.
Definition eos.C:338
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.
Definition eos.C:363
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.
Definition eos.C:385
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition gravastar.C:185
~Gravastar()
Destructor.
Definition gravastar.C:59
double rho_core
Energy density in gravastar's core.
Definition gravastar.h:57
Gravastar(Map &mp_i, int nzet_i, const Eos &eos_i, const double rho_core_i)
Standard constructor.
Definition gravastar.C:51
void equation_of_state()
Allows to computes the proper baryon and energy density, as well as pressure from the enthalpy,...
Definition gravastar.C:67
double * x
Array of values of at the nr collocation points.
Definition grilles.h:209
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Multi-domain array.
Definition mtbl.h:118
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition mtbl.h:221
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Parameter storage.
Definition param.h:125
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
Class for isolated rotating stars.
Definition star_rot.h:86
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.
Definition star_rot.h:150
Scalar b_car
Square of the metric factor B.
Definition star_rot.h:110
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
Definition star_rot.C:640
double omega
Rotation angular velocity ([f_unit] )
Definition star_rot.h:101
virtual double r_circ() const
Circumferential radius.
Scalar uuu
Norm of u_euler.
Definition star_rot.h:121
virtual double z_eqf() const
Forward redshift factor at equator.
Scalar nphi
Metric coefficient .
Definition star_rot.h:113
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 .
Definition star_rot.h:134
Scalar a_car
Square of the metric factor A.
Definition star_rot.h:104
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star_rot.C:297
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.
Definition star_rot.h:160
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Scalar nn
Lapse function N .
Definition star.h:225
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
double ray_eq() const
Coordinate radius at , [r_unit].
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Scalar press
Fluid pressure.
Definition star.h:194
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Scalar ent
Log-enthalpy.
Definition star.h:190
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
double ray_pole() const
Coordinate radius at [r_unit].
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.