LORENE
et_rot_diff.C
1 /*
2  * Methods for class Et_rot_diff.
3  *
4  * (see file et_rot_diff.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2001 Eric Gourgoulhon
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 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 $" ;
31 
32 /*
33  * $Id: et_rot_diff.C,v 1.4 2014/10/13 08:52:57 j_novak Exp $
34  * $Log: et_rot_diff.C,v $
35  * Revision 1.4 2014/10/13 08:52:57 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.3 2004/03/25 10:29:05 j_novak
39  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
40  *
41  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
42  *
43  * All writing/reading to a binary file are now performed according to
44  * the big endian convention, whatever the system is big endian or
45  * small endian, thanks to the functions fwrite_be and fread_be
46  *
47  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
48  * LORENE
49  *
50  * Revision 1.3 2001/10/25 09:20:54 eric
51  * Ajout de la fonction virtuelle display_poly.
52  * Affichage de Max nbar, ener et press.
53  *
54  * Revision 1.2 2001/10/24 16:23:01 eric
55  * *** empty log message ***
56  *
57  * Revision 1.1 2001/10/19 08:18:10 eric
58  * Initial revision
59  *
60  *
61  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff.C,v 1.4 2014/10/13 08:52:57 j_novak Exp $
62  *
63  */
64 
65 // Headers C
66 #include "math.h"
67 
68 // Headers Lorene
69 #include "et_rot_diff.h"
70 #include "eos.h"
71 #include "nbr_spx.h"
72 #include "utilitaires.h"
73 #include "unites.h"
74 
75  //--------------//
76  // Constructors //
77  //--------------//
78 
79 // Standard constructor
80 // --------------------
81 namespace Lorene {
82 Et_rot_diff::Et_rot_diff(Map& mp_i, int nzet_i, bool relat, const Eos& eos_i,
83  double (*frot_i)(double, const Tbl&),
84  double (*primfrot_i)(double, const Tbl&),
85  const Tbl& par_frot_i)
86  : Etoile_rot(mp_i, nzet_i, relat, eos_i),
87  frot(frot_i),
88  primfrot(primfrot_i),
89  par_frot(par_frot_i),
90  omega_field(mp_i),
91  prim_field(mp_i)
92  {
93 
94  // To make sure that omega is not used
95  omega = __infinity ;
96 
97  // Initialization to a static state :
98  omega_field = 0 ;
99  prim_field = 0 ;
100  omega_min = 0 ;
101  omega_max = 0 ;
102 
103 }
104 
105 // Copy constructor
106 // ----------------
108  : Etoile_rot(et),
109  frot(et.frot),
110  primfrot(et.primfrot),
111  par_frot(et.par_frot),
113  omega_min(et.omega_min),
114  omega_max(et.omega_max),
116  {}
117 
118 
119 // Constructor from a file
120 // -----------------------
121 Et_rot_diff::Et_rot_diff(Map& mp_i, const Eos& eos_i, FILE* fich,
122  double (*frot_i)(double, const Tbl&),
123  double (*primfrot_i)(double, const Tbl&) )
124  : Etoile_rot(mp_i, eos_i, fich),
125  frot(frot_i),
126  primfrot(primfrot_i),
127  par_frot(fich),
128  omega_field(mp_i),
129  prim_field(mp_i)
130  {
131 
132  Tenseur omega_field_file(mp, fich) ;
133  omega_field = omega_field_file ;
134  fait_prim_field() ;
135 
136  // omega_min and omega_max are read in the file:
137  fread_be(&omega_min, sizeof(double), 1, fich) ;
138  fread_be(&omega_max, sizeof(double), 1, fich) ;
139 
140 }
141 
142 
143  //------------//
144  // Destructor //
145  //------------//
146 
148 
149 
150  //--------------//
151  // Assignment //
152  //--------------//
153 
154 // Assignment to another Et_rot_diff
155 // ---------------------------------
157 
158  // Assignment of quantities common to all the derived classes of Etoile_rot
159  Etoile_rot::operator=(et) ;
160 
161  // Assignment of proper quantities of class Etoile_rot
162  frot = et.frot ;
163  primfrot = et.primfrot ;
164  par_frot = et.par_frot ;
165  omega_field = et.omega_field ;
166  prim_field = et.prim_field ;
167  omega_min = et.omega_min ;
168  omega_max = et.omega_max ;
169 
170 }
171 
172  //--------------//
173  // Outputs //
174  //--------------//
175 
176 // Save in a file
177 // --------------
178 
179 void Et_rot_diff::sauve(FILE* fich) const {
180 
181  Etoile_rot::sauve(fich) ;
182 
183  par_frot.sauve(fich) ;
184 
185  omega_field.sauve(fich) ;
186 
187  fwrite_be(&omega_min, sizeof(double), 1, fich) ;
188  fwrite_be(&omega_max, sizeof(double), 1, fich) ;
189 
190 }
191 
192 
193 // Printing
194 // --------
195 
196 ostream& Et_rot_diff::operator>>(ostream& ost) const {
197 
198  using namespace Unites ;
199 
200  Etoile_rot::operator>>(ost) ;
201 
202  ost << endl ;
203  ost << "Differentially rotating star" << endl ;
204  ost << "-----------------------------" << endl ;
205 
206  ost << endl << "Parameters of F(Omega) : " << endl ;
207  ost << par_frot << endl ;
208 
209  ost << "Min, Max of Omega/(2pi) : " << omega_min / (2*M_PI) * f_unit
210  << " Hz, " << omega_max / (2*M_PI) * f_unit << " Hz" << endl ;
211  int lsurf = nzet - 1;
212  int nt = mp.get_mg()->get_nt(lsurf) ;
213  int nr = mp.get_mg()->get_nr(lsurf) ;
214  ost << "Central, equatorial value of Omega: "
215  << omega_field()(0, 0, 0, 0) * f_unit << " rad/s, "
216  << omega_field()(nzet-1, 0, nt-1, nr-1) * f_unit << " rad/s" << endl ;
217 
218  ost << "Central, equatorial value of Omega/(2 Pi): "
219  << omega_field()(0, 0, 0, 0) * f_unit / (2*M_PI) << " Hz, "
220  << omega_field()(nzet-1, 0, nt-1, nr-1) * f_unit / (2*M_PI)
221  << " Hz" << endl ;
222 
223  double nbar_max = max( max( nbar() ) ) ;
224  double ener_max = max( max( ener() ) ) ;
225  double press_max = max( max( press() ) ) ;
226  ost << "Max prop. bar. dens. : " << nbar_max
227  << " x 0.1 fm^-3 = " << nbar_max / nbar()(0, 0, 0, 0) << " central"
228  << endl ;
229  ost << "Max prop. ener. dens. (e_max) : " << ener_max
230  << " rho_nuc c^2 = " << ener_max / ener()(0, 0, 0, 0) << " central"
231  << endl ;
232  ost << "Max pressure : " << press_max
233  << " rho_nuc c^2 = " << press_max / press()(0, 0, 0, 0) << " central"
234  << endl ;
235 
236  // Length scale set by the maximum energy density:
237  double len_rho = 1. / sqrt( ggrav * ener_max ) ;
238  ost << endl << "Value of A = par_frot(1) in units of e_max^{1/2} : " <<
239  par_frot(1) / len_rho << endl ;
240 
241  ost << "Value of A / r_eq : " <<
242  par_frot(1) / ray_eq() << endl ;
243 
244  ost << "r_p/r_eq : " << aplat() << endl ;
245  ost << "KEH l^2 = (c/G^2)^{2/3} J^2 e_max^{1/3} M_B^{-10/3} : " <<
246  angu_mom() * angu_mom() / pow(len_rho, 0.6666666666666666)
247  / pow(mass_b(), 3.3333333333333333)
248  / pow(ggrav, 1.3333333333333333) << endl ;
249 
250  ost << "M e_max^{1/2} : " << ggrav * mass_g() / len_rho << endl ;
251 
252  ost << "r_eq e_max^{1/2} : " << ray_eq() / len_rho << endl ;
253 
254  ost << "T/W : " << tsw() << endl ;
255 
256  ost << "Omega_c / e_max^{1/2} : " << get_omega_c() * len_rho << endl ;
257 
258  display_poly(ost) ;
259 
260  return ost ;
261 
262 
263 }
264 
265 // display_poly
266 // ------------
267 
268 void Et_rot_diff::display_poly(ostream& ost) const {
269 
270  using namespace Unites ;
271 
272  Etoile_rot::display_poly( ost ) ;
273 
274  const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( &eos ) ;
275 
276  if (p_eos_poly != 0x0) {
277 
278  double kappa = p_eos_poly->get_kap() ;
279  double gamma = p_eos_poly->get_gam() ; ;
280 
281  // kappa^{n/2}
282  double kap_ns2 = pow( kappa, 0.5 /(gamma-1) ) ;
283 
284  // Polytropic unit of length in terms of r_unit :
285  double r_poly = kap_ns2 / sqrt(ggrav) ;
286 
287  // Polytropic unit of time in terms of t_unit :
288  double t_poly = r_poly ;
289 
290  // Polytropic unit of density in terms of rho_unit :
291  double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
292 
293  ost.precision(10) ;
294  ost << " n_max : " << max( max( nbar() ) ) / rho_poly << endl ;
295  ost << " e_max : " << max( max( ener() ) ) / rho_poly << endl ;
296  ost << " P_min : " << 2.*M_PI / omega_max / t_poly << endl ;
297  ost << " P_max : " << 2.*M_PI / omega_min / t_poly << endl ;
298 
299  int lsurf = nzet - 1;
300  int nt = mp.get_mg()->get_nt(lsurf) ;
301  int nr = mp.get_mg()->get_nr(lsurf) ;
302  ost << " P_eq : " << 2.*M_PI /
303  omega_field()(nzet-1, 0, nt-1, nr-1) / t_poly << endl ;
304 
305  }
306 
307 }
308 
309 
310 
311  //-----------------------//
312  // Miscellaneous //
313  //-----------------------//
314 
315 double Et_rot_diff::funct_omega(double omeg) const {
316 
317  return frot(omeg, par_frot) ;
318 
319 }
320 
321 double Et_rot_diff::prim_funct_omega(double omeg) const {
322 
323  return primfrot(omeg, par_frot) ;
324 
325 }
326 
327 double Et_rot_diff::get_omega_c() const {
328 
329  return omega_field()(0, 0, 0, 0) ;
330 
331 }
332 
333 
334 
335 }
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
void operator=(const Et_rot_diff &)
Assignment to another Et_rot_diff.
Definition: et_rot_diff.C:156
Tenseur prim_field
Field .
Definition: et_rot_diff.h:111
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
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition: etoile_rot.C:690
Lorene prototypes.
Definition: app_hor.h:64
double(* frot)(double, const Tbl &)
Function defining the rotation profile.
Definition: et_rot_diff.h:85
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:190
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double ray_eq() const
Coordinate radius at , [r_unit].
Base class for coordinate mappings.
Definition: map.h:670
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition: etoile.h:1496
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition: et_rot_diff.h:93
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition: etoile_rot.C:404
Tenseur press
Fluid pressure.
Definition: etoile.h:461
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
Definition: et_rot_diff.C:315
double omega_max
Maximum value of .
Definition: et_rot_diff.h:108
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:256
Class for differentially rotating stars.
Definition: et_rot_diff.h:70
void fait_prim_field()
Computes the member prim_field from omga_field .
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition: et_rot_diff.C:268
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile_rot.C:466
virtual double angu_mom() const
Angular momentum.
void sauve(FILE *) const
Save in a file.
Definition: tenseur.C:1325
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:260
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_rot.C:441
Tbl par_frot
Parameters of the function .
Definition: et_rot_diff.h:102
virtual double mass_b() const
Baryon mass.
Polytropic equation of state (relativistic case).
Definition: eos.h:757
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
const Eos & eos
Equation of state of the stellar matter.
Definition: etoile.h:451
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
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
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
double omega
Rotation angular velocity ([f_unit] )
Definition: etoile.h:1501
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
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
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
virtual void sauve(FILE *) const
Save in a file.
Definition: et_rot_diff.C:179
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:460
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: et_rot_diff.C:196
void sauve(FILE *) const
Save in a file.
Definition: tbl.C:326
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
double omega_min
Minimum value of .
Definition: et_rot_diff.h:107
virtual double mass_g() const
Gravitational mass.
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
Definition: et_rot_diff.C:327
virtual ~Et_rot_diff()
Destructor.
Definition: et_rot_diff.C:147
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
Tenseur omega_field
Field .
Definition: et_rot_diff.h:105
virtual double aplat() const
Flatening r_pole/r_eq.
virtual double tsw() const
Ratio T/W.