LORENE
star_global.C
1 /*
2  * Methods of class Star to compute global quantities
3  */
4 
5 /*
6  * Copyright (c) 2004 francois Limousin
7  *
8  * This file is part of LORENE.
9  *
10  * LORENE is free software; you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation; either version 2 of the License, or
13  * (at your option) any later version.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 
27 
28 
29 /*
30  * $Id: star_global.C,v 1.7 2016/12/05 16:18:15 j_novak Exp $
31  * $Log: star_global.C,v $
32  * Revision 1.7 2016/12/05 16:18:15 j_novak
33  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34  *
35  * Revision 1.6 2014/10/13 08:53:39 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.5 2013/04/25 15:46:06 j_novak
39  * Added special treatment in the case np = 1, for type_p = NONSYM.
40  *
41  * Revision 1.4 2009/10/26 10:54:33 j_novak
42  * Added the case of a NONSYM base in theta.
43  *
44  * Revision 1.3 2007/06/21 19:55:09 k_taniguchi
45  * Introduction of a method to compute ray_eq_3pis2().
46  *
47  * Revision 1.2 2004/01/20 15:20:48 f_limousin
48  * First version
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Star/star_global.C,v 1.7 2016/12/05 16:18:15 j_novak Exp $
52  *
53  */
54 
55 // Headers C
56 #include "math.h"
57 
58 // Headers Lorene
59 #include "star.h"
60 
61  //--------------------------//
62  // Stellar surface //
63  //--------------------------//
64 
65 namespace Lorene {
66 const Itbl& Star::l_surf() const {
67 
68  if (p_l_surf == 0x0) { // a new computation is required
69 
70  assert(p_xi_surf == 0x0) ; // consistency check
71 
72  int np = mp.get_mg()->get_np(0) ;
73  int nt = mp.get_mg()->get_nt(0) ;
74 
75  p_l_surf = new Itbl(np, nt) ;
76  p_xi_surf = new Tbl(np, nt) ;
77 
78  double ent0 = 0 ; // definition of the surface
79  double precis = 1.e-15 ;
80  int nitermax = 100 ;
81  int niter ;
82 
83  (ent.get_spectral_va()).equipot(ent0, nzet, precis, nitermax, niter, *p_l_surf,
84  *p_xi_surf) ;
85 
86  }
87 
88  return *p_l_surf ;
89 
90 }
91 
92 const Tbl& Star::xi_surf() const {
93 
94  if (p_xi_surf == 0x0) { // a new computation is required
95 
96  assert(p_l_surf == 0x0) ; // consistency check
97 
98  l_surf() ; // the computation is done by l_surf()
99 
100  }
101 
102  return *p_xi_surf ;
103 
104 }
105 
106 
107  //--------------------------//
108  // Coordinate radii //
109  //--------------------------//
110 
111 double Star::ray_eq() const {
112 
113  if (p_ray_eq == 0x0) { // a new computation is required
114 
115  const Mg3d& mg = *(mp.get_mg()) ;
116 
117  int type_t = mg.get_type_t() ;
118 #ifndef NDEBUG
119  int type_p = mg.get_type_p() ;
120 #endif
121  int nt = mg.get_nt(0) ;
122 
123  assert( (type_t == SYM) || (type_t == NONSYM) ) ;
124  assert( (type_p == SYM) || (type_p == NONSYM) ) ;
125  int k = 0 ;
126  int j = (type_t == SYM ? nt-1 : nt / 2);
127  int l = l_surf()(k, j) ;
128  double xi = xi_surf()(k, j) ;
129  double theta = M_PI / 2 ;
130  double phi = 0 ;
131 
132  p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
133 
134  }
135 
136  return *p_ray_eq ;
137 
138 }
139 
140 
141 double Star::ray_eq_pis2() const {
142 
143  if (p_ray_eq_pis2 == 0x0) { // a new computation is required
144 
145  const Mg3d& mg = *(mp.get_mg()) ;
146 
147  int type_t = mg.get_type_t() ;
148  int type_p = mg.get_type_p() ;
149  int nt = mg.get_nt(0) ;
150  int np = mg.get_np(0) ;
151 
152  int j = (type_t == SYM ? nt-1 : nt / 2);
153  double theta = M_PI / 2 ;
154  double phi = M_PI / 2 ;
155 
156  switch (type_p) {
157 
158  case SYM : {
159  int k = np / 2 ;
160  int l = l_surf()(k, j) ;
161  double xi = xi_surf()(k, j) ;
162  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
163  break ;
164  }
165 
166  case NONSYM : {
167  assert( (np == 1) || (np % 4 == 0) ) ;
168  int k = np / 4 ;
169  int l = l_surf()(k, j) ;
170  double xi = xi_surf()(k, j) ;
171  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
172  break ;
173  }
174 
175  default : {
176  cout << "Star::ray_eq_pis2 : the case type_p = "
177  << type_p << " is not contemplated yet !" << endl ;
178  abort() ;
179  }
180  }
181 
182  }
183 
184  return *p_ray_eq_pis2 ;
185 
186 }
187 
188 
189 double Star::ray_eq_pi() const {
190 
191  if (p_ray_eq_pi == 0x0) { // a new computation is required
192 
193  const Mg3d& mg = *(mp.get_mg()) ;
194 
195  int type_t = mg.get_type_t() ;
196  int type_p = mg.get_type_p() ;
197  int nt = mg.get_nt(0) ;
198  int np = mg.get_np(0) ;
199 
200  assert ( ( type_t == SYM ) || ( type_t == NONSYM ) ) ;
201 
202  switch (type_p) {
203 
204  case SYM : {
205  p_ray_eq_pi = new double( ray_eq() ) ;
206  break ;
207  }
208 
209  case NONSYM : {
210  int k = np / 2 ;
211  int j = (type_t == SYM ? nt-1 : nt/2 ) ;
212  int l = l_surf()(k, j) ;
213  double xi = xi_surf()(k, j) ;
214  double theta = M_PI / 2 ;
215  double phi = M_PI ;
216 
217  p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
218  break ;
219  }
220 
221  default : {
222 
223  cout << "Star::ray_eq_pi : the case type_p = " << type_p << endl ;
224  cout << " is not contemplated yet !" << endl ;
225  abort() ;
226  break ;
227  }
228  }
229 
230  }
231 
232  return *p_ray_eq_pi ;
233 
234 }
235 
236 double Star::ray_eq_3pis2() const {
237 
238  if (p_ray_eq_3pis2 == 0x0) { // a new computation is required
239 
240  const Mg3d& mg = *(mp.get_mg()) ;
241 
242  int type_t = mg.get_type_t() ;
243  int type_p = mg.get_type_p() ;
244  int nt = mg.get_nt(0) ;
245  int np = mg.get_np(0) ;
246 
247  assert( ( type_t == SYM ) || ( type_t == NONSYM ) ) ;
248 
249  int j = (type_t == SYM ? nt-1 : nt/2 );
250  double theta = M_PI / 2 ;
251  double phi = 3. * M_PI / 2 ;
252 
253  switch (type_p) {
254 
255  case SYM : {
256  p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
257  break ;
258  }
259 
260  case NONSYM : {
261  assert( np % 4 == 0 ) ;
262  int k = 3 * np / 4 ;
263  int l = l_surf()(k, j) ;
264  double xi = xi_surf()(k, j) ;
265  p_ray_eq_3pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
266  break ;
267  }
268 
269  default : {
270  cout << "Star::ray_eq_3pis2 : the case type_p = "
271  << type_p << " is not implemented yet !" << endl ;
272  abort() ;
273  }
274  }
275  }
276 
277  return *p_ray_eq_3pis2 ;
278 
279 }
280 
281 double Star::ray_pole() const {
282 
283  if (p_ray_pole == 0x0) { // a new computation is required
284 
285 #ifndef NDEBUG
286  const Mg3d& mg = *(mp.get_mg()) ;
287  int type_t = mg.get_type_t() ;
288 #endif
289  assert( (type_t == SYM) || (type_t == NONSYM) ) ;
290 
291  int k = 0 ;
292  int j = 0 ;
293  int l = l_surf()(k, j) ;
294  double xi = xi_surf()(k, j) ;
295  double theta = 0 ;
296  double phi = 0 ;
297 
298  p_ray_pole = new double( mp.val_r(l, xi, theta, phi) ) ;
299 
300  }
301 
302  return *p_ray_pole ;
303 
304 }
305 
306  //--------------------------//
307  // Baryon mass //
308  //--------------------------//
309 
310 double Star::mass_b() const {
311 
312  if (p_mass_b == 0x0) { // a new computation is required
313 
314  cout <<
315  "Star::mass_b : in the relativistic case, the baryon mass"
316  << endl <<
317  "computation cannot be performed by the base class Star !"
318  << endl ;
319  abort() ;
320  }
321 
322  return *p_mass_b ;
323 
324 }
325 
326  //----------------------------//
327  // Gravitational mass //
328  //----------------------------//
329 
330 double Star::mass_g() const {
331 
332  if (p_mass_g == 0x0) { // a new computation is required
333 
334  cout <<
335  "Star::mass_g : in the relativistic case, the gravitational mass"
336  << endl <<
337  "computation cannot be performed by the base class Star !"
338  << endl ;
339  abort() ;
340  }
341 
342  return *p_mass_g ;
343 
344 }
345 
346 }
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:512
double * p_mass_b
Baryon mass.
Definition: star.h:268
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
double * p_ray_eq_3pis2
Coordinate radius at , .
Definition: star.h:251
double * p_ray_eq_pi
Coordinate radius at , .
Definition: star.h:248
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
Basic integer array class.
Definition: itbl.h:122
Scalar ent
Log-enthalpy.
Definition: star.h:190
double * p_ray_eq
Coordinate radius at , .
Definition: star.h:242
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinate ...
Definition: star_global.C:92
Tbl * p_xi_surf
Description of the stellar surface: 2-D Tbl containing the values of the radial coordinate on the su...
Definition: star.h:266
double * p_ray_eq_pis2
Coordinate radius at , .
Definition: star.h:245
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
double ray_eq_3pis2() const
Coordinate radius at , [r_unit].
Definition: star_global.C:236
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Definition: star_global.C:141
double * p_ray_pole
Coordinate radius at .
Definition: star.h:254
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l on...
Definition: star_global.C:66
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:189
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:281
Multi-domain grid.
Definition: grilles.h:279
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
virtual double mass_b() const =0
Baryon mass.
Definition: star_global.C:310
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
Itbl * p_l_surf
Description of the stellar surface: 2-D Itbl containing the values of the domain index l on the surfa...
Definition: star.h:260
double * p_mass_g
Gravitational mass.
Definition: star.h:269
virtual double mass_g() const =0
Gravitational mass.
Definition: star_global.C:330
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607