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
30char 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// --------------------
81namespace Lorene {
82Et_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),
112 omega_field(et.omega_field),
113 omega_min(et.omega_min),
114 omega_max(et.omega_max),
115 prim_field(et.prim_field)
116 {}
117
118
119// Constructor from a file
120// -----------------------
121Et_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
160
161 // Assignment of proper quantities of class Etoile_rot
162 frot = et.frot ;
163 primfrot = et.primfrot ;
164 par_frot = et.par_frot ;
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
179void 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
196ostream& Et_rot_diff::operator>>(ostream& ost) const {
197
198 using namespace Unites ;
199
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
268void Et_rot_diff::display_poly(ostream& ost) const {
269
270 using namespace Unites ;
271
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
315double Et_rot_diff::funct_omega(double omeg) const {
316
317 return frot(omeg, par_frot) ;
318
319}
320
321double Et_rot_diff::prim_funct_omega(double omeg) const {
322
323 return primfrot(omeg, par_frot) ;
324
325}
326
328
329 return omega_field()(0, 0, 0, 0) ;
330
331}
332
333
334
335}
Polytropic equation of state (relativistic case).
Definition eos.h:757
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition eos_poly.C:256
double get_kap() const
Returns the pressure coefficient (cf.
Definition eos_poly.C:260
Equation of state base class.
Definition eos.h:190
Class for differentially rotating stars.
Definition et_rot_diff.h:70
double omega_min
Minimum value of .
Tenseur prim_field
Field .
virtual void display_poly(ostream &) const
Display in polytropic units.
double omega_max
Maximum value of .
virtual double tsw() const
Ratio T/W.
Tenseur omega_field
Field .
void fait_prim_field()
Computes the member prim_field from omga_field .
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
virtual void sauve(FILE *) const
Save in a file.
virtual ~Et_rot_diff()
Destructor.
Tbl par_frot
Parameters of the function .
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition et_rot_diff.h:93
double prim_funct_omega(double omeg) const
Evaluates the primitive of , where F is the function defining the rotation profile.
void operator=(const Et_rot_diff &)
Assignment to another Et_rot_diff.
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
double(* frot)(double, const Tbl &)
Function defining the rotation profile.
Definition et_rot_diff.h:85
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition etoile.h:1496
double omega
Rotation angular velocity ([f_unit] )
Definition etoile.h:1501
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition etoile_rot.C:404
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile_rot.C:466
virtual double mass_g() const
Gravitational mass.
virtual double aplat() const
Flatening r_pole/r_eq.
virtual double mass_b() const
Baryon mass.
virtual double angu_mom() const
Angular momentum.
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition etoile_rot.C:690
virtual void sauve(FILE *) const
Save in a file.
Definition etoile_rot.C:441
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
const Eos & eos
Equation of state of the stellar matter.
Definition etoile.h:451
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Tenseur press
Fluid pressure.
Definition etoile.h:461
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Basic array class.
Definition tbl.h:161
void sauve(FILE *) const
Save in a file.
Definition tbl.C:326
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
void sauve(FILE *) const
Save in a file.
Definition tenseur.C:1325
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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 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 prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.