gravastar.C
/*
* Methods of class Gravastar
*
* (see file gravastar.h for documentation).
*
*/
/*
* Copyright (c) 2010 Frederic Vincent
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2
* as published by the Free Software Foundation.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char gravastar_C[] = "$Header: /cvsroot/Lorene/C++/Source/Gravastar/gravastar.C,v 1.4 2016/09/19 15:26:23 j_novak Exp $" ;
// C++ headers
//#include <>
// C headers
#include <cmath>
// Lorene headers
#include "gravastar.h"
#include "eos.h"
#include "utilitaires.h"
#include "param.h"
#include "unites.h"
#include "nbr_spx.h"
//---------------------//
// Constructor //
//---------------------//
namespace Lorene {
Gravastar::Gravastar(Map& mpi, int nzet_i, const Eos& eos_i, const double rho_core_i)
: Star_rot(mpi,nzet_i,true,eos_i), rho_core(rho_core_i)
{}
//------------//
// Destructor //
//------------//
Gravastar::~Gravastar(){
del_deriv() ;
}
//----------------------//
// Eq. of state //
//----------------------//
void Gravastar::equation_of_state() {
Scalar ent_eos = ent ;
// Slight rescale of the enthalpy field in case of 2 domains inside the
// star
double epsilon = 1.e-12 ;
const Mg3d* mg = mp.get_mg() ;
int nz = mg->get_nzone() ;
Mtbl xi(mg) ;
xi.set_etat_qcq() ;
for (int l=0; l<nz; l++) {
xi.t[l]->set_etat_qcq() ;
for (int k=0; k<mg->get_np(l); k++) {
for (int j=0; j<mg->get_nt(l); j++) {
for (int i=0; i<mg->get_nr(l); i++) {
xi.set(l,k,j,i) =
mg->get_grille3d(l)->x[i] ;
}
}
}
}
Scalar fact_ent(mp) ;
fact_ent.allocate_all() ;
fact_ent.set_domain(0) = 1 + epsilon * xi(0) * xi(0) ;
fact_ent.set_domain(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
for (int l=nzet; l<nz; l++) {
fact_ent.set_domain(l) = 1 ;
}
if (nzet > 1) {
if (nzet > 2) {
cout << "Gravastar::equation_of_state: not ready yet for nzet > 2 !"
<< endl ;
}
ent_eos = fact_ent * ent_eos ;
ent_eos.std_spectral_base() ;
}
/*
Call to the Eos
There is no Eos inside the core where p=-rho=cst. Thus the most internal domain in treated separately.
Inside the crust, the Eos is called domain by domain.
*/
Scalar tempo(mp) ;
//nbar.set_etat_qcq() ;
nbar = 0 ;
for (int l=1; l<nzet; l++) {
Param par ; // Paramater for multi-domain equation of state
par.add_int(l) ;
tempo = eos.nbar_ent(ent_eos, 1, l, &par) ;
//cout << "tempo eos= " << tempo << endl;
nbar = nbar + tempo ;
}
//ener.set_etat_qcq() ;
ener = rho_core ;
ener.annule(1,nz-1);
for (int l=1; l<nzet; l++) {
Param par ; // Paramater for multi-domain equation of state
par.add_int(l) ;
tempo = eos.ener_ent(ent_eos, 1, l, &par) ;
ener = ener + tempo ;
}
//press.set_etat_qcq() ;
press = -rho_core ;
press.annule(1,nz-1);
for (int l=1; l<nzet; l++) {
Param par ; // Paramater for multi-domain equation of state
par.add_int(l) ;
tempo = eos.press_ent(ent_eos, 1, l, &par) ;
press = press + tempo ;
}
// Set the bases for spectral expansion
nbar.std_spectral_base() ;
ener.std_spectral_base() ;
press.std_spectral_base() ;
// The derived quantities are obsolete
del_deriv() ;
}
// Printing
// --------
ostream& Gravastar::operator>>(ostream& ost) const {
using namespace Unites ;
ost << endl ;
ost << "Number of domains occupied by the star : " << nzet << endl ;
ost << "Equation of state : " << endl ;
ost << eos << endl ;
const Mg3d* mg = mp.get_mg() ;
int l_cr = 0, i_cr = mg->get_nr(l_cr) - 1, j_cr = mg->get_nt(l_cr) - 1, k_cr = 0 ;
ost << endl << "Inner crust enthalpy : " << ent.val_grid_point(l_cr, k_cr, j_cr, i_cr) << " c^2" << endl ;
ost << "Central proper energy density : " << ener.val_grid_point(0,0,0,0)
<< " rho_nuc c^2" << endl ;
ost << "Central pressure : " << press.val_grid_point(0,0,0,0)
<< " rho_nuc c^2" << endl ;
ost << endl ;
ost << "Central lapse N : " << nn.val_grid_point(0,0,0,0) << endl ;
// ost << "Central value of lnq : " << lnq.val_grid_point(0,0,0,0) << endl ;
ost << endl
<< "Coordinate equatorial radius (phi=0) a1 = "
<< ray_eq()/km << " km" << endl ;
ost << "Coordinate equatorial radius (phi=pi/2) a2 = "
<< ray_eq_pis2()/km << " km" << endl ;
ost << "Coordinate equatorial radius (phi=pi): "
<< ray_eq_pi()/km << " km" << endl ;
ost << "Coordinate polar radius a3 = "
<< ray_pole()/km << " km" << endl ;
ost << "Axis ratio a2/a1 = " << ray_eq_pis2() / ray_eq()
<< " a3/a1 = " << ray_pole() / ray_eq() << endl ;
double omega_c = get_omega_c() ;
ost << endl ;
if (omega != __infinity) {
ost << "Uniformly rotating star" << endl ;
ost << "-----------------------" << endl ;
double freq = omega / (2.*M_PI) ;
ost << "Omega : " << omega * f_unit
<< " rad/s f : " << freq * f_unit << " Hz" << endl ;
ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
<< endl ;
}
else {
ost << "Differentially rotating star" << endl ;
ost << "----------------------------" << endl ;
double freq = omega_c / (2.*M_PI) ;
ost << "Central value of Omega : " << omega_c * f_unit
<< " rad/s f : " << freq * f_unit << " Hz" << endl ;
ost << "Central rotation period : " << 1000. / (freq * f_unit) << " ms"
<< endl ;
}
// double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
//ost << "Compactness G M_g /(c^2 R_circ) : " << compact << endl ;
double nphi_c = nphi.val_grid_point(0, 0, 0, 0) ;
if ( (omega_c==0) && (nphi_c==0) ) {
ost << "Central N^phi : " << nphi_c << endl ;
}
else{
ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
}
ost << "Error on the virial identity GRV2 : " << grv2() << endl ;
double xgrv3 = grv3(&ost) ;
ost << "Error on the virial identity GRV3 : " << xgrv3 << endl ;
// cout << "ICI" << endl;
double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
/ double(1.e38) ) ;
//cout << "LA" << endl;
ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
<< endl ;
/*ost << "Q / (M R_circ^2) : "
<< mom_quad() / ( mass_g() * pow( r_circ(), 2. ) ) << endl ;
ost << "c^4 Q / (G^2 M^3) : "
<< mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
<< endl ; */
ost << "Angular momentum J : "
<< angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
<< endl ;
/* ost << "c J / (G M^2) : "
<< angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ; */
if (omega != __infinity) {
double mom_iner = angu_mom() / omega ;
double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
/ double(1.e38) ) ;
ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
<< endl ;
}
ost << "Ratio T/W : " << tsw() << endl ;
ost << "Circumferential equatorial radius R_circ : "
<< r_circ()/km << " km" << endl ;
ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
<< endl ;
ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
//cout << "ICI" << endl;
int lsurf = nzet - 1;
int nt = mp.get_mg()->get_nt(lsurf) ;
int nr = mp.get_mg()->get_nr(lsurf) ;
ost << "Equatorial value of the velocity U: "
<< uuu.val_grid_point(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
ost << "Redshift at the equator (forward) : " << z_eqf() << endl ;
ost << "Redshift at the equator (backward): " << z_eqb() << endl ;
ost << "Redshift at the pole : " << z_pole() << endl ;
ost << "Central value of log(N) : "
<< logn.val_grid_point(0, 0, 0, 0) << endl ;
ost << "Central value of dzeta=log(AN) : "
<< dzeta.val_grid_point(0, 0, 0, 0) << endl ;
if ( (omega_c==0) && (nphi_c==0) ) {
ost << "Central N^phi : " << nphi_c << endl ;
}
else{
ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
}
ost << " ... w_shift (NB: components in the star Cartesian frame) [c] : "
<< endl
<< w_shift(1).val_grid_point(0, 0, 0, 0) << " "
<< w_shift(2).val_grid_point(0, 0, 0, 0) << " "
<< w_shift(3).val_grid_point(0, 0, 0, 0) << endl ;
ost << "Central value of khi_shift [km c] : "
<< khi_shift.val_grid_point(0, 0, 0, 0) / km << endl ;
ost << "Central value of B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
Tbl diff_a_b = diffrel( a_car, b_car ) ;
ost <<
"Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
<< endl ;
for (int l=0; l<diff_a_b.get_taille(); l++) {
ost << diff_a_b(l) << " " ;
}
ost << endl ;
// Approximate formula for R_isco = 6 R_g (1-(2/3)^1.5 j )
// up to the first order in j
/*double jdimless = angu_mom() / ( ggrav * pow(mass_g(), 2.) ) ;
double r_grav = ggrav * mass_g() ;
double r_isco_appr = 6. * r_grav * ( 1. - pow(2./3.,1.5) * jdimless ) ;*/
// Approximate formula for the ISCO frequency
// freq_ms = 6^{-1.5}/2pi/R_g (1+11*6^(-1.5) j )
// up to the first order in j
/* double f_isco_appr = ( 1. + 11. /6. /sqrt(6.) * jdimless ) / r_grav /
(12. * M_PI ) / sqrt(6.) ;*/
/* ost << endl << "Innermost stable circular orbit (ISCO) : " << endl ;
double xr_isco = r_isco(&ost) ;
ost <<" circumferential radius r_isco = "
<< xr_isco / km << " km" << endl ;
ost <<" (approx. 6M + 1st order in j : "
<< r_isco_appr / km << " km)" << endl ;
ost <<" (approx. 6M : "
<< 6. * r_grav / km << " km)" << endl ;
ost <<" orbital frequency f_isco = "
<< f_isco() * f_unit << " Hz" << endl ;
ost <<" (approx. 1st order in j : "
<< f_isco_appr * f_unit << " Hz)" << endl ;*/
return ost ;
}
}
Lorene::Eos
Equation of state base class.
Definition eos.h:190
Lorene::Gravastar::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition gravastar.C:185
Lorene::Gravastar::~Gravastar
~Gravastar()
Destructor.
Definition gravastar.C:59
Lorene::Gravastar::rho_core
double rho_core
Energy density in gravastar's core.
Definition gravastar.h:57
Lorene::Gravastar::Gravastar
Gravastar(Map &mp_i, int nzet_i, const Eos &eos_i, const double rho_core_i)
Standard constructor.
Definition gravastar.C:51
Lorene::Gravastar::equation_of_state
void equation_of_state()
Allows to computes the proper baryon and energy density, as well as pressure from the enthalpy,...
Definition gravastar.C:67
Lorene::Grille3d::x
double * x
Array of values of at the nr collocation points.
Definition grilles.h:209
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_grille3d
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Mtbl
Multi-domain array.
Definition mtbl.h:118
Lorene::Mtbl::set
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition mtbl.h:221
Lorene::Mtbl::t
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Lorene::Mtbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::std_spectral_base
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
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
Lorene::Scalar::allocate_all
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
Lorene::Star_rot::Star_rot
Star_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition star_rot.C:86
Lorene::Star_rot::mom_quad
virtual double mom_quad() const
Quadrupole moment.
Definition star_rot_global.C:521
Lorene::Star_rot::z_pole
virtual double z_pole() const
Redshift factor at North pole.
Definition star_rot_global.C:502
Lorene::Star_rot::z_eqb
virtual double z_eqb() const
Backward redshift factor at equator.
Definition star_rot_global.C:471
Lorene::Star_rot::w_shift
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:150
Lorene::Star_rot::b_car
Scalar b_car
Square of the metric factor B.
Definition star_rot.h:110
Lorene::Star_rot::get_omega_c
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] ).
Definition star_rot.C:640
Lorene::Star_rot::omega
double omega
Rotation angular velocity ([f_unit] ).
Definition star_rot.h:101
Lorene::Star_rot::r_circ
virtual double r_circ() const
Circumferential radius.
Definition star_rot_global.C:342
Lorene::Star_rot::uuu
Scalar uuu
Norm of u_euler.
Definition star_rot.h:121
Lorene::Star_rot::z_eqf
virtual double z_eqf() const
Forward redshift factor at equator.
Definition star_rot_global.C:443
Lorene::Star_rot::nphi
Scalar nphi
Metric coefficient .
Definition star_rot.h:113
Lorene::Star_rot::aplat
virtual double aplat() const
Flatening r_pole/r_eq.
Definition star_rot_global.C:426
Lorene::Star_rot::tsw
virtual double tsw() const
Ratio T/W.
Definition star_rot_global.C:176
Lorene::Star_rot::grv3
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Definition star_rot_global.C:235
Lorene::Star_rot::dzeta
Scalar dzeta
Metric potential .
Definition star_rot.h:134
Lorene::Star_rot::a_car
Scalar a_car
Square of the metric factor A.
Definition star_rot.h:104
Lorene::Star_rot::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star_rot.C:297
Lorene::Star_rot::grv2
virtual double grv2() const
Error on the virial identity GRV2.
Definition star_rot_global.C:204
Lorene::Star_rot::angu_mom
virtual double angu_mom() const
Angular momentum.
Definition star_rot_global.C:150
Lorene::Star_rot::khi_shift
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:160
Lorene::Star::ener
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
Lorene::Star::logn
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Lorene::Star::nn
Scalar nn
Lapse function N .
Definition star.h:225
Lorene::Star::eos
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Lorene::Star::nbar
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
Lorene::Star::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition star_global.C:108
Lorene::Star::ray_eq_pis2
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Definition star_global.C:138
Lorene::Star::press
Scalar press
Fluid pressure.
Definition star.h:194
Lorene::Star::ray_eq_pi
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition star_global.C:186
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition star.h:190
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition star.h:180
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Lorene::Star::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition star_global.C:278
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Unites
Standard units of space, time and mass.