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 char star_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_global.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $" ;
28 
29 /*
30  * $Id: star_global.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $
31  * $Log: star_global.C,v $
32  * Revision 1.6 2014/10/13 08:53:39 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.5 2013/04/25 15:46:06 j_novak
36  * Added special treatment in the case np = 1, for type_p = NONSYM.
37  *
38  * Revision 1.4 2009/10/26 10:54:33 j_novak
39  * Added the case of a NONSYM base in theta.
40  *
41  * Revision 1.3 2007/06/21 19:55:09 k_taniguchi
42  * Introduction of a method to compute ray_eq_3pis2().
43  *
44  * Revision 1.2 2004/01/20 15:20:48 f_limousin
45  * First version
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Star/star_global.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $
49  *
50  */
51 
52 // Headers C
53 #include "math.h"
54 
55 // Headers Lorene
56 #include "star.h"
57 
58  //--------------------------//
59  // Stellar surface //
60  //--------------------------//
61 
62 namespace Lorene {
63 const Itbl& Star::l_surf() const {
64 
65  if (p_l_surf == 0x0) { // a new computation is required
66 
67  assert(p_xi_surf == 0x0) ; // consistency check
68 
69  int np = mp.get_mg()->get_np(0) ;
70  int nt = mp.get_mg()->get_nt(0) ;
71 
72  p_l_surf = new Itbl(np, nt) ;
73  p_xi_surf = new Tbl(np, nt) ;
74 
75  double ent0 = 0 ; // definition of the surface
76  double precis = 1.e-15 ;
77  int nitermax = 100 ;
78  int niter ;
79 
80  (ent.get_spectral_va()).equipot(ent0, nzet, precis, nitermax, niter, *p_l_surf,
81  *p_xi_surf) ;
82 
83  }
84 
85  return *p_l_surf ;
86 
87 }
88 
89 const Tbl& Star::xi_surf() const {
90 
91  if (p_xi_surf == 0x0) { // a new computation is required
92 
93  assert(p_l_surf == 0x0) ; // consistency check
94 
95  l_surf() ; // the computation is done by l_surf()
96 
97  }
98 
99  return *p_xi_surf ;
100 
101 }
102 
103 
104  //--------------------------//
105  // Coordinate radii //
106  //--------------------------//
107 
108 double Star::ray_eq() const {
109 
110  if (p_ray_eq == 0x0) { // a new computation is required
111 
112  const Mg3d& mg = *(mp.get_mg()) ;
113 
114  int type_t = mg.get_type_t() ;
115 #ifndef NDEBUG
116  int type_p = mg.get_type_p() ;
117 #endif
118  int nt = mg.get_nt(0) ;
119 
120  assert( (type_t == SYM) || (type_t == NONSYM) ) ;
121  assert( (type_p == SYM) || (type_p == NONSYM) ) ;
122  int k = 0 ;
123  int j = (type_t == SYM ? nt-1 : nt / 2);
124  int l = l_surf()(k, j) ;
125  double xi = xi_surf()(k, j) ;
126  double theta = M_PI / 2 ;
127  double phi = 0 ;
128 
129  p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
130 
131  }
132 
133  return *p_ray_eq ;
134 
135 }
136 
137 
138 double Star::ray_eq_pis2() const {
139 
140  if (p_ray_eq_pis2 == 0x0) { // a new computation is required
141 
142  const Mg3d& mg = *(mp.get_mg()) ;
143 
144  int type_t = mg.get_type_t() ;
145  int type_p = mg.get_type_p() ;
146  int nt = mg.get_nt(0) ;
147  int np = mg.get_np(0) ;
148 
149  int j = (type_t == SYM ? nt-1 : nt / 2);
150  double theta = M_PI / 2 ;
151  double phi = M_PI / 2 ;
152 
153  switch (type_p) {
154 
155  case SYM : {
156  int k = np / 2 ;
157  int l = l_surf()(k, j) ;
158  double xi = xi_surf()(k, j) ;
159  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
160  break ;
161  }
162 
163  case NONSYM : {
164  assert( (np == 1) || (np % 4 == 0) ) ;
165  int k = np / 4 ;
166  int l = l_surf()(k, j) ;
167  double xi = xi_surf()(k, j) ;
168  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
169  break ;
170  }
171 
172  default : {
173  cout << "Star::ray_eq_pis2 : the case type_p = "
174  << type_p << " is not contemplated yet !" << endl ;
175  abort() ;
176  }
177  }
178 
179  }
180 
181  return *p_ray_eq_pis2 ;
182 
183 }
184 
185 
186 double Star::ray_eq_pi() const {
187 
188  if (p_ray_eq_pi == 0x0) { // a new computation is required
189 
190  const Mg3d& mg = *(mp.get_mg()) ;
191 
192  int type_t = mg.get_type_t() ;
193  int type_p = mg.get_type_p() ;
194  int nt = mg.get_nt(0) ;
195  int np = mg.get_np(0) ;
196 
197  assert ( ( type_t == SYM ) || ( type_t == NONSYM ) ) ;
198 
199  switch (type_p) {
200 
201  case SYM : {
202  p_ray_eq_pi = new double( ray_eq() ) ;
203  break ;
204  }
205 
206  case NONSYM : {
207  int k = np / 2 ;
208  int j = (type_t == SYM ? nt-1 : nt/2 ) ;
209  int l = l_surf()(k, j) ;
210  double xi = xi_surf()(k, j) ;
211  double theta = M_PI / 2 ;
212  double phi = M_PI ;
213 
214  p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
215  break ;
216  }
217 
218  default : {
219 
220  cout << "Star::ray_eq_pi : the case type_p = " << type_p << endl ;
221  cout << " is not contemplated yet !" << endl ;
222  abort() ;
223  break ;
224  }
225  }
226 
227  }
228 
229  return *p_ray_eq_pi ;
230 
231 }
232 
233 double Star::ray_eq_3pis2() const {
234 
235  if (p_ray_eq_3pis2 == 0x0) { // a new computation is required
236 
237  const Mg3d& mg = *(mp.get_mg()) ;
238 
239  int type_t = mg.get_type_t() ;
240  int type_p = mg.get_type_p() ;
241  int nt = mg.get_nt(0) ;
242  int np = mg.get_np(0) ;
243 
244  assert( ( type_t == SYM ) || ( type_t == NONSYM ) ) ;
245 
246  int j = (type_t == SYM ? nt-1 : nt/2 );
247  double theta = M_PI / 2 ;
248  double phi = 3. * M_PI / 2 ;
249 
250  switch (type_p) {
251 
252  case SYM : {
253  p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
254  break ;
255  }
256 
257  case NONSYM : {
258  assert( np % 4 == 0 ) ;
259  int k = 3 * np / 4 ;
260  int l = l_surf()(k, j) ;
261  double xi = xi_surf()(k, j) ;
262  p_ray_eq_3pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
263  break ;
264  }
265 
266  default : {
267  cout << "Star::ray_eq_3pis2 : the case type_p = "
268  << type_p << " is not implemented yet !" << endl ;
269  abort() ;
270  }
271  }
272  }
273 
274  return *p_ray_eq_3pis2 ;
275 
276 }
277 
278 double Star::ray_pole() const {
279 
280  if (p_ray_pole == 0x0) { // a new computation is required
281 
282 #ifndef NDEBUG
283  const Mg3d& mg = *(mp.get_mg()) ;
284  int type_t = mg.get_type_t() ;
285 #endif
286  assert( (type_t == SYM) || (type_t == NONSYM) ) ;
287 
288  int k = 0 ;
289  int j = 0 ;
290  int l = l_surf()(k, j) ;
291  double xi = xi_surf()(k, j) ;
292  double theta = 0 ;
293  double phi = 0 ;
294 
295  p_ray_pole = new double( mp.val_r(l, xi, theta, phi) ) ;
296 
297  }
298 
299  return *p_ray_pole ;
300 
301 }
302 
303  //--------------------------//
304  // Baryon mass //
305  //--------------------------//
306 
307 double Star::mass_b() const {
308 
309  if (p_mass_b == 0x0) { // a new computation is required
310 
311  cout <<
312  "Star::mass_b : in the relativistic case, the baryon mass"
313  << endl <<
314  "computation cannot be performed by the base class Star !"
315  << endl ;
316  abort() ;
317  }
318 
319  return *p_mass_b ;
320 
321 }
322 
323  //----------------------------//
324  // Gravitational mass //
325  //----------------------------//
326 
327 double Star::mass_g() const {
328 
329  if (p_mass_g == 0x0) { // a new computation is required
330 
331  cout <<
332  "Star::mass_g : in the relativistic case, the gravitational mass"
333  << endl <<
334  "computation cannot be performed by the base class Star !"
335  << endl ;
336  abort() ;
337  }
338 
339  return *p_mass_g ;
340 
341 }
342 
343 }
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:495
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:462
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:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
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:89
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:233
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Definition: star_global.C:138
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:63
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:186
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
Multi-domain grid.
Definition: grilles.h:273
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
virtual double mass_b() const =0
Baryon mass.
Definition: star_global.C:307
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
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:327
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601