LORENE
etoile_global.C
1 /*
2  * Methods of class Etoile to compute global quantities
3  */
4 
5 /*
6  * Copyright (c) 2000-2001 Eric Gourgoulhon
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: etoile_global.C,v 1.8 2016/12/05 16:17:54 j_novak Exp $
31  * $Log: etoile_global.C,v $
32  * Revision 1.8 2016/12/05 16:17:54 j_novak
33  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34  *
35  * Revision 1.7 2014/10/13 08:52:59 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.6 2012/08/12 17:48:36 p_cerda
39  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
40  *
41  * Revision 1.5 2005/01/18 22:37:30 k_taniguchi
42  * Modify the method of ray_eq(int kk).
43  *
44  * Revision 1.4 2005/01/18 20:35:46 k_taniguchi
45  * Addition of ray_eq(int kk).
46  *
47  * Revision 1.3 2005/01/17 20:40:56 k_taniguchi
48  * Addition of ray_eq_3pis2().
49  *
50  * Revision 1.2 2003/12/05 14:50:26 j_novak
51  * To suppress some warnings...
52  *
53  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
54  * LORENE
55  *
56  * Revision 1.2 2000/07/21 12:02:14 eric
57  * Etoile::ray_eq_pi() : traitement du cas type_p = SYM.
58  *
59  * Revision 1.1 2000/01/28 17:18:45 eric
60  * Initial revision
61  *
62  *
63  * $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_global.C,v 1.8 2016/12/05 16:17:54 j_novak Exp $
64  *
65  */
66 
67 // Headers C
68 #include "math.h"
69 
70 // Headers Lorene
71 #include "etoile.h"
72 
73  //--------------------------//
74  // Stellar surface //
75  //--------------------------//
76 
77 namespace Lorene {
78 const Itbl& Etoile::l_surf() const {
79 
80  if (p_l_surf == 0x0) { // a new computation is required
81 
82  assert(p_xi_surf == 0x0) ; // consistency check
83 
84  int np = mp.get_mg()->get_np(0) ;
85  int nt = mp.get_mg()->get_nt(0) ;
86 
87  p_l_surf = new Itbl(np, nt) ;
88  p_xi_surf = new Tbl(np, nt) ;
89 
90  double ent0 = 0 ; // definition of the surface
91  double precis = 1.e-15 ;
92  int nitermax = 100 ;
93  int niter ;
94 
95  (ent().va).equipot(ent0, nzet, precis, nitermax, niter, *p_l_surf,
96  *p_xi_surf) ;
97 
98  }
99 
100  return *p_l_surf ;
101 
102 }
103 
104 const Tbl& Etoile::xi_surf() const {
105 
106  if (p_xi_surf == 0x0) { // a new computation is required
107 
108  assert(p_l_surf == 0x0) ; // consistency check
109 
110  l_surf() ; // the computation is done by l_surf()
111 
112  }
113 
114  return *p_xi_surf ;
115 
116 }
117 
118 
119  //--------------------------//
120  // Coordinate radii //
121  //--------------------------//
122 
123 double Etoile::ray_eq() const {
124 
125  if (p_ray_eq == 0x0) { // a new computation is required
126 
127  const Mg3d& mg = *(mp.get_mg()) ;
128 
129  int type_t = mg.get_type_t() ;
130 #ifndef NDEBUG
131  int type_p = mg.get_type_p() ;
132 #endif
133  int nt = mg.get_nt(0) ;
134 
135  if ( type_t == SYM ) {
136  assert( (type_p == SYM) || (type_p == NONSYM) ) ;
137  int k = 0 ;
138  int j = nt-1 ;
139  int l = l_surf()(k, j) ;
140  double xi = xi_surf()(k, j) ;
141  double theta = M_PI / 2 ;
142  double phi = 0 ;
143 
144  p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
145 
146  }
147  else {
148 
149  assert( (type_p == SYM) || (type_p == NONSYM) ) ;
150  int k = 0 ;
151  int j = (nt-1)/2 ;
152  int l = l_surf()(k, j) ;
153  double xi = xi_surf()(k, j) ;
154  double theta = M_PI / 2 ;
155  double phi = 0 ;
156 
157  p_ray_eq = new double( mp.val_r(l, xi, theta, phi) ) ;
158 
159 
160  // cout << "Etoile::ray_eq : the case type_t = " << type_t
161  // << " is not contemplated yet !" << endl ;
162  //abort() ;
163  }
164 
165  }
166 
167  return *p_ray_eq ;
168 
169 }
170 
171 
172 double Etoile::ray_eq_pis2() const {
173 
174  if (p_ray_eq_pis2 == 0x0) { // a new computation is required
175 
176  const Mg3d& mg = *(mp.get_mg()) ;
177 
178  int type_t = mg.get_type_t() ;
179  int type_p = mg.get_type_p() ;
180  int nt = mg.get_nt(0) ;
181  int np = mg.get_np(0) ;
182 
183  if ( type_t == SYM ) {
184 
185  int j = nt-1 ;
186  double theta = M_PI / 2 ;
187  double phi = M_PI / 2 ;
188 
189  switch (type_p) {
190 
191  case SYM : {
192  int k = np / 2 ;
193  int l = l_surf()(k, j) ;
194  double xi = xi_surf()(k, j) ;
195  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
196  break ;
197  }
198 
199  case NONSYM : {
200  assert( np % 4 == 0 ) ;
201  int k = np / 4 ;
202  int l = l_surf()(k, j) ;
203  double xi = xi_surf()(k, j) ;
204  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
205  break ;
206  }
207 
208  default : {
209  cout << "Etoile::ray_eq_pis2 : the case type_p = "
210  << type_p << " is not contemplated yet !" << endl ;
211  abort() ;
212  }
213  }
214 
215  }
216  else {
217 
218  int j = (nt-1)/2 ;
219  double theta = M_PI / 2 ;
220  double phi = M_PI / 2 ;
221 
222  switch (type_p) {
223 
224  case SYM : {
225  int k = np / 2 ;
226  int l = l_surf()(k, j) ;
227  double xi = xi_surf()(k, j) ;
228  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
229  break ;
230  }
231 
232  case NONSYM : {
233  assert( np % 4 == 0 ) ;
234  int k = np / 4 ;
235  int l = l_surf()(k, j) ;
236  double xi = xi_surf()(k, j) ;
237  p_ray_eq_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
238  break ;
239  }
240 
241  default : {
242  cout << "Etoile::ray_eq_pis2 : the case type_p = "
243  << type_p << " is not contemplated yet !" << endl ;
244  abort() ;
245  }
246  }
247 
248 
249 
250  }
251 
252  }
253 
254  return *p_ray_eq_pis2 ;
255 
256 }
257 
258 
259 double Etoile::ray_eq_pi() const {
260 
261  if (p_ray_eq_pi == 0x0) { // a new computation is required
262 
263  const Mg3d& mg = *(mp.get_mg()) ;
264 
265  int type_t = mg.get_type_t() ;
266  int type_p = mg.get_type_p() ;
267  int nt = mg.get_nt(0) ;
268  int np = mg.get_np(0) ;
269 
270  if ( type_t == SYM ) {
271 
272  switch (type_p) {
273 
274  case SYM : {
275  p_ray_eq_pi = new double( ray_eq() ) ;
276  break ;
277  }
278 
279  case NONSYM : {
280  int k = np / 2 ;
281  int j = nt-1 ;
282  int l = l_surf()(k, j) ;
283  double xi = xi_surf()(k, j) ;
284  double theta = M_PI / 2 ;
285  double phi = M_PI ;
286 
287  p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
288  break ;
289  }
290 
291  default : {
292 
293  cout << "Etoile::ray_eq_pi : the case type_t = " << type_t
294  << " and type_p = " << type_p << endl ;
295  cout << " is not contemplated yet !" << endl ;
296  abort() ;
297  break ;
298  }
299  }
300  }else{
301  switch (type_p) {
302 
303  case SYM : {
304  p_ray_eq_pi = new double( ray_eq() ) ;
305  break ;
306  }
307 
308  case NONSYM : {
309  int k = np / 2 ;
310  int j = (nt-1)/2 ;
311  int l = l_surf()(k, j) ;
312  double xi = xi_surf()(k, j) ;
313  double theta = M_PI / 2 ;
314  double phi = M_PI ;
315 
316  p_ray_eq_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
317  break ;
318  }
319 
320  default : {
321 
322  cout << "Etoile::ray_eq_pi : the case type_t = " << type_t
323  << " and type_p = " << type_p << endl ;
324  cout << " is not contemplated yet !" << endl ;
325  abort() ;
326  break ;
327  }
328  }
329 
330  }
331 
332  }
333 
334  return *p_ray_eq_pi ;
335 
336 }
337 
338 double Etoile::ray_eq_3pis2() const {
339 
340  if (p_ray_eq_3pis2 == 0x0) { // a new computation is required
341 
342  const Mg3d& mg = *(mp.get_mg()) ;
343 
344  int type_t = mg.get_type_t() ;
345  int type_p = mg.get_type_p() ;
346  int nt = mg.get_nt(0) ;
347  int np = mg.get_np(0) ;
348 
349  if ( type_t == SYM ) {
350 
351  int j = nt-1 ;
352  double theta = M_PI / 2 ;
353  double phi = 3. * M_PI / 2 ;
354 
355  switch (type_p) {
356 
357  case SYM : {
358  p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
359  break ;
360  }
361 
362  case NONSYM : {
363  assert( np % 4 == 0 ) ;
364  int k = 3 * np / 4 ;
365  int l = l_surf()(k, j) ;
366  double xi = xi_surf()(k, j) ;
367  p_ray_eq_3pis2 = new double( mp.val_r(l,xi,theta,phi) ) ;
368  break ;
369  }
370 
371  default : {
372  cout << "Etoile::ray_eq_3pis2 : the case type_p = "
373  << type_p << " is not contemplated yet !" << endl ;
374  abort() ;
375  }
376  }
377 
378  }
379  else {
380 
381  int j = (nt-1)/2 ;
382  double theta = M_PI / 2 ;
383  double phi = 3. * M_PI / 2 ;
384 
385  switch (type_p) {
386 
387  case SYM : {
388  p_ray_eq_3pis2 = new double( ray_eq_pis2() ) ;
389  break ;
390  }
391 
392  case NONSYM : {
393  assert( np % 4 == 0 ) ;
394  int k = 3 * np / 4 ;
395  int l = l_surf()(k, j) ;
396  double xi = xi_surf()(k, j) ;
397  p_ray_eq_3pis2 = new double( mp.val_r(l,xi,theta,phi) ) ;
398  break ;
399  }
400 
401  default : {
402  cout << "Etoile::ray_eq_3pis2 : the case type_p = "
403  << type_p << " is not contemplated yet !" << endl ;
404  abort() ;
405  }
406  }
407 
408 
409 
410  }
411 
412  }
413 
414  return *p_ray_eq_3pis2 ;
415 
416 }
417 
418 double Etoile::ray_pole() const {
419 
420  if (p_ray_pole == 0x0) { // a new computation is required
421 
422 #ifndef NDEBUG
423  const Mg3d& mg = *(mp.get_mg()) ;
424  int type_t = mg.get_type_t() ;
425 #endif
426  assert( (type_t == SYM) || (type_t == NONSYM) ) ;
427 
428  int k = 0 ;
429  int j = 0 ;
430  int l = l_surf()(k, j) ;
431  double xi = xi_surf()(k, j) ;
432  double theta = 0 ;
433  double phi = 0 ;
434 
435  p_ray_pole = new double( mp.val_r(l, xi, theta, phi) ) ;
436 
437  }
438 
439  return *p_ray_pole ;
440 
441 }
442 
443 double Etoile::ray_eq(int kk) const {
444 
445  const Mg3d& mg = *(mp.get_mg()) ;
446 
447  int type_t = mg.get_type_t() ;
448  int type_p = mg.get_type_p() ;
449  int nt = mg.get_nt(0) ;
450  int np = mg.get_np(0) ;
451 
452  assert( kk >= 0 ) ;
453  assert( kk < np ) ;
454 
455  double ray_eq_kk ;
456  if ( type_t == SYM ) {
457 
458  int j = nt-1 ;
459  double theta = M_PI / 2 ;
460 
461  switch (type_p) {
462 
463  case SYM : {
464  cout << "Etoile::ray_eq(kk) : the case type_p = "
465  << type_p << " is not contemplated yet !" << endl ;
466  abort() ;
467  }
468 
469  case NONSYM : {
470  double phi = 2. * kk * M_PI / np ;
471  int l = l_surf()(kk, j) ;
472  double xi = xi_surf()(kk, j) ;
473  ray_eq_kk = mp.val_r(l,xi,theta,phi) ;
474  break ;
475  }
476 
477  default : {
478  cout << "Etoile::ray_eq(kk) : the case type_p = "
479  << type_p << " is not contemplated yet !" << endl ;
480  abort() ;
481  }
482  }
483 
484  }
485  else {
486 
487  int j = (nt-1)/2 ;
488  double theta = M_PI / 2 ;
489 
490  switch (type_p) {
491 
492  case SYM : {
493  cout << "Etoile::ray_eq(kk) : the case type_p = "
494  << type_p << " is not contemplated yet !" << endl ;
495  abort() ;
496  }
497 
498  case NONSYM : {
499  double phi = 2. * kk * M_PI / np ;
500  int l = l_surf()(kk, j) ;
501  double xi = xi_surf()(kk, j) ;
502  ray_eq_kk = mp.val_r(l,xi,theta,phi) ;
503  break ;
504  }
505 
506  default : {
507  cout << "Etoile::ray_eq(kk) : the case type_p = "
508  << type_p << " is not contemplated yet !" << endl ;
509  abort() ;
510  }
511  }
512 
513 
514 
515 
516 
517 
518  }
519 
520  return ray_eq_kk ;
521 }
522 
523 
524  //--------------------------//
525  // Baryon mass //
526  //--------------------------//
527 
528 double Etoile::mass_b() const {
529 
530  if (p_mass_b == 0x0) { // a new computation is required
531 
532  if (relativistic) {
533  cout <<
534  "Etoile::mass_b : in the relativistic case, the baryon mass"
535  << endl <<
536  "computation cannot be performed by the base class Etoile !"
537  << endl ;
538  abort() ;
539  }
540 
541  assert(nbar.get_etat() == ETATQCQ) ;
542  p_mass_b = new double( nbar().integrale() ) ;
543 
544  }
545 
546  return *p_mass_b ;
547 
548 }
549 
550  //----------------------------//
551  // Gravitational mass //
552  //----------------------------//
553 
554 double Etoile::mass_g() const {
555 
556  if (p_mass_g == 0x0) { // a new computation is required
557 
558  if (relativistic) {
559  cout <<
560  "Etoile::mass_g : in the relativistic case, the gravitational mass"
561  << endl <<
562  "computation cannot be performed by the base class Etoile !"
563  << endl ;
564  abort() ;
565  }
566 
567  p_mass_g = new double( mass_b() ) ; // in the Newtonian case
568  // M_g = M_b
569 
570  }
571 
572  return *p_mass_g ;
573 
574 }
575 
576 }
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:512
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinate...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
virtual double mass_g() const
Gravitational mass.
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
double ray_eq() const
Coordinate radius at , [r_unit].
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
double * p_ray_pole
Coordinate radius at .
Definition: etoile.h:536
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.
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:462
virtual double mass_b() const
Baryon mass.
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
Tbl * p_xi_surf
Description of the stellar surface: 2-D Tbl containing the values of the radial coordinate on the su...
Definition: etoile.h:548
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
double * p_ray_eq
Coordinate radius at , .
Definition: etoile.h:524
double * p_mass_g
Gravitational mass.
Definition: etoile.h:551
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:435
double * p_ray_eq_pis2
Coordinate radius at , .
Definition: etoile.h:527
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:440
Multi-domain grid.
Definition: grilles.h:279
double * p_ray_eq_3pis2
Coordinate radius at , .
Definition: etoile.h:533
double ray_pole() const
Coordinate radius at [r_unit].
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:460
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l o...
Definition: etoile_global.C:78
double * p_ray_eq_pi
Coordinate radius at , .
Definition: etoile.h:530
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
double ray_eq_3pis2() const
Coordinate radius at , [r_unit].
double * p_mass_b
Baryon mass.
Definition: etoile.h:550
Itbl * p_l_surf
Description of the stellar surface: 2-D Itbl containing the values of the domain index l on the surfa...
Definition: etoile.h:542