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