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
27char 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
74namespace Lorene {
75const 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
101const 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
120double 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
169double 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
256double 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
335double 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
415double 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
440double 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
525double 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
551double 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}
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double * p_mass_b
Baryon mass.
Definition etoile.h:547
double * p_mass_g
Gravitational mass.
Definition etoile.h:548
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
double ray_eq_pi() const
Coordinate radius at , [r_unit].
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
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
Map & mp
Mapping associated with the star.
Definition etoile.h:429
double * p_ray_eq_3pis2
Coordinate radius at , .
Definition etoile.h:530
double * p_ray_eq_pi
Coordinate radius at , .
Definition etoile.h:527
virtual double mass_b() const
Baryon mass.
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
double * p_ray_pole
Coordinate radius at .
Definition etoile.h:533
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l ...
double ray_eq_3pis2() const
Coordinate radius at , [r_unit].
virtual double mass_g() const
Gravitational mass.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
double * p_ray_eq_pis2
Coordinate radius at , .
Definition etoile.h:524
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinat...
double * p_ray_eq
Coordinate radius at , .
Definition etoile.h:521
double ray_pole() const
Coordinate radius at [r_unit].
Basic integer array class.
Definition itbl.h:122
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
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
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Basic array class.
Definition tbl.h:161
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Lorene prototypes.
Definition app_hor.h:64