LORENE
star_rot_dirac_diff.C
1 /*
2  * Methods of class Star_rot_Dirac_diff
3  *
4  * (see file star_rot_dirac_diff.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Motoyuki Saijo
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 version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Header: /cvsroot/Lorene/C++/Source/Star/star_rot_dirac_diff.C,v 1.5 2016/12/05 16:18:15 j_novak Exp $
32  *
33  */
34 
35 
36 // C headers
37 #include <cmath>
38 #include <cassert>
39 
40 // Lorene headers
41 #include "star_rot_dirac_diff.h"
42 #include "unites.h"
43 #include "utilitaires.h"
44 
45 
46  //--------------//
47  // Constructors //
48  //--------------//
49 
50 // Standard constructor
51 //-------------------------
52 namespace Lorene {
53 Star_rot_Dirac_diff::Star_rot_Dirac_diff(Map& mpi, int nzet_i, const Eos& eos_i,
54  double (*frot_i)(double, const Tbl&),
55  double (*primfrot_i)(double, const Tbl&),
56  const Tbl& par_frot_i)
57  :Star_rot_Dirac(mpi, nzet_i, eos_i),
58  frot(frot_i),
59  primfrot(primfrot_i),
60  par_frot(par_frot_i),
61  omega_field(mpi),
62  prim_field(mpi)
63  {
64 
65  // Initialization to a static state :
66  omega_field = 0 ;
67  prim_field = 0 ;
68  omega_min = 0 ;
69  omega_max = 0 ;
70 
71 }
72 
73 
74 // Copy constructor
75 //-----------------
77  : Star_rot_Dirac(star),
78  frot(star.frot),
79  primfrot(star.primfrot),
80  par_frot(star.par_frot),
81  omega_field(star.omega_field),
82  omega_min(star.omega_min),
83  omega_max(star.omega_max),
84  prim_field(star.prim_field)
85  {}
86 
87 
88 //Constructor from a file //## to be more general...
89 //------------------------
90 Star_rot_Dirac_diff::Star_rot_Dirac_diff(Map& mpi, const Eos& eos_i, FILE* fich,
91  double (*frot_i)(double, const Tbl&),
92  double (*primfrot_i)(double, const Tbl&) )
93  : Star_rot_Dirac(mpi, eos_i, fich),
94  frot(frot_i),
95  primfrot(primfrot_i),
96  par_frot(fich),
97  omega_field(mpi),
98  prim_field(mpi)
99  {
100 
101  Scalar omega_field_file(mp, *(mp.get_mg()), fich) ;
102  omega_field = omega_field_file ;
103  fait_prim_field() ;
104 
105  // omega_min and omega_max are read in the file:
106  fread_be(&omega_min, sizeof(double), 1, fich) ;
107  fread_be(&omega_max, sizeof(double), 1, fich) ;
108 
109 }
110 
111 
112  //------------//
113  // Destructor //
114  //------------//
115 
117 
118 
119  //---------------//
120  // Assignment //
121  //---------------//
122 
123 // Assignment to another Star_rot_Dirac_diff
124 // ------------------------------------
125 
127 
128  // Assignment of quantities common to all the derived classes of
129  // Star_rot_Dirac
131 
132  // Assignment of proper quantities of class Star_rot_Dirac
133  frot = star.frot ;
134  primfrot = star.primfrot ;
135  par_frot = star.par_frot ;
136  omega_field = star.omega_field ;
137  prim_field = star.prim_field ;
138  omega_min = star.omega_min ;
139  omega_max = star.omega_max ;
140 
141 }
142 
143 
144  //-----------//
145  // Outputs //
146  //-----------//
147 
148 // Save in a file
149 // --------------
150 
151 void Star_rot_Dirac_diff::sauve(FILE* fich) const {
152 
153  Star::sauve(fich) ;
154 
155  par_frot.sauve(fich) ;
156 
157  omega_field.sauve(fich) ;
158 
159  fwrite_be(&omega_min, sizeof(double), 1, fich) ;
160  fwrite_be(&omega_max, sizeof(double), 1, fich) ;
161 
162  // What else to save? //## to be more general ...
163 
164 }
165 
166 
167 // Printing
168 // ---------
169 
170 ostream& Star_rot_Dirac_diff::operator>>(ostream& ost) const {
171 
172  using namespace Unites ;
173 
174  Star::operator>>(ost) ;
175 
176  ost << "Differentially rotating star in Dirac gauge" << '\n';
177 
178  // Only differentially rotating star for the moment....
179  ost << '\n';
180  ost << "Differentially rotating star" << '\n';
181  ost << "-----------------------" << '\n';
182 
183  ost << '\n'<< "Parameters of F(Omega) : " << '\n';
184  ost << par_frot << '\n';
185 
186  ost << "Min, Max of Omega/(2pi) : " << omega_min / (2*M_PI) * f_unit
187  << " Hz, " << omega_max / (2*M_PI) * f_unit << " Hz" << '\n';
188  int lsurf = nzet - 1;
189  int nt = mp.get_mg()->get_nt(lsurf) ;
190  int nr = mp.get_mg()->get_nr(lsurf) ;
191  ost << "Central, equatorial value of Omega: "
192  << omega_field.val_grid_point(0, 0, 0, 0) * f_unit
193  << " rad/s, "
194  << omega_field.val_grid_point(nzet-1, 0, nt-1, nr-1) * f_unit
195  << " rad/s" << '\n';
196 
197  ost << "Central, equatorial value of Omega/(2 Pi): "
198  << omega_field.val_grid_point(0, 0, 0, 0) * f_unit / (2*M_PI)
199  << " Hz, "
200  << omega_field.val_grid_point(nzet-1, 0, nt-1, nr-1) *
201  f_unit / (2*M_PI)
202  << " Hz" << '\n';
203 
204  ost << "Error on the virial identity GRV2 : " << '\n';
205  ost << "GRV2 = " << grv2() << '\n';
206  ost << "Error on the virial identity GRV3 : " << '\n';
207  ost << "GRV3 = " << grv3() << '\n';
208 
209  ost << "Angular momentum J : "
210  << angu_mom()/( qpig / (4*M_PI) *msol*msol) << " G M_sol^2 / c"
211  << '\n';
212  ost << "c J / (G M^2) : "
213  << angu_mom()/( qpig / (4*M_PI) * pow(mass_g(), 2.) ) << '\n';
214 
215  if (omega != 0.) {
216  double mom_iner = angu_mom() / omega ;
217  double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
218  / double(1.e38) ) ;
219  ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
220  << '\n';
221  }
222 
223  ost << "Ratio T/W : " << tsw() << '\n';
224 
225 // ost << "Omega_c / e_max^{1/2} : " << get_omega_c() * len_rho << '\n';
226 
227  ost << "Circumferential equatorial radius R_circ : "
228  << r_circ()/km << " km" << '\n';
229  ost << "Circumferential polar radius Rp_circ : "
230  << rp_circ()/km << " km" << endl ;
231  ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
232  << endl ;
233  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
234  ost << "Ellipticity sqrt(1-(Rp_circ/R_circ)^2) : " << ellipt() << endl ;
235  double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
236  ost << "Compaction parameter M_g / R_circ : " << compact << '\n';
237 
238 
239  // More to come here.....
240 
241  return ost ;
242 
243 }
244 
245 
246  //-----------------------//
247  // Miscellaneous //
248  //-----------------------//
249 
250 double Star_rot_Dirac_diff::funct_omega(double omeg) const {
251 
252  return frot(omeg, par_frot) ;
253 
254 }
255 
256 double Star_rot_Dirac_diff::prim_funct_omega(double omeg) const {
257 
258  return primfrot(omeg, par_frot) ;
259 
260 }
261 
263 
264  return omega_field.val_grid_point(0, 0, 0, 0) ;
265 
266 }
267 }
void fait_prim_field()
Computes the member prim_field from omega_field .
virtual double grv3() const
Error on the virial identity GRV3.
virtual double tsw() const
Ratio T/W.
void operator=(const Star_rot_Dirac_diff &)
Assignment to another Star_rot_Dirac_diff.
double omega
Rotation angular velocity ([f_unit] )
Map & mp
Mapping associated with the star.
Definition: star.h:180
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:209
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
virtual double r_circ() const
Circumferential equatorial radius.
void operator=(const Star_rot_Dirac &)
Assignment to another Star_rot_Dirac.
virtual void sauve(FILE *) const
Save in a file.
Tbl par_frot
Parameters of the function .
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
virtual double angu_mom() const
Angular momentum.
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:692
virtual double mass_g() const
Gravitational mass.
virtual double ellipt() const
Ellipticity e.
virtual double grv2() const
Error on the virial identity GRV2.
double omega_max
Maximum value of .
double(* frot)(double, const Tbl &)
Function defining the rotation profile.
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:73
Star_rot_Dirac_diff(Map &mp_i, int nzet_i, const Eos &eos_i, double(*frot_i)(double, const Tbl &), double(*primfrot_i)(double, const Tbl &), const Tbl &par_frot_i)
Standard constructor.
Class for relativistic rotating stars in Dirac gauge and maximal slicing.
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
double prim_funct_omega(double omeg) const
Evaluates the primitive of , where F is the function defining the rotation profile.
virtual double aplat() const
Flattening r_pole/r_eq.
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:72
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual ~Star_rot_Dirac_diff()
Destructor.
Class for relativistic differentially rotating stars in Dirac gauge and maximal slicing.
virtual double rp_circ() const
Circumferential polar radius.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: star.C:420
void sauve(FILE *) const
Save in a file.
Definition: tbl.C:329
virtual void sauve(FILE *) const
Save in a file.
Definition: star.C:400
double omega_min
Minimum value of .
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
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.