etoile_global.C
/*
* Methods of class Etoile to compute global quantities
*/
/*
* 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_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_global.C,v 1.7 2014/10/13 08:52:59 j_novak Exp $" ;
/*
* $Id: etoile_global.C,v 1.7 2014/10/13 08:52:59 j_novak Exp $
* $Log: etoile_global.C,v $
* Revision 1.7 2014/10/13 08:52:59 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2012/08/12 17:48:36 p_cerda
* Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
*
* Revision 1.5 2005/01/18 22:37:30 k_taniguchi
* Modify the method of ray_eq(int kk).
*
* Revision 1.4 2005/01/18 20:35:46 k_taniguchi
* Addition of ray_eq(int kk).
*
* Revision 1.3 2005/01/17 20:40:56 k_taniguchi
* Addition of ray_eq_3pis2().
*
* Revision 1.2 2003/12/05 14:50:26 j_novak
* To suppress some warnings...
*
* Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
* LORENE
*
* Revision 1.2 2000/07/21 12:02:14 eric
* Etoile::ray_eq_pi() : traitement du cas type_p = SYM.
*
* Revision 1.1 2000/01/28 17:18:45 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_global.C,v 1.7 2014/10/13 08:52:59 j_novak Exp $
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "etoile.h"
//--------------------------//
// Stellar surface //
//--------------------------//
namespace Lorene {
const Itbl& Etoile::l_surf() const {
if (p_l_surf == 0x0) { // a new computation is required
assert(p_xi_surf == 0x0) ; // consistency check
int np = mp.get_mg()->get_np(0) ;
int nt = mp.get_mg()->get_nt(0) ;
p_l_surf = new Itbl(np, nt) ;
p_xi_surf = new Tbl(np, nt) ;
double ent0 = 0 ; // definition of the surface
double precis = 1.e-15 ;
int nitermax = 100 ;
int niter ;
(ent().va).equipot(ent0, nzet, precis, nitermax, niter, *p_l_surf,
*p_xi_surf) ;
}
return *p_l_surf ;
}
const Tbl& Etoile::xi_surf() const {
if (p_xi_surf == 0x0) { // a new computation is required
assert(p_l_surf == 0x0) ; // consistency check
l_surf() ; // the computation is done by l_surf()
}
return *p_xi_surf ;
}
//--------------------------//
// Coordinate radii //
//--------------------------//
double Etoile::ray_eq() const {
if (p_ray_eq == 0x0) { // a new computation is required
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
#ifndef NDEBUG
int type_p = mg.get_type_p() ;
#endif
int nt = mg.get_nt(0) ;
if ( type_t == SYM ) {
assert( (type_p == SYM) || (type_p == NONSYM) ) ;
int k = 0 ;
int j = nt-1 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
double theta = M_PI / 2 ;
double phi = 0 ;
p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
}
else {
assert( (type_p == SYM) || (type_p == NONSYM) ) ;
int k = 0 ;
int j = (nt-1)/2 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
double theta = M_PI / 2 ;
double phi = 0 ;
p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
// cout << "Etoile::ray_eq : the case type_t = " << type_t
// << " is not contemplated yet !" << endl ;
//abort() ;
}
}
return *p_ray_eq ;
}
double Etoile::ray_eq_pis2() const {
if (p_ray_eq_pis2 == 0x0) { // a new computation is required
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
int type_p = mg.get_type_p() ;
int nt = mg.get_nt(0) ;
int np = mg.get_np(0) ;
if ( type_t == SYM ) {
int j = nt-1 ;
double theta = M_PI / 2 ;
double phi = M_PI / 2 ;
switch (type_p) {
case SYM : {
int k = np / 2 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
case NONSYM : {
assert( np % 4 == 0 ) ;
int k = np / 4 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_pis2 : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
else {
int j = (nt-1)/2 ;
double theta = M_PI / 2 ;
double phi = M_PI / 2 ;
switch (type_p) {
case SYM : {
int k = np / 2 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
case NONSYM : {
assert( np % 4 == 0 ) ;
int k = np / 4 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_pis2 : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
}
return *p_ray_eq_pis2 ;
}
double Etoile::ray_eq_pi() const {
if (p_ray_eq_pi == 0x0) { // a new computation is required
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
int type_p = mg.get_type_p() ;
int nt = mg.get_nt(0) ;
int np = mg.get_np(0) ;
if ( type_t == SYM ) {
switch (type_p) {
case SYM : {
p_ray_eq_pi = new double( ray_eq() ) ;
break ;
}
case NONSYM : {
int k = np / 2 ;
int j = nt-1 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
double theta = M_PI / 2 ;
double phi = M_PI ;
p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_pi : the case type_t = " << type_t
<< " and type_p = " << type_p << endl ;
cout << " is not contemplated yet !" << endl ;
abort() ;
break ;
}
}
}else{
switch (type_p) {
case SYM : {
p_ray_eq_pi = new double( ray_eq() ) ;
break ;
}
case NONSYM : {
int k = np / 2 ;
int j = (nt-1)/2 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
double theta = M_PI / 2 ;
double phi = M_PI ;
p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_pi : the case type_t = " << type_t
<< " and type_p = " << type_p << endl ;
cout << " is not contemplated yet !" << endl ;
abort() ;
break ;
}
}
}
}
return *p_ray_eq_pi ;
}
double Etoile::ray_eq_3pis2() const {
if (p_ray_eq_3pis2 == 0x0) { // a new computation is required
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
int type_p = mg.get_type_p() ;
int nt = mg.get_nt(0) ;
int np = mg.get_np(0) ;
if ( type_t == SYM ) {
int j = nt-1 ;
double theta = M_PI / 2 ;
double phi = 3. * M_PI / 2 ;
switch (type_p) {
case SYM : {
p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
break ;
}
case NONSYM : {
assert( np % 4 == 0 ) ;
int k = 3 * np / 4 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_3pis2 = new double( mp.val_r(l,xi,theta,phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_3pis2 : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
else {
int j = (nt-1)/2 ;
double theta = M_PI / 2 ;
double phi = 3. * M_PI / 2 ;
switch (type_p) {
case SYM : {
p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
break ;
}
case NONSYM : {
assert( np % 4 == 0 ) ;
int k = 3 * np / 4 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
p_ray_eq_3pis2 = new double( mp.val_r(l,xi,theta,phi) ) ;
break ;
}
default : {
cout << "Etoile::ray_eq_3pis2 : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
}
return *p_ray_eq_3pis2 ;
}
double Etoile::ray_pole() const {
if (p_ray_pole == 0x0) { // a new computation is required
#ifndef NDEBUG
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
#endif
assert( (type_t == SYM) || (type_t == NONSYM) ) ;
int k = 0 ;
int j = 0 ;
int l = l_surf()(k, j) ;
double xi = xi_surf()(k, j) ;
double theta = 0 ;
double phi = 0 ;
p_ray_pole = new double( mp.val_r(l, xi, theta, phi) ) ;
}
return *p_ray_pole ;
}
double Etoile::ray_eq(int kk) const {
const Mg3d& mg = *(mp.get_mg()) ;
int type_t = mg.get_type_t() ;
int type_p = mg.get_type_p() ;
int nt = mg.get_nt(0) ;
int np = mg.get_np(0) ;
assert( kk >= 0 ) ;
assert( kk < np ) ;
double ray_eq_kk ;
if ( type_t == SYM ) {
int j = nt-1 ;
double theta = M_PI / 2 ;
switch (type_p) {
case SYM : {
cout << "Etoile::ray_eq(kk) : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
case NONSYM : {
double phi = 2. * kk * M_PI / np ;
int l = l_surf()(kk, j) ;
double xi = xi_surf()(kk, j) ;
ray_eq_kk = mp.val_r(l,xi,theta,phi) ;
break ;
}
default : {
cout << "Etoile::ray_eq(kk) : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
else {
int j = (nt-1)/2 ;
double theta = M_PI / 2 ;
switch (type_p) {
case SYM : {
cout << "Etoile::ray_eq(kk) : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
case NONSYM : {
double phi = 2. * kk * M_PI / np ;
int l = l_surf()(kk, j) ;
double xi = xi_surf()(kk, j) ;
ray_eq_kk = mp.val_r(l,xi,theta,phi) ;
break ;
}
default : {
cout << "Etoile::ray_eq(kk) : the case type_p = "
<< type_p << " is not contemplated yet !" << endl ;
abort() ;
}
}
}
return ray_eq_kk ;
}
//--------------------------//
// Baryon mass //
//--------------------------//
double Etoile::mass_b() const {
if (p_mass_b == 0x0) { // a new computation is required
if (relativistic) {
cout <<
"Etoile::mass_b : in the relativistic case, the baryon mass"
<< endl <<
"computation cannot be performed by the base class Etoile !"
<< endl ;
abort() ;
}
assert(nbar.get_etat() == ETATQCQ) ;
p_mass_b = new double( nbar().integrale() ) ;
}
return *p_mass_b ;
}
//----------------------------//
// Gravitational mass //
//----------------------------//
double Etoile::mass_g() const {
if (p_mass_g == 0x0) { // a new computation is required
if (relativistic) {
cout <<
"Etoile::mass_g : in the relativistic case, the gravitational mass"
<< endl <<
"computation cannot be performed by the base class Etoile !"
<< endl ;
abort() ;
}
p_mass_g = new double( mass_b() ) ; // in the Newtonian case
// M_g = M_b
}
return *p_mass_g ;
}
}
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::p_mass_b
double * p_mass_b
Baryon mass.
Definition etoile.h:547
Lorene::Etoile::p_mass_g
double * p_mass_g
Gravitational mass.
Definition etoile.h:548
Lorene::Etoile::p_xi_surf
Tbl * p_xi_surf
Description of the stellar surface: 2-D Tbl containing the values of the radial coordinate on the su...
Definition etoile.h:545
Lorene::Etoile::ray_eq_pi
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:256
Lorene::Etoile::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:120
Lorene::Etoile::nbar
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
Lorene::Etoile::p_l_surf
Itbl * p_l_surf
Description of the stellar surface: 2-D Itbl containing the values of the domain index l on the surfa...
Definition etoile.h:539
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Etoile::p_ray_eq_3pis2
double * p_ray_eq_3pis2
Coordinate radius at , .
Definition etoile.h:530
Lorene::Etoile::p_ray_eq_pi
double * p_ray_eq_pi
Coordinate radius at , .
Definition etoile.h:527
Lorene::Etoile::mass_b
virtual double mass_b() const
Baryon mass.
Definition etoile_global.C:525
Lorene::Etoile::ray_eq_pis2
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:169
Lorene::Etoile::p_ray_pole
double * p_ray_pole
Coordinate radius at .
Definition etoile.h:533
Lorene::Etoile::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 ...
Definition etoile_global.C:75
Lorene::Etoile::ray_eq_3pis2
double ray_eq_3pis2() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:335
Lorene::Etoile::mass_g
virtual double mass_g() const
Gravitational mass.
Definition etoile_global.C:551
Lorene::Etoile::relativistic
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Lorene::Etoile::p_ray_eq_pis2
double * p_ray_eq_pis2
Coordinate radius at , .
Definition etoile.h:524
Lorene::Etoile::ent
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case).
Definition etoile.h:457
Lorene::Etoile::xi_surf
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinat...
Definition etoile_global.C:101
Lorene::Etoile::p_ray_eq
double * p_ray_eq
Coordinate radius at , .
Definition etoile.h:521
Lorene::Etoile::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition etoile_global.C:415
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
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_type_p
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
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::Tbl
Basic array class.
Definition tbl.h:161
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Lorene::integrale
virtual Tbl * integrale(const Cmp &) const =0
Computes the integral over all space of a Cmp .