strot_dirac_global.C
/*
* Methods for computing global quantities within the class Star_rot_Dirac
*
* (see file star.h for documentation).
*
*/
/*
* Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
*
* 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 strot_dirac_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $" ;
/*
* $Id: strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
* $Log: strot_dirac_global.C,v $
* Revision 1.13 2014/10/13 08:53:40 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.12 2014/10/06 15:13:18 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.11 2010/11/17 11:21:52 j_novak
* Corrected minor problems in the case nz=1 and nt=1.
*
* Revision 1.10 2010/10/22 08:08:40 j_novak
* Removal of the method Star_rot_dirac::lambda_grv2() and call to the C++ version of integrale2d.
*
* Revision 1.9 2009/10/26 10:54:33 j_novak
* Added the case of a NONSYM base in theta.
*
* Revision 1.8 2008/05/30 08:27:38 j_novak
* New global quantities rp_circ and ellipt (circumferential polar coordinate and
* ellipticity).
*
* Revision 1.7 2008/02/18 13:51:20 j_novak
* Change of a dzpuis
*
* Revision 1.6 2005/03/25 14:11:49 j_novak
* In the 1D case, GRV2 returns -1 (because of a problem in integral2d).
*
* Revision 1.5 2005/02/17 17:30:42 f_limousin
* Change the name of some quantities to be consistent with other classes
* (for instance nnn is changed to nn, shift to beta, beta to lnq...)
*
* Revision 1.4 2005/02/09 13:36:01 lm_lin
*
* Add the calculations of GRV2, T/W, R_circ, and flattening.
*
* Revision 1.3 2005/02/02 10:11:24 j_novak
* Better calculation of the GRV3 identity.
*
*
*
* $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_global.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
*
*/
// C headers
#include <cmath>
#include <cassert>
// Lorene headers
#include "star_rot_dirac.h"
#include "unites.h"
#include "utilitaires.h"
#include "proto.h"
//-------------------------------------------------------//
// Baryonic mass //
// //
// Note: In Lorene units, mean particle mass is unity //
//-------------------------------------------------------//
namespace Lorene {
double Star_rot_Dirac::mass_b() const {
if (p_mass_b == 0x0) { // a new computation is required
Scalar dens = sqrt( gamma.determinant() ) * gam_euler * nbar ;
dens.std_spectral_base() ;
p_mass_b = new double( dens.integrale() ) ;
}
return *p_mass_b ;
}
//---------------------------------------------------------//
// Gravitational mass //
// //
// Note: This is the Komar mass for stationary and //
// asymptotically flat spacetime (see, eg, Wald) //
//---------------------------------------------------------//
double Star_rot_Dirac::mass_g() const {
if (p_mass_g == 0x0) { // a new computation is required
Scalar j_source = 2.* contract(contract(gamma.cov(), 0, j_euler, 0),
0, beta, 0) ;
Scalar dens = sqrt( gamma.determinant() ) *
( nn * (ener_euler + s_euler) - j_source ) ;
dens.std_spectral_base() ;
p_mass_g = new double( dens.integrale() ) ;
}
return *p_mass_g ;
}
//--------------------------------------//
// Angular momentum //
// //
// Komar-type integral (see, eg, Wald) //
//--------------------------------------//
double Star_rot_Dirac::angu_mom() const {
if (p_angu_mom == 0x0) { // a new computation is required
// phi_kill = axial killing vector
Vector phi_kill(mp, CON, mp.get_bvect_spher()) ;
phi_kill.set(1).set_etat_zero() ;
phi_kill.set(2).set_etat_zero() ;
phi_kill.set(3) = 1. ;
phi_kill.set(3).std_spectral_base() ;
phi_kill.set(3).mult_rsint() ;
Scalar j_source = contract(contract(gamma.cov(), 0, j_euler, 0),
0, phi_kill, 0) ;
Scalar dens = sqrt( gamma.determinant() ) * j_source ;
dens.std_spectral_base() ;
p_angu_mom = new double( dens.integrale() ) ;
}
return *p_angu_mom ;
}
//---------------------//
// T/W //
//---------------------//
double Star_rot_Dirac::tsw() const {
if (p_tsw == 0x0) { // a new computation is required
double tcin = 0.5 * omega * angu_mom() ;
Scalar dens = sqrt( gamma.determinant() ) * gam_euler * ener ;
dens.std_spectral_base() ;
double mass_p = dens.integrale() ;
p_tsw = new double( tcin / ( mass_p + tcin - mass_g() ) ) ;
}
return *p_tsw ;
}
//--------------------------------------------------------------//
// GRV2 //
// cf. Eq. (28) of Bonazzola & Gourgoulhon CQG, 11, 1775 (1994) //
// //
//--------------------------------------------------------------//
double Star_rot_Dirac::grv2() const {
using namespace Unites ;
if (p_grv2 == 0x0) { // a new computation is required
// determinant of the 2-metric k_{ab}
Scalar k_det = gamma.cov()(1,1)*gamma.cov()(2,2) -
gamma.cov()(1,2)*gamma.cov()(1,2) ;
//**
// sou_m = 8\pi T_{\mu\nu} m^{\mu}m^{\nu}
// => sou_m = 8\pi [ (E+P) U^2 + P ], where v^2 = v_i v^i
//
//**
Scalar sou_m = 2 * qpig * ( (ener_euler + press)*v2 + press ) ;
sou_m = sqrt( k_det )*sou_m ;
sou_m.std_spectral_base() ;
// This is the term 3k_a k^a.
Scalar sou_q = 3 *( taa(1,3) * aa(1,3) + taa(2,3)*aa(2,3) ) ;
// This is the term \nu_{|| a}\nu^{|| a}.
//
Scalar sou_tmp = gamma.con()(1,1) * logn.dsdr() * logn.dsdr() ;
Scalar term_2 = 2 * gamma.con()(1,2) * logn.dsdr() * logn.dsdt() ;
term_2.div_r_dzpuis(4) ;
Scalar term_3 = gamma.con()(2,2) * logn.dsdt() * logn.dsdt() ;
term_3.div_r_dzpuis(2) ;
term_3.div_r_dzpuis(4) ;
sou_tmp += term_2 + term_3 ;
// Source of the gravitational part
sou_q -= sou_tmp ;
sou_q = sqrt( k_det )*sou_q ;
sou_q.std_spectral_base() ;
p_grv2 = new double( double(1) + integrale2d(sou_m)/integrale2d(sou_q) ) ;
}
return *p_grv2 ;
}
//-------------------------------------------------------------//
// GRV3 //
// cf. Eq. (29) of Gourgoulhon & Bonazzola CQG, 11, 443 (1994) //
//-------------------------------------------------------------//
double Star_rot_Dirac::grv3() const {
using namespace Unites ;
if (p_grv3 == 0x0) { // a new computation is required
// Gravitational term
// -------------------
Scalar sou_q = 0.75*aa_quad - contract(logn.derive_cov(gamma), 0,
logn.derive_con(gamma), 0) ;
Tensor t_tmp = contract(gamma.connect().get_delta(), 2,
gamma.connect().get_delta(), 0) ;
Scalar tmp_1 = 0.25* contract( gamma.con(), 0, 1,
contract(t_tmp, 0, 3), 0, 1 ) ;
Scalar tmp_2 = 0.25* contract( gamma.con(), 0, 1,
contract( contract( gamma.connect().get_delta(), 0, 1),
0, gamma.connect().get_delta(), 0), 0, 1) ;
sou_q = sou_q + tmp_1 - tmp_2 ;
sou_q = sqrt( gamma.determinant() ) * sou_q ;
sou_q.std_spectral_base() ;
double int_grav = sou_q.integrale() ;
// Matter term
// --------------
Scalar sou_m = qpig*s_euler ;
sou_m = sqrt( gamma.determinant() ) * sou_m ;
sou_m.std_spectral_base() ;
double int_mat = sou_m.integrale() ;
p_grv3 = new double( (int_grav + int_mat) / int_mat ) ;
}
return *p_grv3 ;
}
//--------------------//
// R_circ //
//--------------------//
double Star_rot_Dirac::r_circ() const {
if (p_r_circ ==0x0) { // a new computation is required
// Index of the point at phi=0, theta=pi/2 at the surface of the star:
const Mg3d* mg = mp.get_mg() ;
int l_b = nzet - 1 ;
int i_b = mg->get_nr(l_b) - 1 ;
int j_b = (mg->get_type_t() == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b) / 2) ;
int k_b = 0 ;
double gamma_phi = gamma.cov()(3,3).val_grid_point(l_b, k_b, j_b, i_b) ;
p_r_circ = new double( sqrt( gamma_phi ) * ray_eq() ) ;
}
return *p_r_circ ;
}
//--------------------//
// R_circ //
//--------------------//
double Star_rot_Dirac::rp_circ() const {
if (p_rp_circ ==0x0) { // a new computation is required
const Mg3d& mg = *mp.get_mg() ;
int nz = mg.get_nzone() ;
assert(nz>2) ;
int np = mg.get_np(0) ;
if (np != 1) {
cout << "The polar circumferential radius is only well defined\n"
<< "with np = 1!" << endl ;
abort() ;
}
int nt = mg.get_nt(0) ;
Sym_tensor gam = gamma.cov() ;
const Coord& tet = mp.tet ;
Scalar rrr(mp) ;
rrr.annule_hard() ;
rrr.annule(nzet,nz-1) ;
double phi = 0 ;
for (int j=0; j<nt; j++) {
double theta = (+tet)(0, 0, j, 0) ;
double erre =
mp.val_r(l_surf()(0,j), xi_surf()(0, j), theta, phi) ;
for (int lz=0; lz<nzet; lz++) {
int nrz = mg.get_nr(lz) ;
for (int i=0; i<nrz; i++) {
rrr.set_grid_point(lz, 0, j, i) = erre ;
}
}
}
rrr.std_spectral_base() ;
Scalar drrr = rrr.dsdt() ;
Scalar fff(mp) ;
fff.annule_hard() ;
fff.annule(nzet,nz-1) ;
for (int j=0; j<nt; j++) {
double theta = (+tet)(0, 0, j, 0) ;
int ls = l_surf()(0, j) ;
double xs = xi_surf()(0, j) ;
double grr = gam(1,1).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
double grt = gam(1,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
double gtt = gam(2,2).get_spectral_va().val_point_jk(ls, xs, j, 0) ;
double rr = mp.val_r(ls, xs, theta, phi) ;
double dr = drrr.get_spectral_va().val_point_jk(ls, xs, j, 0) ;
for (int i=0; i< mg.get_nr(nzet-1); i++) {
fff.set_grid_point(nzet-1, 0, j, i)
= sqrt(grr*dr*dr + 2*grt*rr*dr + gtt*rr*rr) ;
}
}
fff.std_spectral_base() ;
fff.set_spectral_va().coef() ;
p_rp_circ = new double((*fff.get_spectral_va().c_cf)(nzet-1, 0, 0, 0)) ;
}
return *p_rp_circ ;
}
//--------------------------//
// Flattening //
//--------------------------//
double Star_rot_Dirac::aplat() const {
return ray_pole() / ray_eq() ;
}
//--------------------------//
// Ellipticity //
//--------------------------//
double Star_rot_Dirac::ellipt() const {
return sqrt(1. - (rp_circ()*rp_circ())/(r_circ()*r_circ())) ;
}
}
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_type_t
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition grilles.h:485
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::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition scalar_r_manip.C:147
Lorene::Scalar::dsdt
const Scalar & dsdt() const
Returns of *this .
Definition scalar_deriv.C:208
Lorene::Scalar::integrale
double integrale() const
Computes the integral over all space of *this .
Definition scalar_integ.C:61
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_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:604
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version).
Definition scalar.h:601
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
Lorene::Scalar::mult_rsint
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Definition scalar_r_manip.C:281
Lorene::Scalar::set_grid_point
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
Lorene::Star_rot_Dirac::aa
Sym_tensor aa
Definition star_rot_dirac.h:87
Lorene::Star_rot_Dirac::mass_g
virtual double mass_g() const
Gravitational mass.
Definition strot_dirac_global.C:117
Lorene::Star_rot_Dirac::ellipt
virtual double ellipt() const
Ellipticity e.
Definition strot_dirac_global.C:428
Lorene::Star_rot_Dirac::grv3
virtual double grv3() const
Error on the virial identity GRV3.
Definition strot_dirac_global.C:273
Lorene::Star_rot_Dirac::angu_mom
virtual double angu_mom() const
Angular momentum.
Definition strot_dirac_global.C:143
Lorene::Star_rot_Dirac::p_grv3
double * p_grv3
Error on the virial identity GRV3.
Definition star_rot_dirac.h:108
Lorene::Star_rot_Dirac::r_circ
virtual double r_circ() const
Circumferential equatorial radius.
Definition strot_dirac_global.C:331
Lorene::Star_rot_Dirac::omega
double omega
Rotation angular velocity ([f_unit] ).
Definition star_rot_dirac.h:59
Lorene::Star_rot_Dirac::p_tsw
double * p_tsw
Ratio T/W.
Definition star_rot_dirac.h:109
Lorene::Star_rot_Dirac::p_r_circ
double * p_r_circ
Circumferential equatorial radius.
Definition star_rot_dirac.h:110
Lorene::Star_rot_Dirac::tsw
virtual double tsw() const
Ratio T/W.
Definition strot_dirac_global.C:177
Lorene::Star_rot_Dirac::j_euler
Vector j_euler
Momentum density 3-vector with respect to the Eulerian observer.
Definition star_rot_dirac.h:77
Lorene::Star_rot_Dirac::taa
Sym_tensor taa
Definition star_rot_dirac.h:88
Lorene::Star_rot_Dirac::rp_circ
virtual double rp_circ() const
Circumferential polar radius.
Definition strot_dirac_global.C:357
Lorene::Star_rot_Dirac::v2
Scalar v2
Definition star_rot_dirac.h:78
Lorene::Star_rot_Dirac::p_grv2
double * p_grv2
Error on the virial identity GRV2.
Definition star_rot_dirac.h:107
Lorene::Star_rot_Dirac::mass_b
virtual double mass_b() const
Baryonic mass.
Definition strot_dirac_global.C:94
Lorene::Star_rot_Dirac::aa_quad
Scalar aa_quad
Definition star_rot_dirac.h:89
Lorene::Star_rot_Dirac::p_angu_mom
double * p_angu_mom
Angular momentum.
Definition star_rot_dirac.h:106
Lorene::Star_rot_Dirac::p_rp_circ
double * p_rp_circ
Circumferential polar radius.
Definition star_rot_dirac.h:111
Lorene::Star_rot_Dirac::aplat
virtual double aplat() const
Flattening r_pole/r_eq.
Definition strot_dirac_global.C:418
Lorene::Star_rot_Dirac::grv2
virtual double grv2() const
Error on the virial identity GRV2.
Definition strot_dirac_global.C:204
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::l_surf
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l on...
Definition star_global.C:63
Lorene::Star::p_mass_b
double * p_mass_b
Baryon mass.
Definition star.h:268
Lorene::Star::nbar
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
Lorene::Star::ener_euler
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
Lorene::Star::xi_surf
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinate ...
Definition star_global.C:89
Lorene::Star::gam_euler
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
Lorene::Star::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition star_global.C:108
Lorene::Star::s_euler
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Lorene::Star::p_mass_g
double * p_mass_g
Gravitational mass.
Definition star.h:269
Lorene::Star::gamma
Metric gamma
3-metric
Definition star.h:235
Lorene::Star::press
Scalar press
Fluid pressure.
Definition star.h:194
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::beta
Vector beta
Shift vector.
Definition star.h:228
Lorene::Star::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition star_global.C:278
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Valeur::val_point_jk
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Definition valeur.C:900
Lorene::Valeur::c_cf
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
Lorene::Valeur::coef
void coef() const
Computes the coeffcients of *this.
Definition valeur_coef.C:148
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Lorene::tet
Coord tet
coordinate centered on the grid
Definition map.h:719
Unites
Standard units of space, time and mass.