etoile_rot.C
/*
* Methods for the class Etoile_rot
*
* (see file etoile.h for documentation)
*/
/*
* Copyright (c) 2000-2001 Eric Gourgoulhon
*
* 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 as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* 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 etoile_rot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $" ;
/*
* $Id: etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $
* $Log: etoile_rot.C,v $
* Revision 1.7 2015/12/03 14:17:24 j_novak
* Check added for the computation of area (thanks S. Koeppel).
*
* Revision 1.6 2015/06/10 14:37:44 a_sourie
* Corrected the formula for the quadrupole.
*
* Revision 1.5 2014/10/13 08:52:59 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2004/03/25 10:29:07 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.3 2001/12/06 15:11:43 jl_zdunik
* Introduction of the new function f_eq() in the class Etoile_rot
*
* Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
*
* All writing/reading to a binary file are now performed according to
* the big endian convention, whatever the system is big endian or
* small endian, thanks to the functions fwrite_be and fread_be
*
* Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
* LORENE
*
* Revision 2.17 2001/10/24 15:36:20 eric
* Ajout de la fonction display_poly.
*
* Revision 2.16 2001/10/16 14:49:02 eric
* Appel de get_omega_c() pour avoir la valeur centrale de Omega.
* Affichage different si rotation differentielle.
*
* Revision 2.15 2001/09/13 08:32:01 eric
* Ajout du facteur de compacite M/R dans l'affichage.
*
* Revision 2.14 2001/06/20 14:20:56 novak
* Appel a Etoile_rot::set_der0x0 dans del_deriv (au lieu de set_der0x0
* tout court).
*
* Revision 2.13 2001/03/26 09:30:58 jlz
* New members p_espec_isco and p_lspec_isco.
*
* Revision 2.12 2000/11/20 21:42:02 eric
* Appel de fait_nphi() dans le constructeur par lecture de fichier.
*
* Revision 2.11 2000/11/18 23:18:30 eric
* Modifs affichage.
*
* Revision 2.10 2000/11/18 21:09:57 eric
* Ajout des membres p_r_isco et p_f_isco.
*
* Revision 2.9 2000/11/07 16:33:08 eric
* Modif affichage.
*
* Revision 2.8 2000/10/12 15:37:01 eric
* Ajout de la fonction fait_nphi().
*
* Revision 2.7 2000/09/18 16:15:12 eric
* Ajout du membre tkij.
*
* Revision 2.6 2000/08/31 15:38:00 eric
* Bases spectrales standards pour bbb et b_car dans le constructeur
* standard (initialisation a la metrique plate).
*
* Revision 2.5 2000/08/31 11:25:45 eric
* Ajout des membres tnphi et ak_car.
*
* Revision 2.4 2000/08/25 12:28:29 eric
* Modif affichage.
*
* Revision 2.3 2000/08/18 14:01:59 eric
* Ajout de partial_display
*
* Revision 2.2 2000/08/17 12:40:04 eric
* *** empty log message ***
*
* Revision 2.1 2000/07/21 16:31:26 eric
* *** empty log message ***
*
* Revision 1.1 2000/07/20 15:32:37 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "etoile.h"
#include "eos.h"
#include "nbr_spx.h"
#include "utilitaires.h"
#include "unites.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Etoile_rot::Etoile_rot(Map& mpi, int nzet_i, bool relat, const Eos& eos_i)
: Etoile(mpi, nzet_i, relat, eos_i),
bbb(mpi),
b_car(mpi),
nphi(mpi),
tnphi(mpi),
uuu(mpi),
logn(logn_auto),
nuf(mpi),
nuq(mpi),
dzeta(beta_auto),
tggg(mpi),
w_shift(mpi, 1, CON, mp.get_bvect_cart()),
khi_shift(mpi),
tkij(mpi, 2, COV, mp.get_bvect_cart()),
ak_car(mpi),
ssjm1_nuf(mpi),
ssjm1_nuq(mpi),
ssjm1_dzeta(mpi),
ssjm1_tggg(mpi),
ssjm1_khi(mpi),
ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart())
{
// Initialization to a static state :
omega = 0 ;
uuu = 0 ;
// Initialization to a flat metric :
bbb = 1 ;
bbb.set_std_base() ;
b_car = 1 ;
b_car.set_std_base() ;
nphi = 0 ;
tnphi = 0 ;
nuf = 0 ;
nuq = 0 ;
tggg = 0 ;
w_shift.set_etat_qcq() ;
for (int i=0; i<3; i++) {
w_shift.set(i) = 0 ;
}
khi_shift.set_etat_qcq() ;
khi_shift.set() = 0 ;
tkij.set_etat_zero() ;
ak_car = 0 ;
ssjm1_nuf = 0 ;
ssjm1_nuq = 0 ;
ssjm1_dzeta = 0 ;
ssjm1_tggg = 0 ;
ssjm1_khi = 0 ;
ssjm1_wshift.set_etat_qcq() ;
for (int i=0; i<3; i++) {
ssjm1_wshift.set(i) = 0 ;
}
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
}
// Copy constructor
// ----------------
Etoile_rot::Etoile_rot(const Etoile_rot& et)
: Etoile(et),
bbb(et.bbb),
b_car(et.b_car),
nphi(et.nphi),
tnphi(et.tnphi),
uuu(et.uuu),
logn(logn_auto),
nuf(et.nuf),
nuq(et.nuq),
dzeta(beta_auto),
tggg(et.tggg),
w_shift(et.w_shift),
khi_shift(et.khi_shift),
tkij(et.tkij),
ak_car(et.ak_car),
ssjm1_nuf(et.ssjm1_nuf),
ssjm1_nuq(et.ssjm1_nuq),
ssjm1_dzeta(et.ssjm1_dzeta),
ssjm1_tggg(et.ssjm1_tggg),
ssjm1_khi(et.ssjm1_khi),
ssjm1_wshift(et.ssjm1_wshift)
{
omega = et.omega ;
// Pointers of derived quantities initialized to zero :
set_der_0x0() ;
}
// Constructor from a file
// -----------------------
Etoile_rot::Etoile_rot(Map& mpi, const Eos& eos_i, FILE* fich)
: Etoile(mpi, eos_i, fich),
bbb(mpi),
b_car(mpi),
nphi(mpi),
tnphi(mpi),
uuu(mpi),
logn(logn_auto),
nuf(mpi),
nuq(mpi),
dzeta(beta_auto),
tggg(mpi),
w_shift(mpi, 1, CON, mp.get_bvect_cart()),
khi_shift(mpi),
tkij(mpi, 2, COV, mp.get_bvect_cart()),
ak_car(mpi),
ssjm1_nuf(mpi),
ssjm1_nuq(mpi),
ssjm1_dzeta(mpi),
ssjm1_tggg(mpi),
ssjm1_khi(mpi),
ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart())
{
// Etoile parameters
// -----------------
// omega is read in the file:
fread_be(&omega, sizeof(double), 1, fich) ;
// Read of the saved fields:
// ------------------------
Tenseur nuf_file(mp, fich) ;
nuf = nuf_file ;
Tenseur nuq_file(mp, fich) ;
nuq = nuq_file ;
Tenseur tggg_file(mp, fich) ;
tggg = tggg_file ;
Tenseur w_shift_file(mp, mp.get_bvect_cart(), fich) ;
w_shift = w_shift_file ;
Tenseur khi_shift_file(mp, fich) ;
khi_shift = khi_shift_file ;
fait_shift() ; // constructs shift from w_shift and khi_shift
fait_nphi() ; // constructs N^phi from (N^x,N^y,N^z)
Cmp ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ;
ssjm1_nuf = ssjm1_nuf_file ;
Cmp ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ;
ssjm1_nuq = ssjm1_nuq_file ;
Cmp ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ;
ssjm1_dzeta = ssjm1_dzeta_file ;
Cmp ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ;
ssjm1_tggg = ssjm1_tggg_file ;
Cmp ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
ssjm1_khi = ssjm1_khi_file ;
Tenseur ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
ssjm1_wshift = ssjm1_wshift_file ;
// All other fields are initialized to zero :
// ----------------------------------------
bbb = 0 ;
b_car = 0 ;
uuu = 0 ;
// Pointers of derived quantities initialized to zero
// --------------------------------------------------
set_der_0x0() ;
}
//------------//
// Destructor //
//------------//
Etoile_rot::~Etoile_rot(){
del_deriv() ;
}
//----------------------------------//
// Management of derived quantities //
//----------------------------------//
void Etoile_rot::del_deriv() const {
Etoile::del_deriv() ;
if (p_angu_mom != 0x0) delete p_angu_mom ;
if (p_tsw != 0x0) delete p_tsw ;
if (p_grv2 != 0x0) delete p_grv2 ;
if (p_grv3 != 0x0) delete p_grv3 ;
if (p_r_circ != 0x0) delete p_r_circ ;
if (p_area != 0x0) delete p_area ;
if (p_aplat != 0x0) delete p_aplat ;
if (p_z_eqf != 0x0) delete p_z_eqf ;
if (p_z_eqb != 0x0) delete p_z_eqb ;
if (p_z_pole != 0x0) delete p_z_pole ;
if (p_mom_quad != 0x0) delete p_mom_quad ;
if (p_mom_quad_old != 0x0) delete p_mom_quad_old ;
if (p_mom_quad_Bo != 0x0) delete p_mom_quad_Bo ;
if (p_r_isco != 0x0) delete p_r_isco ;
if (p_f_isco != 0x0) delete p_f_isco ;
if (p_lspec_isco != 0x0) delete p_lspec_isco ;
if (p_espec_isco != 0x0) delete p_espec_isco ;
if (p_f_eq != 0x0) delete p_f_eq ;
Etoile_rot::set_der_0x0() ;
}
void Etoile_rot::set_der_0x0() const {
Etoile::set_der_0x0() ;
p_angu_mom = 0x0 ;
p_tsw = 0x0 ;
p_grv2 = 0x0 ;
p_grv3 = 0x0 ;
p_r_circ = 0x0 ;
p_area = 0x0 ;
p_aplat = 0x0 ;
p_z_eqf = 0x0 ;
p_z_eqb = 0x0 ;
p_z_pole = 0x0 ;
p_mom_quad = 0x0 ;
p_mom_quad_old = 0x0 ;
p_mom_quad_Bo = 0x0 ;
p_r_isco = 0x0 ;
p_f_isco = 0x0 ;
p_lspec_isco = 0x0 ;
p_espec_isco = 0x0 ;
p_f_eq = 0x0 ;
}
void Etoile_rot::del_hydro_euler() {
Etoile::del_hydro_euler() ;
del_deriv() ;
}
//--------------//
// Assignment //
//--------------//
// Assignment to another Etoile_rot
// --------------------------------
void Etoile_rot::operator=(const Etoile_rot& et) {
// Assignment of quantities common to all the derived classes of Etoile
Etoile::operator=(et) ;
// Assignement of proper quantities of class Etoile_rot
omega = et.omega ;
bbb = et.bbb ;
b_car = et.b_car ;
nphi = et.nphi ;
tnphi = et.tnphi ;
uuu = et.uuu ;
nuf = et.nuf ;
nuq = et.nuq ;
tggg = et.tggg ;
w_shift = et.w_shift ;
khi_shift = et.khi_shift ;
tkij = et.tkij ;
ak_car = et.ak_car ;
ssjm1_nuf = et.ssjm1_nuf ;
ssjm1_nuq = et.ssjm1_nuq ;
ssjm1_dzeta = et.ssjm1_dzeta ;
ssjm1_tggg = et.ssjm1_tggg ;
ssjm1_khi = et.ssjm1_khi ;
ssjm1_wshift = et.ssjm1_wshift ;
del_deriv() ; // Deletes all derived quantities
}
//--------------//
// Outputs //
//--------------//
// Save in a file
// --------------
void Etoile_rot::sauve(FILE* fich) const {
Etoile::sauve(fich) ;
fwrite_be(&omega, sizeof(double), 1, fich) ;
nuf.sauve(fich) ;
nuq.sauve(fich) ;
tggg.sauve(fich) ;
w_shift.sauve(fich) ;
khi_shift.sauve(fich) ;
ssjm1_nuf.sauve(fich) ;
ssjm1_nuq.sauve(fich) ;
ssjm1_dzeta.sauve(fich) ;
ssjm1_tggg.sauve(fich) ;
ssjm1_khi.sauve(fich) ;
ssjm1_wshift.sauve(fich) ;
}
// Printing
// --------
ostream& Etoile_rot::operator>>(ostream& ost) const {
using namespace Unites ;
Etoile::operator>>(ost) ;
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 nphi_c = nphi()(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 : " << endl ;
ost << "GRV2 = " << grv2() << endl ;
ost << "Error on the virial identity GRV3 : " << endl ;
double xgrv3 = grv3(&ost) ;
ost << "GRV3 = " << xgrv3 << endl ;
double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
/ double(1.e38) ) ;
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 ;
if (mp.get_mg()->get_np(0) == 1) {
ost << "Surface area : " << area()/(km*km) << " km^2" << endl ;
ost << "Mean radius R_mean : "
<< mean_radius()/km << " km" << endl ;
} else {
ost <<
"Skipping surface statements due to number of points in phi direction np == 1"
<< endl;
}
ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
<< endl ;
ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
ost << "Compaction parameter M_g / R_circ : " << compact << 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()(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()(0, 0, 0, 0) << endl ;
ost << "Central value of dzeta=log(AN) : "
<< dzeta()(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(0)(0, 0, 0, 0) << " "
<< w_shift(1)(0, 0, 0, 0) << " "
<< w_shift(2)(0, 0, 0, 0) << endl ;
ost << "Central value of khi_shift [km c] : "
<< khi_shift()(0, 0, 0, 0) / km << endl ;
ost << "Central value of B^2 : " << b_car()(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 ;
}
void Etoile_rot::partial_display(ostream& ost) const {
using namespace Unites ;
double omega_c = get_omega_c() ;
double freq = omega_c / (2.*M_PI) ;
ost << "Central Omega : " << omega_c * f_unit
<< " rad/s f : " << freq * f_unit << " Hz" << endl ;
ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
<< endl ;
ost << endl << "Central enthalpy : " << ent()(0,0,0,0) << " c^2" << endl ;
ost << "Central proper baryon density : " << nbar()(0,0,0,0)
<< " x 0.1 fm^-3" << endl ;
ost << "Central proper energy density : " << ener()(0,0,0,0)
<< " rho_nuc c^2" << endl ;
ost << "Central pressure : " << press()(0,0,0,0)
<< " rho_nuc c^2" << endl ;
ost << "Central value of log(N) : "
<< logn()(0, 0, 0, 0) << endl ;
ost << "Central lapse N : " << nnn()(0,0,0,0) << endl ;
ost << "Central value of dzeta=log(AN) : "
<< dzeta()(0, 0, 0, 0) << endl ;
ost << "Central value of A^2 : " << a_car()(0,0,0,0) << endl ;
ost << "Central value of B^2 : " << b_car()(0,0,0,0) << endl ;
double nphi_c = nphi()(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 ;
}
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()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
ost << endl
<< "Coordinate equatorial radius r_eq = "
<< ray_eq()/km << " km" << endl ;
ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
}
double Etoile_rot::get_omega_c() const {
return omega ;
}
// display_poly
// ------------
void Etoile_rot::display_poly(ostream& ost) const {
using namespace Unites ;
const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( &eos ) ;
if (p_eos_poly != 0x0) {
double kappa = p_eos_poly->get_kap() ;
double gamma = p_eos_poly->get_gam() ; ;
// kappa^{n/2}
double kap_ns2 = pow( kappa, 0.5 /(gamma-1) ) ;
// Polytropic unit of length in terms of r_unit :
double r_poly = kap_ns2 / sqrt(ggrav) ;
// Polytropic unit of time in terms of t_unit :
double t_poly = r_poly ;
// Polytropic unit of mass in terms of m_unit :
double m_poly = r_poly / ggrav ;
// Polytropic unit of angular momentum in terms of j_unit :
double j_poly = r_poly * r_poly / ggrav ;
// Polytropic unit of density in terms of rho_unit :
double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
ost.precision(10) ;
ost << endl << "Quantities in polytropic units : " << endl ;
ost << "==============================" << endl ;
ost << " ( r_poly = " << r_poly / km << " km )" << endl ;
ost << " n_c : " << nbar()(0, 0, 0, 0) / rho_poly << endl ;
ost << " e_c : " << ener()(0, 0, 0, 0) / rho_poly << endl ;
ost << " Omega_c : " << get_omega_c() * t_poly << endl ;
ost << " P_c : " << 2.*M_PI / get_omega_c() / t_poly << endl ;
ost << " M_bar : " << mass_b() / m_poly << endl ;
ost << " M : " << mass_g() / m_poly << endl ;
ost << " J : " << angu_mom() / j_poly << endl ;
ost << " r_eq : " << ray_eq() / r_poly << endl ;
ost << " R_circ : " << r_circ() / r_poly << endl ;
}
}
//-------------------------//
// Computational routines //
//-------------------------//
void Etoile_rot::fait_shift() {
Tenseur d_khi = khi_shift.gradient() ;
if (d_khi.get_etat() == ETATQCQ) {
d_khi.dec2_dzpuis() ; // divide by r^2 in the external compactified
// domain
}
// x_k dW^k/dx_i
Tenseur x_d_w = skxk( w_shift.gradient() ) ;
x_d_w.dec_dzpuis() ;
double lambda = double(1) / double(3) ;
// The final computation is done component by component because
// d_khi and x_d_w are covariant comp. whereas w_shift is
// contravariant
shift.set_etat_qcq() ;
for (int i=0; i<3; i++) {
shift.set(i) = (lambda+2)/2./(lambda+1) * w_shift(i)
- (lambda/2./(lambda+1)) * (d_khi(i) + x_d_w(i)) ;
}
shift.set_triad( *(w_shift.get_triad()) ) ;
}
void Etoile_rot::fait_nphi() {
if ( shift.get_etat() == ETATZERO ) {
tnphi = 0 ;
nphi = 0 ;
return ;
}
assert( shift.get_etat() == ETATQCQ ) ;
// Computation of tnphi
// --------------------
tnphi.set_etat_qcq() ;
mp.comp_p_from_cartesian( shift(0), shift(1), tnphi.set() ) ;
// Computation of nphi
// -------------------
nphi = tnphi ;
(nphi.set()).div_rsint() ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Eos_poly
Polytropic equation of state (relativistic case).
Definition eos.h:757
Lorene::Eos_poly::get_gam
double get_gam() const
Returns the adiabatic index (cf. Eq. (3)).
Definition eos_poly.C:256
Lorene::Eos_poly::get_kap
double get_kap() const
Returns the pressure coefficient (cf.
Definition eos_poly.C:260
Lorene::Eos
Equation of state base class.
Definition eos.h:190
Lorene::Etoile_rot::ssjm1_wshift
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition etoile.h:1625
Lorene::Etoile_rot::uuu
Tenseur uuu
Norm of u_euler.
Definition etoile.h:1518
Lorene::Etoile_rot::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition etoile_rot.C:335
Lorene::Etoile_rot::omega
double omega
Rotation angular velocity ([f_unit] ).
Definition etoile.h:1501
Lorene::Etoile_rot::p_z_pole
double * p_z_pole
Redshift factor at North pole.
Definition etoile.h:1640
Lorene::Etoile_rot::mom_quad
virtual double mom_quad() const
Quadrupole moment.
Definition et_rot_global.C:600
Lorene::Etoile_rot::p_z_eqf
double * p_z_eqf
Forward redshift factor at equator.
Definition etoile.h:1638
Lorene::Etoile_rot::operator=
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition etoile_rot.C:404
Lorene::Etoile_rot::r_circ
virtual double r_circ() const
Circumferential radius.
Definition et_rot_global.C:390
Lorene::Etoile_rot::set_der_0x0
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition etoile_rot.C:364
Lorene::Etoile_rot::logn
Tenseur & logn
Metric potential = logn_auto.
Definition etoile.h:1521
Lorene::Etoile_rot::p_mom_quad_old
double * p_mom_quad_old
Part of the quadrupole moment.
Definition etoile.h:1642
Lorene::Etoile_rot::ssjm1_dzeta
Cmp ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta .
Definition etoile.h:1603
Lorene::Etoile_rot::nuq
Tenseur nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition etoile.h:1531
Lorene::Etoile_rot::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile_rot.C:466
Lorene::Etoile_rot::p_z_eqb
double * p_z_eqb
Backward redshift factor at equator.
Definition etoile.h:1639
Lorene::Etoile_rot::mass_g
virtual double mass_g() const
Gravitational mass.
Definition et_rot_global.C:157
Lorene::Etoile_rot::f_isco
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition et_rot_isco.C:267
Lorene::Etoile_rot::p_r_circ
double * p_r_circ
Circumferential radius.
Definition etoile.h:1635
Lorene::Etoile_rot::khi_shift
Tenseur khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition etoile.h:1560
Lorene::Etoile_rot::aplat
virtual double aplat() const
Flatening r_pole/r_eq.
Definition et_rot_global.C:496
Lorene::Etoile_rot::grv3
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Definition et_rot_global.C:285
Lorene::Etoile_rot::Etoile_rot
Etoile_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition etoile_rot.C:137
Lorene::Etoile_rot::tsw
virtual double tsw() const
Ratio T/W.
Definition et_rot_global.C:220
Lorene::Etoile_rot::tggg
Tenseur tggg
Metric potential .
Definition etoile.h:1537
Lorene::Etoile_rot::nuf
Tenseur nuf
Part of the Metric potential = logn generated by the matter terms.
Definition etoile.h:1526
Lorene::Etoile_rot::ssjm1_tggg
Cmp ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition etoile.h:1608
Lorene::Etoile_rot::nphi
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Lorene::Etoile_rot::mass_b
virtual double mass_b() const
Baryon mass.
Definition et_rot_global.C:125
Lorene::Etoile_rot::r_isco
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition et_rot_isco.C:84
Lorene::Etoile_rot::bbb
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
Lorene::Etoile_rot::p_area
double * p_area
Surface area.
Definition etoile.h:1636
Lorene::Etoile_rot::ak_car
Tenseur ak_car
Scalar .
Definition etoile.h:1586
Lorene::Etoile_rot::~Etoile_rot
virtual ~Etoile_rot()
Destructor.
Definition etoile_rot.C:325
Lorene::Etoile_rot::dzeta
Tenseur & dzeta
Metric potential = beta_auto.
Definition etoile.h:1534
Lorene::Etoile_rot::mean_radius
virtual double mean_radius() const
Mean radius.
Definition et_rot_global.C:472
Lorene::Etoile_rot::p_grv3
double * p_grv3
Error on the virial identity GRV3.
Definition etoile.h:1634
Lorene::Etoile_rot::p_grv2
double * p_grv2
Error on the virial identity GRV2.
Definition etoile.h:1633
Lorene::Etoile_rot::p_aplat
double * p_aplat
Flatening r_pole/r_eq.
Definition etoile.h:1637
Lorene::Etoile_rot::ssjm1_nuf
Cmp ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition etoile.h:1592
Lorene::Etoile_rot::p_mom_quad_Bo
double * p_mom_quad_Bo
Part of the quadrupole moment.
Definition etoile.h:1643
Lorene::Etoile_rot::grv2
virtual double grv2() const
Error on the virial identity GRV2.
Definition et_rot_global.C:254
Lorene::Etoile_rot::fait_nphi
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition etoile_rot.C:781
Lorene::Etoile_rot::p_mom_quad
double * p_mom_quad
Quadrupole moment.
Definition etoile.h:1641
Lorene::Etoile_rot::p_f_eq
double * p_f_eq
Orbital frequency at the equator.
Definition etoile.h:1650
Lorene::Etoile_rot::del_hydro_euler
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition etoile_rot.C:389
Lorene::Etoile_rot::area
virtual double area() const
Surface area.
Definition et_rot_global.C:420
Lorene::Etoile_rot::p_angu_mom
double * p_angu_mom
Angular momentum.
Definition etoile.h:1631
Lorene::Etoile_rot::ssjm1_khi
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:1616
Lorene::Etoile_rot::angu_mom
virtual double angu_mom() const
Angular momentum.
Definition et_rot_global.C:188
Lorene::Etoile_rot::tkij
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition etoile.h:1567
Lorene::Etoile_rot::z_eqf
virtual double z_eqf() const
Forward redshift factor at equator.
Definition et_rot_global.C:513
Lorene::Etoile_rot::display_poly
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition etoile_rot.C:690
Lorene::Etoile_rot::p_f_isco
double * p_f_isco
Orbital frequency of the ISCO.
Definition etoile.h:1645
Lorene::Etoile_rot::p_lspec_isco
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition etoile.h:1649
Lorene::Etoile_rot::p_tsw
double * p_tsw
Ratio T/W.
Definition etoile.h:1632
Lorene::Etoile_rot::b_car
Tenseur b_car
Square of the metric factor B.
Definition etoile.h:1507
Lorene::Etoile_rot::z_eqb
virtual double z_eqb() const
Backward redshift factor at equator.
Definition et_rot_global.C:545
Lorene::Etoile_rot::tnphi
Tenseur tnphi
Component of the shift vector.
Definition etoile.h:1515
Lorene::Etoile_rot::p_espec_isco
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition etoile.h:1647
Lorene::Etoile_rot::p_r_isco
double * p_r_isco
Circumferential radius of the ISCO.
Definition etoile.h:1644
Lorene::Etoile_rot::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition etoile_rot.C:441
Lorene::Etoile_rot::get_omega_c
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] ).
Definition etoile_rot.C:680
Lorene::Etoile_rot::ssjm1_nuq
Cmp ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition etoile.h:1598
Lorene::Etoile_rot::z_pole
virtual double z_pole() const
Redshift factor at North pole.
Definition et_rot_global.C:580
Lorene::Etoile_rot::w_shift
Tenseur w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition etoile.h:1550
Lorene::Etoile_rot::partial_display
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition etoile_rot.C:630
Lorene::Etoile_rot::fait_shift
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition etoile_rot.C:748
Lorene::Etoile::set_der_0x0
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition etoile.C:396
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::operator=
void operator=(const Etoile &)
Assignment to another Etoile.
Definition etoile.C:430
Lorene::Etoile::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:120
Lorene::Etoile::nnn
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Lorene::Etoile::logn_auto
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:484
Lorene::Etoile::nbar
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
Lorene::Etoile::eos
const Eos & eos
Equation of state of the stellar matter.
Definition etoile.h:451
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Etoile::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition etoile.C:378
Lorene::Etoile::ener
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Lorene::Etoile::Etoile
Etoile(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition etoile.C:143
Lorene::Etoile::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile.C:511
Lorene::Etoile::press
Tenseur press
Fluid pressure.
Definition etoile.h:461
Lorene::Etoile::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition etoile.C:483
Lorene::Etoile::shift
Tenseur shift
Total shift vector.
Definition etoile.h:512
Lorene::Etoile::ent
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case).
Definition etoile.h:457
Lorene::Etoile::del_hydro_euler
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition etoile.C:410
Lorene::Etoile::beta_auto
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:506
Lorene::Etoile::a_car
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::Tenseur::dec2_dzpuis
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
Lorene::Tenseur::dec_dzpuis
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
Lorene::Tenseur::get_etat
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
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::fread_be
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
Lorene::fwrite_be
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.
Definition fwrite_be.C:70
Lorene::skxk
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Definition tenseur_operateur.C:580
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::get_bvect_cart
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:139
Unites
Standard units of space, time and mass.