et_rot_diff.C
/*
* Methods for class Et_rot_diff.
*
* (see file et_rot_diff.h for documentation).
*
*/
/*
* Copyright (c) 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 et_rot_diff_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff.C,v 1.4 2014/10/13 08:52:57 j_novak Exp $" ;
/*
* $Id: et_rot_diff.C,v 1.4 2014/10/13 08:52:57 j_novak Exp $
* $Log: et_rot_diff.C,v $
* Revision 1.4 2014/10/13 08:52:57 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.3 2004/03/25 10:29:05 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* 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 1.3 2001/10/25 09:20:54 eric
* Ajout de la fonction virtuelle display_poly.
* Affichage de Max nbar, ener et press.
*
* Revision 1.2 2001/10/24 16:23:01 eric
* *** empty log message ***
*
* Revision 1.1 2001/10/19 08:18:10 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff.C,v 1.4 2014/10/13 08:52:57 j_novak Exp $
*
*/
// Headers C
#include "math.h"
// Headers Lorene
#include "et_rot_diff.h"
#include "eos.h"
#include "nbr_spx.h"
#include "utilitaires.h"
#include "unites.h"
//--------------//
// Constructors //
//--------------//
// Standard constructor
// --------------------
namespace Lorene {
Et_rot_diff::Et_rot_diff(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i,
double (*frot_i)(double, const Tbl&),
double (*primfrot_i)(double, const Tbl&),
const Tbl& par_frot_i)
: Etoile_rot(mp_i, nzet_i, relat, eos_i),
frot(frot_i),
primfrot(primfrot_i),
par_frot(par_frot_i),
omega_field(mp_i),
prim_field(mp_i)
{
// To make sure that omega is not used
omega = __infinity ;
// Initialization to a static state :
omega_field = 0 ;
prim_field = 0 ;
omega_min = 0 ;
omega_max = 0 ;
}
// Copy constructor
// ----------------
Et_rot_diff::Et_rot_diff(const Et_rot_diff& et)
: Etoile_rot(et),
frot(et.frot),
primfrot(et.primfrot),
par_frot(et.par_frot),
omega_field(et.omega_field),
omega_min(et.omega_min),
omega_max(et.omega_max),
prim_field(et.prim_field)
{}
// Constructor from a file
// -----------------------
Et_rot_diff::Et_rot_diff(Map& mp_i, const Eos& eos_i, FILE* fich,
double (*frot_i)(double, const Tbl&),
double (*primfrot_i)(double, const Tbl&) )
: Etoile_rot(mp_i, eos_i, fich),
frot(frot_i),
primfrot(primfrot_i),
par_frot(fich),
omega_field(mp_i),
prim_field(mp_i)
{
Tenseur omega_field_file(mp, fich) ;
omega_field = omega_field_file ;
fait_prim_field() ;
// omega_min and omega_max are read in the file:
fread_be(&omega_min, sizeof(double), 1, fich) ;
fread_be(&omega_max, sizeof(double), 1, fich) ;
}
//------------//
// Destructor //
//------------//
Et_rot_diff::~Et_rot_diff(){}
//--------------//
// Assignment //
//--------------//
// Assignment to another Et_rot_diff
// ---------------------------------
void Et_rot_diff::operator=(const Et_rot_diff& et) {
// Assignment of quantities common to all the derived classes of Etoile_rot
Etoile_rot::operator=(et) ;
// Assignment of proper quantities of class Etoile_rot
frot = et.frot ;
primfrot = et.primfrot ;
par_frot = et.par_frot ;
omega_field = et.omega_field ;
prim_field = et.prim_field ;
omega_min = et.omega_min ;
omega_max = et.omega_max ;
}
//--------------//
// Outputs //
//--------------//
// Save in a file
// --------------
void Et_rot_diff::sauve(FILE* fich) const {
Etoile_rot::sauve(fich) ;
par_frot.sauve(fich) ;
omega_field.sauve(fich) ;
fwrite_be(&omega_min, sizeof(double), 1, fich) ;
fwrite_be(&omega_max, sizeof(double), 1, fich) ;
}
// Printing
// --------
ostream& Et_rot_diff::operator>>(ostream& ost) const {
using namespace Unites ;
Etoile_rot::operator>>(ost) ;
ost << endl ;
ost << "Differentially rotating star" << endl ;
ost << "-----------------------------" << endl ;
ost << endl << "Parameters of F(Omega) : " << endl ;
ost << par_frot << endl ;
ost << "Min, Max of Omega/(2pi) : " << omega_min / (2*M_PI) * f_unit
<< " Hz, " << omega_max / (2*M_PI) * f_unit << " Hz" << endl ;
int lsurf = nzet - 1;
int nt = mp.get_mg()->get_nt(lsurf) ;
int nr = mp.get_mg()->get_nr(lsurf) ;
ost << "Central, equatorial value of Omega: "
<< omega_field()(0, 0, 0, 0) * f_unit << " rad/s, "
<< omega_field()(nzet-1, 0, nt-1, nr-1) * f_unit << " rad/s" << endl ;
ost << "Central, equatorial value of Omega/(2 Pi): "
<< omega_field()(0, 0, 0, 0) * f_unit / (2*M_PI) << " Hz, "
<< omega_field()(nzet-1, 0, nt-1, nr-1) * f_unit / (2*M_PI)
<< " Hz" << endl ;
double nbar_max = max( max( nbar() ) ) ;
double ener_max = max( max( ener() ) ) ;
double press_max = max( max( press() ) ) ;
ost << "Max prop. bar. dens. : " << nbar_max
<< " x 0.1 fm^-3 = " << nbar_max / nbar()(0, 0, 0, 0) << " central"
<< endl ;
ost << "Max prop. ener. dens. (e_max) : " << ener_max
<< " rho_nuc c^2 = " << ener_max / ener()(0, 0, 0, 0) << " central"
<< endl ;
ost << "Max pressure : " << press_max
<< " rho_nuc c^2 = " << press_max / press()(0, 0, 0, 0) << " central"
<< endl ;
// Length scale set by the maximum energy density:
double len_rho = 1. / sqrt( ggrav * ener_max ) ;
ost << endl << "Value of A = par_frot(1) in units of e_max^{1/2} : " <<
par_frot(1) / len_rho << endl ;
ost << "Value of A / r_eq : " <<
par_frot(1) / ray_eq() << endl ;
ost << "r_p/r_eq : " << aplat() << endl ;
ost << "KEH l^2 = (c/G^2)^{2/3} J^2 e_max^{1/3} M_B^{-10/3} : " <<
angu_mom() * angu_mom() / pow(len_rho, 0.6666666666666666)
/ pow(mass_b(), 3.3333333333333333)
/ pow(ggrav, 1.3333333333333333) << endl ;
ost << "M e_max^{1/2} : " << ggrav * mass_g() / len_rho << endl ;
ost << "r_eq e_max^{1/2} : " << ray_eq() / len_rho << endl ;
ost << "T/W : " << tsw() << endl ;
ost << "Omega_c / e_max^{1/2} : " << get_omega_c() * len_rho << endl ;
display_poly(ost) ;
return ost ;
}
// display_poly
// ------------
void Et_rot_diff::display_poly(ostream& ost) const {
using namespace Unites ;
Etoile_rot::display_poly( ost ) ;
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 density in terms of rho_unit :
double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
ost.precision(10) ;
ost << " n_max : " << max( max( nbar() ) ) / rho_poly << endl ;
ost << " e_max : " << max( max( ener() ) ) / rho_poly << endl ;
ost << " P_min : " << 2.*M_PI / omega_max / t_poly << endl ;
ost << " P_max : " << 2.*M_PI / omega_min / t_poly << endl ;
int lsurf = nzet - 1;
int nt = mp.get_mg()->get_nt(lsurf) ;
int nr = mp.get_mg()->get_nr(lsurf) ;
ost << " P_eq : " << 2.*M_PI /
omega_field()(nzet-1, 0, nt-1, nr-1) / t_poly << endl ;
}
}
//-----------------------//
// Miscellaneous //
//-----------------------//
double Et_rot_diff::funct_omega(double omeg) const {
return frot(omeg, par_frot) ;
}
double Et_rot_diff::prim_funct_omega(double omeg) const {
return primfrot(omeg, par_frot) ;
}
double Et_rot_diff::get_omega_c() const {
return omega_field()(0, 0, 0, 0) ;
}
}
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::Et_rot_diff::omega_min
double omega_min
Minimum value of .
Definition et_rot_diff.h:107
Lorene::Et_rot_diff::prim_field
Tenseur prim_field
Field .
Definition et_rot_diff.h:111
Lorene::Et_rot_diff::display_poly
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition et_rot_diff.C:268
Lorene::Et_rot_diff::omega_max
double omega_max
Maximum value of .
Definition et_rot_diff.h:108
Lorene::Et_rot_diff::tsw
virtual double tsw() const
Ratio T/W.
Definition et_rot_diff_global.C:67
Lorene::Et_rot_diff::omega_field
Tenseur omega_field
Field .
Definition et_rot_diff.h:105
Lorene::Et_rot_diff::fait_prim_field
void fait_prim_field()
Computes the member prim_field from omga_field .
Definition et_rot_diff_faitomeg.C:199
Lorene::Et_rot_diff::Et_rot_diff
Et_rot_diff(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i, double(*frot_i)(double, const Tbl &), double(*primfrot_i)(double, const Tbl &), const Tbl &par_frot_i)
Standard constructor.
Definition et_rot_diff.C:82
Lorene::Et_rot_diff::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition et_rot_diff.C:179
Lorene::Et_rot_diff::~Et_rot_diff
virtual ~Et_rot_diff()
Destructor.
Definition et_rot_diff.C:147
Lorene::Et_rot_diff::par_frot
Tbl par_frot
Parameters of the function .
Definition et_rot_diff.h:102
Lorene::Et_rot_diff::primfrot
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition et_rot_diff.h:93
Lorene::Et_rot_diff::prim_funct_omega
double prim_funct_omega(double omeg) const
Evaluates the primitive of , where F is the function defining the rotation profile.
Definition et_rot_diff.C:321
Lorene::Et_rot_diff::operator=
void operator=(const Et_rot_diff &)
Assignment to another Et_rot_diff.
Definition et_rot_diff.C:156
Lorene::Et_rot_diff::get_omega_c
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] ).
Definition et_rot_diff.C:327
Lorene::Et_rot_diff::funct_omega
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
Definition et_rot_diff.C:315
Lorene::Et_rot_diff::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition et_rot_diff.C:196
Lorene::Et_rot_diff::frot
double(* frot)(double, const Tbl &)
Function defining the rotation profile.
Definition et_rot_diff.h:85
Lorene::Etoile_rot::omega
double omega
Rotation angular velocity ([f_unit] ).
Definition etoile.h:1501
Lorene::Etoile_rot::operator=
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition etoile_rot.C:404
Lorene::Etoile_rot::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile_rot.C:466
Lorene::Etoile_rot::mass_g
virtual double mass_g() const
Gravitational mass.
Definition et_rot_global.C:157
Lorene::Etoile_rot::aplat
virtual double aplat() const
Flatening r_pole/r_eq.
Definition et_rot_global.C:496
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::mass_b
virtual double mass_b() const
Baryon mass.
Definition et_rot_global.C:125
Lorene::Etoile_rot::angu_mom
virtual double angu_mom() const
Angular momentum.
Definition et_rot_global.C:188
Lorene::Etoile_rot::display_poly
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition etoile_rot.C:690
Lorene::Etoile_rot::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition etoile_rot.C:441
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
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::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::ener
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Lorene::Etoile::press
Tenseur press
Fluid pressure.
Definition etoile.h:461
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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
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.